Quantification of Urban Heat Island-Induced Contribution to Advance in Spring Phenology: A Case Study in Hangzhou, China

Plant phenology is one of the key regulators of ecosystem processes, which are sensitive to environmental change. The acceleration of urbanization in recent years has produced substantial impacts on vegetation phenology over urban areas, such as the local warming induced by the urban heat island effect. However, quantitative contributions of the difference of land surface temperature (LST) between urban and rural (ΔLST) and other factors to the difference of spring phenology (i.e., the start of growing season, SOS) between urban and rural (ΔSOS) were rarely reported. Therefore, the objective of this study is to explore impacts of urbanization on SOS and distinguish corresponding contributions. Using Hangzhou, a typical subtropical metropolis, as the study area, vegetation index-based phenology data (MCD12Q2 and MYD13Q1 EVI) and land surface temperature data (MYD11A2 LST) from 2006–2018 were adopted to analyze the urban–rural gradient in phenology characteristics through buffers. Furthermore, we exploratively quantified the contributions of the ΔLST to the ΔSOS based on a temperature contribution separation model. We found that there was a negative coupling between SOS and LST in over 90% of the vegetated areas in Hangzhou. At the sample-point scale, SOS was weakly, but significantly, negatively correlated with LST at the daytime (R2 = 0.2 and p < 0.01 in rural; R2 = 0.14 and p < 0.05 in urban) rather than that at nighttime. Besides, the ΔSOS dominated by the ΔLST contributed more than 70% of the total ΔSOS. We hope this study could help to deepen the understanding of responses of urban ecosystem to intensive human activities.


Introduction
Plant phenology is the time of a certain growth event in the growth cycle, such as germination, branching, leafing, flowering, fruiting, defoliation and dormancy [1][2][3][4]. It directly or indirectly regulates several processes of plant growth, such as carbon and water cycle, playing a crucial role in the earth system [5,6]. Adapting to seasonal changes of the environment, plants show a growth rhythm, which is sensitive to environmental change [7,8]. As one of the most critical factors affecting plant phenology, an increase in temperature can promote the activity of enzymes, thereby prolonging plant development. Specifically, an increase in spring temperature promotes the release of plant dormancy in spring, and generally extends the growth cycle of plants [9][10][11][12][13][14].
Urbanization is an important feature of world development today, and it is one of the main causes of global environmental change in the 21st century. The acceleration of urbanization in recent years has produced substantial impacts on plant phenology over both urban areas and their rural surroundings [15][16][17][18][19]. This is mostly associated to the local warming effect induced by the urban heat island effect, which resulted from the increase in impervious surface percentage and anthropogenic heat emissions [20][21][22][23]. Moreover, it is as well as through the fertilization effect induced by the increase in the concentration of carbon dioxide (CO 2 ), nitrogen oxides (NO x ), and other atmospheric trace gases over urban surfaces [24][25][26]. These changes affect urban environments that plants depend on, and have impacts on the growth of plants, thereby changing the plant phenology [27,28].
At present, many studies have paid attention to impacts of urbanization on the change of plant phenology [27][28][29][30][31][32]. There are two methods to explore the impacts above: the historical comparison method and the urban-rural comparison method. The historical comparison method compared the phenology before and after urbanization, which was mainly for fast-developing cities [31]. However, due to the difficulty of obtaining long time series data, the historical comparison method is greatly restricted. The urban-rural comparison method used the data of the urban and the rural at the same time to explore the impact of urbanization on phenology, which is a method of changing space for time. The second one has been widely used, because of the great advantages in large-scale observations of remote sensing data [15][16][17][18][19]32]. Meng et al. investigated the urban and rural phenology of the of 85 giant cities in the continental United States from 2001 to 2014, and the results showed that the start of growing season (SOS) in the urban was 6 days earlier than that in the rural [33]. Wohlfahrt et al. found that with the acceleration of urbanization, the SOS advanced and the senescence delayed in the urban areas where the temperature rises [34]. Hu et al. used the Enhanced Vegetation Index (EVI) to explore the spatio-temporal changes of plant phenology and its response to land surface temperature (LST) in Northeast China [35]. The results showed that the LST was significantly negatively correlated with the SOS. Recently, most current studies focused on varieties of plant phenology and influences of temperature on plant phenology under urbanization, but did not quantitatively evaluate the contribution of the temperature differences to the phenological differences between the urban and the rural. That is, the quantitative contribution of the local warming induced by the urban heat island effect (the difference of LST between urban and rural, ∆LST) to the difference of spring phenology (SOS) between urban and rural (∆SOS) was less understood in past research.
With the development of statistical methods, it was possible to distinguish the influence of different factors. Li et al. used a statistical method to quantify the contribution of cooling and water supply to the yield benefits due to irrigation. They found that 16% of irrigation yield increase was due to irrigation cooling while the rest (84%) is due to water supply and other factors [36]. Besides, Zhao et al. also used a statistical method to quantifying the impacts of urbanization on vegetation growth. They found that the growth enhancement offset about 40% of direct loss of vegetation productivity caused by replacing productive vegetated surfaces with nonproductive impervious surfaces [16,30]. Based on the studies above, a statistical model was used to carry out this study.
Therefore, the objective of this study is to explore impacts of urbanization on SOS and exploratively distinguish contributions of local warming induced by the urban heat island effect (∆LST) and other factors to the difference of spring phenology between urban and rural (∆SOS). Hangzhou, a typical subtropical metropolis, was selected as the study area. Specifically, the spatial differences and inter-annual changes of the phenology in the urban and the rural were compared through a gradient analysis method using satellite-based phenology and LST data from 2006 to 2018. Then, the coupling relationship between phenology and temperature were investigated. After that, taking typical forest grid cells in the urban and the rural areas of Hangzhou as test samples, the local SOS was extracted using a remote sensing vegetation index from 2006 to 2018, and the difference of responses of SOS to LST between the urban and rural was explored. Finally, we exploratively distinguish quantitative contributions of the ∆LST and other factors to the ∆SOS.

