2.2.2. 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].

#### 2.2.3. 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–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<sup>2</sup> = 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.

**Figure 2.** Correlation of SOS between MCD12Q2 and EVI. Black line denotes linear regression lines. The DOY denotes the day of year.

#### *2.3. 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–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:

$$SOS\_{ural} = f\_{ural}(T\_{ural}) \tag{1}$$

$$\text{SOS}\_{\text{urban}} = f\_{\text{urban}}(T\_{\text{urban}}) \tag{2}$$

$$T\_{\text{contribution}} = f\_{\text{ural}}(T\_{\text{ural}}) - f\_{\text{ural}}(T\_{\text{urban}}) \tag{3}$$

$$Other\_{\text{contribution}} = f\_{\text{ural}}(T\_{\text{urban}}) - f\_{\text{urban}}(T\_{\text{urban}}) \tag{4}$$

$$T\_{percent} = T\_{contribution} / (T\_{contribution} + Other\_{contribution}) \tag{5}$$

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 *Tcontribution* and *Othercontribution* denote the contributions of the ΔLST and other factors to the ΔSOS, respectively; and the *Tpercent* 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). (3) However, in the same year, the average LST and SOS of the sample point in the urban was located at point D. There was a phenological difference C-D Equation (4) from point C, which refers to the influence of other environmental factors in the urban except temperature. (4) Furthermore, the percentage of the contribution of the ΔLST to the ΔSOS was calculated by Equation (5).

**Figure 3.** Temperature contribution separation model. The red/blue solid lines denote the linear regression line of the urban/rural sample point data. Points A and D denote the data of sample points in the rural and the urban in the same year, points B and C denote the predicted values of the model, and gray dotted lines denote the contributions of different factors. The DOY denotes the day of year.