Study Area
Hangzhou (118 • 21 -120 • 31 E, 29 • 11 -30 • 33 N) is the capital city of Zhejiang Province, whose GDP ranks among the top 10 in China, with 8.133 million in urban population, 2.227 million in rural population and an urbanization rate of 78.5% in 2019. It is located in the north of Zhejiang Province with a subtropical humid monsoon climate. As for temperature, it is lowest in January (average of 3-5 • C) and highest in July (average of 28-29 • C) with an annual average of 15.3-17 • C. The extreme maximum and minimum temperatures in Hangzhou reached 42.9 • C (31 July 1971) and −15 • C (5 January 1977). For precipitation, the annual average is 1100-1600 mm with rainy days of 130-160 days/year. There are two rainy seasons throughout the year. The first is the plum flood season from May to June, with an average rainfall of 350-500 mm, accounting for 25-31% of the year. The second rainy season is the typhoon rainy season from August to September, with an average rainfall of 120-220 mm, accounting for 8-13%. The forest coverage rate is over 64.77% (about 10,900 km 2 ), dominated by evergreen broad-leaved forests and deciduous broad-leaved forests. In this study, the land cover map with 10 m-spatial resolution in 2017 from Gong Peng Research Group of Tsinghua University was aggregated to pixels of 500 m × 500 m to extract forest areas for the following analysis. To assure both a certain level of homogeneity in the land cover type and an adequate number of pixels for a meaningful analysis, only the pixels of 500 m where the forest type was over 75% were included in this study. Besides, the multi-temporal dataset of global urban boundaries of 2018 was used to divide the scope of the urban and the rural of Hangzhou Figure 1. This dataset is derived from the Global Artificial Impervious Area-GAIA, released by Gong Peng Research Group of Tsinghua University [37]. Then, 5 test areas of deciduous broad-leaved forest were selected in the urban and the rural (i.e., 10 km away from the urban area) of Hangzhou, respectively. In each test area, 2 sample points, a total of 20 sample points, were extracted ( Figure 1). Besides, the Google map, latitude and longitude of the test areas in the urban and rural were showed in Table 1.           MODIS MYD11A2 LST product was used in this study, including the LST during the daytime and the nighttime, with a spatial resolution of 1000 m and a temporal resolution of 8 days. The LST of Aqua satellite is observed at 1:30 and 13:30 local solar time, the lowest and highest temperature of the day, which is more representative than that of Terra (monitored at 10:30 and 22:30) in the study of urban heat island. Therefore, the LST data of the Aqua satellite were used in this study [38,39]. The MODIS Reprojection Tool (MRT) was used to process the original images of LST, and so did the following data. Then, they were extracted in light of the study area and resampled to 500 m to be consistent with the phenology data. In addition, to explore the quantitative contributions of the ΔLST to the ΔSOS under urbanization, we collected 1000-m spatial resolution LST data of daytime and nighttime from 2006 to 2018 according to the coordinates of forest sample points.

MODIS Phenology
This paper tried to explore the difference of plant phenology between the urban and the rural of Hangzhou and its temporal and spatial patterns. The SOS from the MCD12Q2 (i.e., a MODIS phenology dataset) across the study area during 2006-2018 was used in this study. The spatial and time resolution of the dataset is 500 m and 1 year, respectively. The phenological events are derived from time series of MODIS 2-band Enhanced Vegetation Index (EVI2), which are fitted by QA/QC-weighted penalized cubic smoothing spline. The phenology data in some high-latitude regions and some semi-arid and arid environments exhibiting low-amplitude EVI2 variation, having uncertainty, while the data of Hangzhou in mid-latitude regions is relatively stable [40]. Besides, these data have been validated with field observations [41] and have been widely used and approved in previous studies [17,18,42].

Enhanced Vegetation Index
The contributions of the ΔLST to the ΔSOS were used to explore at a finer spatial resolution. Here, the Aqua MODIS 16-day EVI data (MYD13Q1, 250 m × 250 m), which matched the LST data from the same satellite (Aqua) to reduce uncertainties, was used to extract the phenology. Previous studies have suggested that EVI could accurately reflect MODIS MYD11A2 LST product was used in this study, including the LST during the daytime and the nighttime, with a spatial resolution of 1000 m and a temporal resolution of 8 days. The LST of Aqua satellite is observed at 1:30 and 13:30 local solar time, the lowest and highest temperature of the day, which is more representative than that of Terra (monitored at 10: 30 and 22:30) in the study of urban heat island. Therefore, the LST data of the Aqua satellite were used in this study [38,39]. The MODIS Reprojection Tool (MRT) was used to process the original images of LST, and so did the following data. Then, they were extracted in light of the study area and resampled to 500 m to be consistent with the phenology data. In addition, to explore the quantitative contributions of the ΔLST to the ΔSOS under urbanization, we collected 1000-m spatial resolution LST data of daytime and nighttime from 2006 to 2018 according to the coordinates of forest sample points.

MODIS Phenology
This paper tried to explore the difference of plant phenology between the urban and the rural of Hangzhou and its temporal and spatial patterns. The SOS from the MCD12Q2 (i.e., a MODIS phenology dataset) across the study area during 2006-2018 was used in this study. The spatial and time resolution of the dataset is 500 m and 1 year, respectively. The phenological events are derived from time series of MODIS 2-band Enhanced Vegetation Index (EVI2), which are fitted by QA/QC-weighted penalized cubic smoothing spline. The phenology data in some high-latitude regions and some semi-arid and arid environments exhibiting low-amplitude EVI2 variation, having uncertainty, while the data of Hangzhou in mid-latitude regions is relatively stable [40]. Besides, these data have been validated with field observations [41] and have been widely used and approved in previous studies [17,18,42].

Enhanced Vegetation Index
The contributions of the ΔLST to the ΔSOS were used to explore at a finer spatial resolution. Here, the Aqua MODIS 16-day EVI data (MYD13Q1, 250 m × 250 m), which matched the LST data from the same satellite (Aqua) to reduce uncertainties, was used to extract the phenology. Previous studies have suggested that EVI could accurately reflect MODIS MYD11A2 LST product was used in this study, including the LST during the daytime and the nighttime, with a spatial resolution of 1000 m and a temporal resolution of 8 days. The LST of Aqua satellite is observed at 1:30 and 13:30 local solar time, the lowest and highest temperature of the day, which is more representative than that of Terra (monitored at 10: 30 and 22:30) in the study of urban heat island. Therefore, the LST data of the Aqua satellite were used in this study [38,39]. The MODIS Reprojection Tool (MRT) was used to process the original images of LST, and so did the following data. Then, they were extracted in light of the study area and resampled to 500 m to be consistent with the phenology data. In addition, to explore the quantitative contributions of the ΔLST to the ΔSOS under urbanization, we collected 1000-m spatial resolution LST data of daytime and nighttime from 2006 to 2018 according to the coordinates of forest sample points.

MODIS Phenology
This paper tried to explore the difference of plant phenology between the urban and the rural of Hangzhou and its temporal and spatial patterns. The SOS from the MCD12Q2 (i.e., a MODIS phenology dataset) across the study area during 2006-2018 was used in this study. The spatial and time resolution of the dataset is 500 m and 1 year, respectively. The phenological events are derived from time series of MODIS 2-band Enhanced Vegetation Index (EVI2), which are fitted by QA/QC-weighted penalized cubic smoothing spline. The phenology data in some high-latitude regions and some semi-arid and arid environments exhibiting low-amplitude EVI2 variation, having uncertainty, while the data of Hangzhou in mid-latitude regions is relatively stable [40]. Besides, these data have been validated with field observations [41] and have been widely used and approved in previous studies [17,18,42].

Enhanced Vegetation Index
The contributions of the ΔLST to the ΔSOS were used to explore at a finer spatial resolution. Here, the Aqua MODIS 16-day EVI data (MYD13Q1, 250 m × 250 m), which matched the LST data from the same satellite (Aqua) to reduce uncertainties, was used to extract the phenology. Previous studies have suggested that EVI could accurately reflect MODIS MYD11A2 LST product was used in this study, including the LST during the daytime and the nighttime, with a spatial resolution of 1000 m and a temporal resolution of 8 days. The LST of Aqua satellite is observed at 1:30 and 13:30 local solar time, the lowest and highest temperature of the day, which is more representative than that of Terra (monitored at 10: 30 and 22:30) in the study of urban heat island. Therefore, the LST data of the Aqua satellite were used in this study [38,39]. The MODIS Reprojection Tool (MRT) was used to process the original images of LST, and so did the following data. Then, they were extracted in light of the study area and resampled to 500 m to be consistent with the phenology data. In addition, to explore the quantitative contributions of the ΔLST to the ΔSOS under urbanization, we collected 1000-m spatial resolution LST data of daytime and nighttime from 2006 to 2018 according to the coordinates of forest sample points.

MODIS Phenology
This paper tried to explore the difference of plant phenology between the urban and the rural of Hangzhou and its temporal and spatial patterns. The SOS from the MCD12Q2 (i.e., a MODIS phenology dataset) across the study area during 2006-2018 was used in this study. The spatial and time resolution of the dataset is 500 m and 1 year, respectively. The phenological events are derived from time series of MODIS 2-band Enhanced Vegetation Index (EVI2), which are fitted by QA/QC-weighted penalized cubic smoothing spline. The phenology data in some high-latitude regions and some semi-arid and arid environments exhibiting low-amplitude EVI2 variation, having uncertainty, while the data of Hangzhou in mid-latitude regions is relatively stable [40]. Besides, these data have been validated with field observations [41] and have been widely used and approved in previous studies [17,18,42].

Enhanced Vegetation Index
The contributions of the ΔLST to the ΔSOS were used to explore at a finer spatial resolution. Here, the Aqua MODIS 16-day EVI data (MYD13Q1, 250 m × 250 m), which matched the LST data from the same satellite (Aqua) to reduce uncertainties, was used to extract the phenology. Previous studies have suggested that EVI could accurately reflect MODIS MYD11A2 LST product was used in this study, including the LST during the daytime and the nighttime, with a spatial resolution of 1000 m and a temporal resolution of 8 days. The LST of Aqua satellite is observed at 1:30 and 13:30 local solar time, the lowest and highest temperature of the day, which is more representative than that of Terra (monitored at 10:30 and 22:30) in the study of urban heat island. Therefore, the LST data of the Aqua satellite were used in this study [38,39]. The MODIS Reprojection Tool (MRT) was used to process the original images of LST, and so did the following data. Then, they were extracted in light of the study area and resampled to 500 m to be consistent with the phenology data. In addition, to explore the quantitative contributions of the ΔLST to the ΔSOS under urbanization, we collected 1000-m spatial resolution LST data of daytime and nighttime from 2006 to 2018 according to the coordinates of forest sample points.

MODIS Phenology
This paper tried to explore the difference of plant phenology between the urban and the rural of Hangzhou and its temporal and spatial patterns. The SOS from the MCD12Q2 (i.e., a MODIS phenology dataset) across the study area during 2006-2018 was used in this study. The spatial and time resolution of the dataset is 500 m and 1 year, respectively. The phenological events are derived from time series of MODIS 2-band Enhanced Vegetation Index (EVI2), which are fitted by QA/QC-weighted penalized cubic smoothing spline. The phenology data in some high-latitude regions and some semi-arid and arid environments exhibiting low-amplitude EVI2 variation, having uncertainty, while the data of Hangzhou in mid-latitude regions is relatively stable [40]. Besides, these data have been validated with field observations [41] and have been widely used and approved in previous studies [17,18,42].

Enhanced Vegetation Index
The contributions of the ΔLST to the ΔSOS were used to explore at a finer spatial resolution. Here, the Aqua MODIS 16-day EVI data (MYD13Q1, 250 m × 250 m), which matched the LST data from the same satellite (Aqua) to reduce uncertainties, was used to extract the phenology. Previous studies have suggested that EVI could accurately reflect     MODIS MYD11A2 LST product was used in this study, including the LST during the daytime and the nighttime, with a spatial resolution of 1000 m and a temporal resolution of 8 days. The LST of Aqua satellite is observed at 1:30 and 13:30 local solar time, the lowest and highest temperature of the day, which is more representative than that of Terra (monitored at 10:30 and 22:30) in the study of urban heat island. Therefore, the LST data of the Aqua satellite were used in this study [38,39]. The MODIS Reprojection Tool (MRT) was used to process the original images of LST, and so did the following data. Then, they were extracted in light of the study area and resampled to 500 m to be consistent with the phenology data. In addition, to explore the quantitative contributions of the ∆LST to the ∆SOS under urbanization, we collected 1000-m spatial resolution LST data of daytime and nighttime from 2006 to 2018 according to the coordinates of forest sample points.

MODIS Phenology
This paper tried to explore the difference of plant phenology between the urban and the rural of Hangzhou and its temporal and spatial patterns. The SOS from the MCD12Q2 (i.e., a MODIS phenology dataset) across the study area during 2006-2018 was used in this study. The spatial and time resolution of the dataset is 500 m and 1 year, respectively. The phenological events are derived from time series of MODIS 2-band Enhanced Vegetation Index (EVI2), which are fitted by QA/QC-weighted penalized cubic smoothing spline. The phenology data in some high-latitude regions and some semi-arid and arid environments exhibiting low-amplitude EVI2 variation, having uncertainty, while the data of Hangzhou in mid-latitude regions is relatively stable [40]. Besides, these data have been validated with field observations [41] and have been widely used and approved in previous studies [17,18,42].

Enhanced Vegetation Index
The contributions of the ∆LST to the ∆SOS were used to explore at a finer spatial resolution. Here, the Aqua MODIS 16-day EVI data (MYD13Q1, 250 m × 250 m), which matched the LST data from the same satellite (Aqua) to reduce uncertainties, was used to extract the phenology. Previous studies have suggested that EVI could accurately reflect the growth status of vegetation, and effectively extract phenology at both regional and local scales [43][44][45][46][47]. In this study, the asymmetric Gaussian function, which is widely used in curve fitting and phenological extraction, while a 20% threshold was used to fit the EVI curve and extract the SOS for each selected forest sample point from 2006 to 2018 [6,7].
Notably, both the MODIS phenology product and EVI-derived phenology were used in this study. On the one hand, to explore the changes in plant phenology caused by the urban heat island effect, we needed to focus on the forest, which is severely affected by urbanization. At this time, the plant phenology data with a resolution of 500 m cannot satisfy the demand, so the EVI data with a resolution of 250 m at a smaller scale was used to improve the resolution and reduce the influence of mixed pixels. On the other hand, MCD12Q2 uses EVI2 data to extract plant phenology, while EVI2 data lacks the blue band, different from EVI [48]. Therefore, the more widely used and robust EVI data was used for extracting more accurate plant phenology. Considering the reasons above, we selected EVI-based SOS for the research of sample-scale.
Besides, the SOS extracted from MCD12Q2 phenology data and EVI data showed a linear correlation (R 2 = 0.54) at the forest sample area in Hangzhou Figure 2. First, the original data and methods to extract SOS were different. The one used the EVI data and the method of asymmetric Gaussian function with a 20% threshold, while the other used the EVI2 data and the method of QA/QC-weighted penalized cubic smoothing spline. Although the correlation of SOS derived from MCD12Q2 and EVI was not satisfied, the relationship was statistically significant (p < 0.05), which could support the consistency between them. Second, in this study, we aimed to explore the quantitative contributions of the difference of land surface temperature between urban and rural and other factors to the difference of spring phenology (∆SOS) under urbanization. The relative ∆SOS was effective instead of the absolute value of SOS. Therefore, despite the difference in the absolute SOS of the two data, they had a statistically significant linear relationship and there was also a certain relationship between the ∆SOS. Besides, the RMSE of the two data was 5 days, which was smaller than the average ∆SOS of >9 days. In summary, the two data of SOS had a certain consistency in this study.
inal data and methods to extract SOS were different. The one used the EVI data and the method of asymmetric Gaussian function with a 20% threshold, while the other used the EVI2 data and the method of QA/QC-weighted penalized cubic smoothing spline. Although the correlation of SOS derived from MCD12Q2 and EVI was not satisfied, the relationship was statistically significant (p < 0.05), which could support the consistency between them. Second, in this study, we aimed to explore the quantitative contributions of the difference of land surface temperature between urban and rural and other factors to the difference of spring phenology (ΔSOS) under urbanization. The relative ΔSOS was effective instead of the absolute value of SOS. Therefore, despite the difference in the absolute SOS of the two data, they had a statistically significant linear relationship and there was also a certain relationship between the ΔSOS. Besides, the RMSE of the two data was 5 days, which was smaller than the average ΔSOS of >9 days. In summary, the two data of SOS had a certain consistency in this study.

Temperature Contribution Separation Model
The rapid development of urbanization dramatically changes the environments which terrestrial ecosystem depended on. Compared with the rural surroundings, there are differences in temperature, photoperiod and atmosphere conditions, having a certain impact on plant phenology. A large number of studies have shown that the acceleration of urbanization in recent years produced substantial impacts on plant phenology over

Temperature Contribution Separation Model
The rapid development of urbanization dramatically changes the environments which terrestrial ecosystem depended on. Compared with the rural surroundings, there are differences in temperature, photoperiod and atmosphere conditions, having a certain impact on plant phenology. A large number of studies have shown that the acceleration of urbanization in recent years produced substantial impacts on plant phenology over both urban areas and their rural surroundings [15][16][17][18]33]. Therefore, in order to distinguish the contributions of ∆LST and other factors between urban and rural to the difference of spring phenology (∆SOS), we followed the statistical method of quantifying the contributions of cooling and water supply to the yield benefits due to irrigation of Li et al. [36], establishing a temperature contribution separation model based on the laboratory of the rural and urban of Hangzhou.
Firstly, we performed regression analysis on SOS and average LST during the daytime in spring (February, March, April) from 2006 to 2018 in the rural and the urban Equations (1) and (2), respectively. Secondly, Equations (3)-(5) were used to distinguish the contributions of the ∆LST and other factors to the ∆SOS: Other T percent = T contribution /(T contribution + Other contribution ) where the subscripts of rural and urban denote the corresponding parameters of the rural and the urban, respectively; the SOS denotes the start of the growing season; the T denotes the average LST during the daytime in spring; the f denotes the regression relationship between LST and SOS. The T contribution and Other contribution denote the contributions of the ∆LST and other factors to the ∆SOS, respectively; and the T percent denotes the percentage of the contribution of the ∆LST.
The temperature contribution separation model was shown in Figure 3. (1) In the figure, the blue and red lines represent the regression relationships between SOS and LST in the rural and the urban Equations (1) and (2), respectively. Points A and D represent the average LST and SOS of the sample points in the same year of the rural and urban, respectively. Line B-D is the difference of the SOS between the rural and the urban.
(2) We supposed that the rural surroundings were heated to reach the LST of the urban in the same year. Then, the SOS of the rural (point A) in that year moved to point C according to Equation (3). At this time, the two points A and C were the phenological state only in different LST, and line B-C was the phenological difference only when the LST rose Equation (3), which refers to the influence of the urban heat island effect (∆LST).

Statistical Analysis
In this study, the multi-year average phenology of the SOS of 2006-2018 were calculated to compare the differences of phenology between that in the urban and the rural of Hangzhou, which were extracted from the MCD12Q2 phenology dataset. Besides, a buffer analysis method was adopted to compare the difference of SOS in urban-rural gradient and its relationship with LST in more detail [7,28]. First, the urban boundary of Hangzhou in 2018 was derived from the global urban boundary dataset [12]. Second, the circular buffer zones outside the urban boundary were drawn every 2 km, which were 2, 4, 6, 8, 10, 12, 14, 16, 18 and 20 km, respectively. Finally, the average of SOS and LST in each buffer zone were calculated to explore the relationship between SOS and LST with the distance away from the urban.
Furthermore, the two-factor combination mapping is a very intuitive visualization method that can express the coupling relationship between two variables. The two factors that are incomparable numerically can be compared in a hierarchical manner, and the different levels of the two factors are matched in pairs to different combinations, representing the different coupling relationships of the two factors. In this study, to explore the coupling relationship between SOS and LST, the natural breakpoint classification method was used to divide the daily average LST and SOS in spring of 2018 into three levels: low, medium and high. Additionally, we used the two-factor combination mapping method to display the different coupling relationships, including low LST-low SOS, low LST-medium SOS, low LST-high SOS, medium LST-low SOS, medium LST-medium SOS, medium LST-high SOS, high LST-low SOS, high LST-medium SOS, and high LST-high SOS.

Statistical Analysis
In this study, the multi-year average phenology of the SOS of 2006-2018 were calculated to compare the differences of phenology between that in the urban and the rural of Hangzhou, which were extracted from the MCD12Q2 phenology dataset. Besides, a buffer analysis method was adopted to compare the difference of SOS in urban-rural gradient and its relationship with LST in more detail [7,28]. First, the urban boundary of Hangzhou in 2018 was derived from the global urban boundary dataset [12]. Second, the circular buffer zones outside the urban boundary were drawn every 2 km, which were 2, 4, 6, 8, 10, 12, 14, 16, 18 and 20 km, respectively. Finally, the average of SOS and LST in each buffer zone were calculated to explore the relationship between SOS and LST with the distance away from the urban.
Furthermore, the two-factor combination mapping is a very intuitive visualization method that can express the coupling relationship between two variables. The two factors that are incomparable numerically can be compared in a hierarchical manner, and the different levels of the two factors are matched in pairs to different combinations, representing the different coupling relationships of the two factors. In this study, to explore the coupling relationship between SOS and LST, the natural breakpoint classification method was used to divide the daily average LST and SOS in spring of 2018 into three levels: low, medium and high. Additionally, we used the two-factor combination mapping method to display the different coupling relationships, including low LST-low SOS, low LST-medium SOS, low LST-high SOS, medium LST-low SOS, medium LST-medium SOS, medium LST-high SOS, high LST-low SOS, high LST-medium SOS, and high LST-high SOS.

Spatial Patterns of the LST and SOS in Hangzhou
The average LST in spring of 2018 was utilized to explore the difference of LST between that in the urban and the rural Figure 4. The results showed that there was a significant spatial heterogeneity in the LST of Hangzhou, showing a gradient of high in the urban and low in the rural. In terms of spatial distribution, the area with a LST greater than 23.5 • C accounted for 10.9%, mainly distributed in the urban area; 16.6%, 32.9% and 26.6% of the area with a LST of 21.5-23.5, 19.5-21.5 and 17.5-19.5 • C distributed in the middle area of Hangzhou, respectively. The area with a LST less than 17.5 • C (accounting for 13%) was mainly distributed in the northern and southern edges of Hangzhou. Moreover, as shown in the inset chart in Figure 4, the LST followed a generally decreasing urban-rural gradient, showing a significant urban heat island effect that the LST was highest (24 • C) in the urban and lowest (18.6 • C) in the rural. In the range of 0-6 km from the urban, the LST decreased fast (0.72 • C/km), and the decrease tended to slow down (0.08 • C/km) after 6 km.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 19 for 13%) was mainly distributed in the northern and southern edges of Hangzhou. Moreover, as shown in the inset chart in Figure 4, the LST followed a generally decreasing urban-rural gradient, showing a significant urban heat island effect that the LST was highest (24 °C) in the urban and lowest (18.6 °C) in the rural. In the range of 0-6 km from the urban, the LST decreased fast (0.72 °C/km), and the decrease tended to slow down (0.08 °C/km) after 6 km. As shown in Figure 5, the annual average SOS of Hangzhou from 2006 to 2018 showed a significant spatial heterogeneity. The SOS was earlier in the urban and later in the rural. In terms of spatial distribution, the area with the SOS less than 76 day of year (termed DOY) accounted for 8.9%, mainly distributed in the urban and the area within 2 km from the urban, located in the east and south of Hangzhou. About 23.1%, 31.1% and 26.6% of the area with SOS of 76-83, 83-88 and 88-94 DOY were distributed in the middle area of Hangzhou, respectively. The area with SOS more than 94 DOY (accounting for 10.3%) mainly distributed in the northern edges of Hangzhou, in the mountainous areas with higher elevations. Besides, as shown in the inset chart in Figure 5, the SOS followed a generally increasing urban-rural gradient, that is, from urban (79 DOY) to rural (87 DOY), the SOS was continuously delayed by 8 days. In the range of 0-6 km from the urban, the SOS increased fast (1.02 days/km), and the increase tended to slow down (0.14 As shown in Figure 5, the annual average SOS of Hangzhou from 2006 to 2018 showed a significant spatial heterogeneity. The SOS was earlier in the urban and later in the rural. In terms of spatial distribution, the area with the SOS less than 76 day of year (termed DOY) accounted for 8.9%, mainly distributed in the urban and the area within 2 km from the urban, located in the east and south of Hangzhou. About 23.1%, 31.1% and 26.6% of the area with SOS of 76-83, 83-88 and 88-94 DOY were distributed in the middle area of Hangzhou, respectively. The area with SOS more than 94 DOY (accounting for 10.3%) mainly distributed in the northern edges of Hangzhou, in the mountainous areas with higher elevations. Besides, as shown in the inset chart in Figure 5, the SOS followed a generally increasing urban-rural gradient, that is, from urban (79 DOY) to rural (87 DOY), the SOS was continuously delayed by 8 days. In the range of 0-6 km from the urban, the SOS increased fast (1.02 days/km), and the increase tended to slow down (0.14 days/km) after 6 km. We further explored the spatial difference of SOS between that in the urban and the rural of Hangzhou from 2006 to 2018, and we found that the results from each year had little significant fluctuations. Therefore, in order to avoid information redundancy and excessively long images, we displayed the results every 4 years (2006, 2010, 2014 and 2018) Figure 6. The results showed that although the absolute value of SOS varied from year to year, the spatial differentiation of SOS yearly was consistent with the annual average SOS from 2006 to 2018 in Figure 5. They both showed a significant spatial heterogeneity that the SOS was smaller in the eastern and southern area of Hangzhou and larger in the northern marginal area. As shown in the inset chart in Figure 6, the SOS of the urban was 9, 9, 6 and 6 days earlier than the rural in 2006, 2010, 2014 and 2018, respectively. In addition, the SOS followed a generally increasing urban-rural gradient. In 2006, 2010, 2014 and 2018, the SOS increased fast (1.25, 1.05, 0.83 and 0.93 days/km, respectively) within the range of 0-6 km from the urban, while it tended to be stable (0.07, 0.21, 0.05 and 0.04 days/km, respectively) after 6 km.
Combined with the analysis of the average LST in spring of 2018 and SOS across Hangzhou above, the difference in LST and SOS presented an opposite state and change trend. That is, the LST tended to be higher in the urban and lower in the rural, while SOS tended to be earlier in the urban and later in the rural. The LST followed a generally decreasing urban-rural gradient, while the opposite occurred for SOS, but both the LST and SOS varied greatly within the range of 0-6 km and then tended to be stable. In summary, the LST and SOS showed a negative correlation, and the coupling relationship between them would be further explored below. We further explored the spatial difference of SOS between that in the urban and the rural of Hangzhou from 2006 to 2018, and we found that the results from each year had little significant fluctuations. Therefore, in order to avoid information redundancy and excessively long images, we displayed the results every 4 years (2006, 2010, 2014 and 2018) Figure 6. The results showed that although the absolute value of SOS varied from year to year, the spatial differentiation of SOS yearly was consistent with the annual average SOS from 2006 to 2018 in Figure 5. They both showed a significant spatial heterogeneity that the SOS was smaller in the eastern and southern area of Hangzhou and larger in the northern marginal area. As shown in the inset chart in Figure 6, the SOS of the urban was 9, 9, 6 and 6 days earlier than the rural in 2006, 2010, 2014 and 2018, respectively. In addition, the SOS followed a generally increasing urban-rural gradient. In 2006, 2010, 2014 and 2018, the SOS increased fast (1.25, 1.05, 0.83 and 0.93 days/km, respectively) within the range of 0-6 km from the urban, while it tended to be stable (0.07, 0.21, 0.05 and 0.04 days/km, respectively) after 6 km.
Combined with the analysis of the average LST in spring of 2018 and SOS across Hangzhou above, the difference in LST and SOS presented an opposite state and change trend. That is, the LST tended to be higher in the urban and lower in the rural, while SOS tended to be earlier in the urban and later in the rural. The LST followed a generally decreasing urban-rural gradient, while the opposite occurred for SOS, but both the LST and SOS varied greatly within the range of 0-6 km and then tended to be stable. In summary, the LST and SOS showed a negative correlation, and the coupling relationship between them would be further explored below. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 19

Relationship between LST and SOS in Hangzhou
The results above showed that the spatial distribution and urban-rural gradient of the average LST in spring and SOS showed an opposite trend. Therefore, it could be inferred that there was negative correlation between LST and SOS. Figure 7 shows the coupling relationship between LST and SOS, and the 94.6% of the area conformed to the inference above. The LST and SOS showed a significant negative correlation accounted for 53.9%, more than a half of the total area. Among them, the low LST-high SOS (accounting for 6.2%) mainly distributed in the northern edge of Hangzhou, the outer rural farthest from urban. The medium LST-medium SOS (45.3%) was mainly distributed in the middle part of Hangzhou, with a moderate distance from urban. The high LST-low SOS (2.4%) was mainly distributed in urban and within 2 km from urban. At the same time, there was 40.7% of the area that LST and SOS showed a weaker negative correlation, including low LST-medium SOS (10.9%), medium LST-low SOS (12.2%), medium LST-high SOS (10.1%), and high LST-medium SOS (7.5%), mainly distributed in the area between urban and outer rural. Besides, there was 5.4% of the area that exhibited a contrary relationship to the inference. That is, the LST and SOS showed a significantly positive correlations, which were low LST-low SOS (3.3%) and high LST-high SOS (2.1%). It might be related to the threshold of LST and SOS for the classification. In addition, in areas with low, moderate and high LST, 83%, 67% and 82% of the SOS has a medium-high, moderate and medium-low distribution, respectively. Overall, the results above confirmed the inference that spring LST was negatively correlated with SOS.

Relationship between LST and SOS in Hangzhou
The results above showed that the spatial distribution and urban-rural gradient of the average LST in spring and SOS showed an opposite trend. Therefore, it could be inferred that there was negative correlation between LST and SOS. Figure 7 shows the coupling relationship between LST and SOS, and the 94.6% of the area conformed to the inference above. The LST and SOS showed a significant negative correlation accounted for 53.9%, more than a half of the total area. Among them, the low LST-high SOS (accounting for 6.2%) mainly distributed in the northern edge of Hangzhou, the outer rural farthest from urban. The medium LST-medium SOS (45.3%) was mainly distributed in the middle part of Hangzhou, with a moderate distance from urban. The high LST-low SOS (2.4%) was mainly distributed in urban and within 2 km from urban. At the same time, there was 40.7% of the area that LST and SOS showed a weaker negative correlation, including low LST-medium SOS (10.9%), medium LST-low SOS (12.2%), medium LST-high SOS (10.1%), and high LST-medium SOS (7.5%), mainly distributed in the area between urban and outer rural. Besides, there was 5.4% of the area that exhibited a contrary relationship to the inference. That is, the LST and SOS showed a significantly positive correlations, which were low LST-low SOS (3.3%) and high LST-high SOS (2.1%). It might be related to the threshold of LST and SOS for the classification. In addition, in areas with low, moderate and high LST, 83%, 67% and 82% of the SOS has a medium-high, moderate and medium-low distribution, respectively. Overall, the results above confirmed the inference that spring LST was negatively correlated with SOS. In order to further verify the inference above, the annual average SOS of 2006-2018 and its change trend at different levels of LST was calculated to explore the relationship between LST and SOS. As shown in Figure 8a, the annual average SOS at low, medium, and high LST were 88.9, 85.8, and 85.0 DOY, respectively. The SOS continently decreased with the increase of the LST. That is, the spring phenology continued to advance. As In order to further verify the inference above, the annual average SOS of 2006-2018 and its change trend at different levels of LST was calculated to explore the relationship between LST and SOS. As shown in Figure 8a, the annual average SOS at low, medium, and high LST were 88.9, 85.8, and 85.0 DOY, respectively. The SOS continently decreased with the increase of the LST. That is, the spring phenology continued to advance. As shown in Figure 8b, the change trend of SOS at low, medium, and high LST were −5.1, −3.9, and −2 days/10 years, respectively. The downward trend of SOS decreased with the increase of LST. That is, the rate of advancement of phenology continuously slowed down. In general, the inference established that there was a negative correlation between SOS and LST, and it developed in a consistent direction, showing a trend of convergence.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 19 shown in Figure 8b, the change trend of SOS at low, medium, and high LST were −5.1, −3.9, and −2 days/10 years, respectively. The downward trend of SOS decreased with the increase of LST. That is, the rate of advancement of phenology continuously slowed down. In general, the inference established that there was a negative correlation between SOS and LST, and it developed in a consistent direction, showing a trend of convergence.

Relationship between LST and SOS at Sample Points
The relationship between LST and SOS at the sample point scale was further employed Figure 9. The LST in the urban of Hangzhou was significantly higher than that in the rural, exhibiting a significant urban heat island effect. In terms of the LST during the daytime in spring, the difference between that in the urban (23.0 ± 1.2 °C) and the rural (19.7 ± 1.0 °C) was significant, reaching 3.3 ± 1.0 °C. From 2006 to 2018, the difference continued to increase, with an average annual increase of 0.2 °C (p < 0.01). For the LST during the nighttime in spring, the difference between that in the urban (8.2 ± 1.0 °C) and the rural (7.0 ± 0.7 °C) reached 1.2 ± 0.3 °C. Compared with the daytime, the urban heat island effect was weaker at night, reduced by 2.1 °C. Integrating the temperature during the daytime and nighttime Figure 9c, the difference of the daily average temperature between that in the urban (15.6 ± 0.8 °C) and the rural (13.3 ± 0.8 °C) reached 2.3 ± 0.5 °C, and it continued to increase from 2006 to 2018 (Slope = 0.09 °C/years, p < 0.01). To further compare the difference of SOS between that in the urban and the rural Figure 9d, the SOS of urban (79.3 ± 5.6 DOY) was 7.4 ± 2.7 days earlier than that in the rural (86.7 ± 5.0 DOY). It indicated that the plant phenology changed significantly under different environmental backgrounds in the urban and the rural. Besides, it was worth noting that, except 2012 (1.5 days) and 2015 (12.9 days), the difference of SOS between that in the urban and the rural was relatively stable from 2006 to 2018, with a difference of 7.5 ± 1.6 days, and it was consistent with the difference of 6-9 days in Hangzhou Figure 6.
As shown in Figure 10, the LST during the daytime rather than nighttime showed a statistically significantly negative correlation with SOS both in the urban and the rural. In terms of the LST during the daytime Figure 10a, the sample points of LST-SOS in the urban and the rural distributed significantly separately, with LST of 19-28 °C and 16-23 °C, and SOS of 40-100 DOY and 50-110 DOY in the urban and the rural, respectively. Besides, the SOS in the rural was statistically significantly negatively correlated with the LST (Slope = −1.64 days/°C; p < 0.01). The correlation in urban had a degree of significance (Slope = −1.01 days/°C, p <0.05), and the change trend of the SOS with LST increasing in the urban was smaller than that in the rural. For the LST during the nighttime Figure 10b, the sample points of LST-SOS in the rural and the urban had poor separation and high similarity. There was less statistically significant correlation between SOS and LST (p >

Relationship between LST and SOS at Sample Points
The relationship between LST and SOS at the sample point scale was further employed Figure 9. The LST in the urban of Hangzhou was significantly higher than that in the rural, exhibiting a significant urban heat island effect. In terms of the LST during the daytime in spring, the difference between that in the urban (23.0 ± 1.2 • C) and the rural (19.7 ± 1.0 • C) was significant, reaching 3.3 ± 1.0 • C. From 2006 to 2018, the difference continued to increase, with an average annual increase of 0.2 • C (p < 0.01). For the LST during the nighttime in spring, the difference between that in the urban (8.2 ± 1.0 • C) and the rural (7.0 ± 0.7 • C) reached 1.2 ± 0.3 • C. Compared with the daytime, the urban heat island effect was weaker at night, reduced by 2.1 • C. Integrating the temperature during the daytime and nighttime Figure 9c, the difference of the daily average temperature between that in the urban (15.6 ± 0.8 • C) and the rural (13.3 ± 0.8 • C) reached 2.3 ± 0.5 • C, and it continued to increase from 2006 to 2018 (Slope = 0.09 • C/years, p < 0.01). To further compare the difference of SOS between that in the urban and the rural Figure 9d, the SOS of urban (79.3 ± 5.6 DOY) was 7.4 ± 2.7 days earlier than that in the rural (86.7 ± 5.0 DOY). It indicated that the plant phenology changed significantly under different environmental backgrounds in the urban and the rural. Besides, it was worth noting that, except 2012 (1.5 days) and 2015 (12.9 days), the difference of SOS between that in the urban and the rural was relatively stable from 2006 to 2018, with a difference of 7.5 ± 1.6 days, and it was consistent with the difference of 6-9 days in Hangzhou Figure 6.
As shown in Figure 10, the LST during the daytime rather than nighttime showed a statistically significantly negative correlation with SOS both in the urban and the rural. In terms of the LST during the daytime Figure 10a, the sample points of LST-SOS in the urban and the rural distributed significantly separately, with LST of 19-28 • C and 16-23 • C, and SOS of 40-100 DOY and 50-110 DOY in the urban and the rural, respectively. Besides, the SOS in the rural was statistically significantly negatively correlated with the LST (Slope = −1.64 days/ • C; p < 0.01). The correlation in urban had a degree of significance (Slope = −1.01 days/ • C, p < 0.05), and the change trend of the SOS with LST increasing in the urban was smaller than that in the rural. For the LST during the nighttime Figure 10b, the sample points of LST-SOS in the rural and the urban had poor separation and high similarity. There was less statistically significant correlation between SOS and LST (p > 0.05), indicating that the LST during the nighttime had little effect on SOS. However, the fitted trends all showed that the SOS decreased with the increasing LST, which was consistent with the response of SOS to LST in Figure 8. 0.05), indicating that the LST during the nighttime had little effect on SOS. However, the fitted trends all showed that the SOS decreased with the increasing LST, which was consistent with the response of SOS to LST in Figure 8.   0.05), indicating that the LST during the nighttime had little effect on SOS. However, the fitted trends all showed that the SOS decreased with the increasing LST, which was consistent with the response of SOS to LST in Figure 8.

Relative Contributions of ∆LST to ∆SOS at Sample Points
The temperature contribution separation model was utilized to explore the contributions of the ∆LST to the ∆SOS under urbanization. The difference of SOS under urbanization between predicted by the model (7.3 ± 1.3 days) and observed through data (7.4 ± 2.7 days) was −0.16 ± 1.4 days Figure 11a. For the results predicted by the model Figure 11b, the ∆LST played a significant role in the advance of SOS under urbanization and the proportion of its contributions was relatively stable, despite the interannual fluctuations. Besides, we found that the ∆SOS dominated by the ∆LST contributed 72 ± 13.3% (5.3 ± 1.7 of 7.3 ± 1.3 days) to the difference of SOS between that in the urban and the rural Figure 11a. The advance of 28 ± 13.3% (1.9 ± 0.7 days) of SOS under urbanization was dominated by other factors such as photoperiod and air pollution. Overall, the local warming effect induced by the urban heat island effect produced substantial impacts on plant phenology under urbanization, but the impact of other factors also cannot be ignored.

Relative Contributions of ΔLST to ΔSOS at Sample Points
The temperature contribution separation model was utilized to explore the contributions of the ΔLST to the ΔSOS under urbanization. The difference of SOS under urbanization between predicted by the model (7.3 ± 1.3 days) and observed through data (7.4 ± 2.7 days) was −0.16 ± 1.4 days Figure 11a. For the results predicted by the model Figure 11b, the ΔLST played a significant role in the advance of SOS under urbanization and the proportion of its contributions was relatively stable, despite the interannual fluctuations. Besides, we found that the ΔSOS dominated by the ΔLST contributed 72 ± 13.3% (5.3 ± 1.7 of 7.3 ± 1.3 days) to the difference of SOS between that in the urban and the rural Figure 11a. The advance of 28 ± 13.3% (1.9 ± 0.7 days) of SOS under urbanization was dominated by other factors such as photoperiod and air pollution. Overall, the local warming effect induced by the urban heat island effect produced substantial impacts on plant phenology under urbanization, but the impact of other factors also cannot be ignored.

Discussion
Plant phenology is one of the key regulators of ecosystem processes, which is sensitive to environmental change. The acceleration of urbanization in recent years has produced substantial impacts on vegetation phenology over urban areas, such as the local warming induced by the urban heat island effect. This study explored impacts of urbanization on SOS and distinguished contributions of ΔLST and other factors to ΔSOS based on a temperature contribution separation model. We found that the SOS was negatively correlated with the daytime LST in spring, and the ΔSOS dominated by the ΔLST and other factors contributed 72% and 28% to the ΔSOS, respectively. Previous studies showed that there were lots of aspects besides temperature were different between urban and rural under urbanization, which had a certain impact on plant phenology. (1) The land cover changes under urbanization changed the soil properties extremely in the urban, affecting the relationship between plants, water and nutrients [49]. (2) There were high concentration of greenhouse gases (such as CO2) and major pollutants (such as NO, NO2, CO, SO2 and particulates with a diameter of 10μm or less) in the urban, resulting from the emissions produced by factories and automobiles [50]. In this regard, many studies showed that pollutants in the urban environment could cause the advance or delay of plant phenology [51][52][53]. (3) Due to the increasing artificial light caused by the human activities at night in the urban, the growth of plants was seriously influenced [54,55]. Therefore, there were various differences between the urban areas and their rural surroundings, not only in terms of temperature, but also in other aspects that affect plant phenology.

Discussion
Plant phenology is one of the key regulators of ecosystem processes, which is sensitive to environmental change. The acceleration of urbanization in recent years has produced substantial impacts on vegetation phenology over urban areas, such as the local warming induced by the urban heat island effect. This study explored impacts of urbanization on SOS and distinguished contributions of ∆LST and other factors to ∆SOS based on a temperature contribution separation model. We found that the SOS was negatively correlated with the daytime LST in spring, and the ∆SOS dominated by the ∆LST and other factors contributed 72% and 28% to the ∆SOS, respectively. Previous studies showed that there were lots of aspects besides temperature were different between urban and rural under urbanization, which had a certain impact on plant phenology. (1) The land cover changes under urbanization changed the soil properties extremely in the urban, affecting the relationship between plants, water and nutrients [49]. (2) There were high concentration of greenhouse gases (such as CO 2 ) and major pollutants (such as NO, NO 2 , CO, SO 2 and particulates with a diameter of 10µm or less) in the urban, resulting from the emissions produced by factories and automobiles [50]. In this regard, many studies showed that pollutants in the urban environment could cause the advance or delay of plant phenology [51][52][53]. (3) Due to the increasing artificial light caused by the human activities at night in the urban, the growth of plants was seriously influenced [54,55]. Therefore, there were various differences between the urban areas and their rural surroundings, not only in terms of temperature, but also in other aspects that affect plant phenology.
The LST in spring from February to April was selected to explore the relationship between LST and SOS in this study. Previous studies showed that meteorological parameters such as temperature and precipitation in a period of time before the phenological event were important determinants affecting the occurrence of phenology. The period of time is significant to the study of the relationship between phenology and climate [56,57], which called preseason duration. Polgar et al. found that the temperature in late winter and spring or preseason temperature played an important role in the occurrence of SOS [58]. Zhou et al. and Jia et al. found that the LST showed a statistically significant correlation with SOS (p < 0.05, R 2 > 0.8) [7,15]. In addition, different temperature indicators were used to explore its effect on the SOS, such as daily maximum temperature and diurnal temperature difference. (1) Piao et al. showed that the SOS was more sensitive to the preseason daily maximum temperature in the northern hemisphere, and 68% of the European Union and 83% of the United States had a preseason duration of 0-3 months [59]. (2) The results of Huang et al. showed that 77.2% of the northern hemisphere, the SOS had the strongest correlation with the average diurnal temperature difference in the preseason period of 1-3 months [60]. In general, the results above indicated that the temperature in spring was relevant to SOS with the preseason duration of 0-3 months.
At the same time, there were some limitations in this study. (1) As shown in Figure 10, the LST during the daytime showed a statistically significantly negative correlation with SOS both in the urban (p < 0.05) and the rural (p < 0.01), while the R 2 was low. The low R 2 might be related to the impurity of the data where existed many mixed pixels in the urban sample. To explore the changes in plant phenology caused by the urban heat island effect, we needed to focus on forest where is severely affected by urbanization. However, the MODIS EVI data with a resolution of 250 m was the data with long time series, the highest resolution and we could currently obtain. Therefore, the problem of mixed pixels inevitably existed in the urban, which affected the correlation between SOS and LST. Although the R 2 was relatively low, they were both statistically significant (p < 0.05), which was meaningful in a certain degree. In the future, data with higher resolutions should be used to reduce the uncertainty caused by the data and make results more reliable. (2) In this study, all types of forest in Hangzhou were used to explore the impact of urbanization on plant phenology. Previous studies found that the phenology and the response to urbanization varied in the different types of vegetation. However, due to the limitation of the accuracy of data, we could only exclude shrub, farmland and grassland. Relatively uniform and stable forests were extracted as the study object to weaken the impact of different vegetation types to a certain extent. Reliably, previous studies utilized all vegetation types for the study and got reliable results, having a certain significance [7,33]. In future research, we hope to more finely distinguish the vegetation types and improve the accuracy of the study, in order to obtain more reliable results. (3) As mentioned above, there were lots of aspects besides temperature that were different between urban and rural under urbanization, which had a certain impact on plant phenology. Although other factors were taken into consideration in this study, no specific data analysis was carried out on other factors such as greenhouse gases. In future studies, a more quantitative and detailed discussion on the effects of other factors on vegetation phenology should be advanced. With the development of high-quality data and online data processing platforms (e.g., Google Earth Engine), it is feasible to apply the methodology of our present work to other study areas even to a global scale.

Conclusions
Focused on Hangzhou, a typical subtropical metropolis, to investigate the impact of urbanization on plant phenology. Vegetation index-based phenology data and land surface temperature data were adopted to analyze the urban-rural gradient in phenology characteristics and the contributions of ∆LST to the ∆SOS under urbanization. We found that there was a negative coupling between SOS and LST in over 90% of the vegetated areas in Hangzhou. At the sample-point scale, SOS was weakly, but significantly, negatively correlated with LST at the daytime (R 2 = 0.2 and p < 0.01 in rural; R 2 = 0.14 and p < 0.05 in urban) rather than that at nighttime. Besides, the ∆SOS dominated by the ∆LST contributed more than 70% of the total ∆SOS. We consider that the achievements of this study will provide quantitative evidence for the impact of urbanization on the plant phenology, and help to deepen the understanding of urban ecosystem adaptation under intensive human activities.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The study did not report any data.