Characterizing Spatial Patterns of the Response Rate of Vegetation Green-Up Dates to Land Surface Temperature in Beijing, China (2001–2019)

: The phenology indicator of vegetation green-up dates (GUD) is prone to being affected by changes in temperature. However, the inﬂuencing degree of urbanization-induced temperature warming on vegetation GUDs among different vegetation species along the urban-rural gradient remains inadequately described. In this study, based on the long-term (2001–2019) satellite-derived vegetation GUDs and nighttime land surface temperature (LST) of forests, grasslands, and croplands along the urban-rural gradient with Beijing (China) as a case study area, the responses of vegetation GUDs to temperature changes were quantitatively analyzed, taking into account the vegetation types and distances away from the urban domain. The results show that (1) long-term GUDs and LST are signiﬁcantly negatively correlated, characterized by a weaker signiﬁcant correlation near the urban area when compared with its surrounding areas, with the greatest absolute linear correlation coefﬁcients ( r ) happening at rings 32 km ( r max = − 0.93, forests), 20 km and 48 km ( r max = − 0.83, grasslands), and 34 km ( r max = − 0.82, croplands), respectively; (2) the magnitude of change in GUDs over the past 19 year (2001–2019) are signiﬁcantly positively correlated with these in LST near the urban area, demonstrating a distance-decay trend, with the greatest advance in GUDs occurring at the ring nearest the urban area, by about 20 days (forests), 24.5 days (grasslands), and 15.6 days (croplands), respectively; (3) the spatial pattern of the response rate of GUDs change to LST change (days K − 1 ) also showed a declining trend with distance, with GUD advanced by 6.8 days K − 1 (forests), 7.5 days K − 1 (grasslands), and 4.9 days K − 1 (croplands) at the closest ring to the urban, decreasing to about 2.3 days K − 1 (48 km) , 4.1 days K − 1 (18 km), and 1 day K − 1 (18 km), respectively, indicating a notable inﬂuence of temperature warming on vegetation GUDs near the urban domains.

Temperature, especially the heat accumulation in the preceding period of time (preseason), has been long recognized as a main driver of plant growth [5,14,15], playing a decisive role in determining vegetation phenology [16,17].In the past few decades, China experienced extremely fast urbanization, and the trend is expected to continue in the coming decades [18][19][20][21], which has caused and will further lead to environment change.The accelerating urbanization has resulted in substantial impacts on the natural environment over both urban areas and their rural surroundings [4,[22][23][24][25][26].One of the resultant effects is featured with higher temperature in urban areas than the surrounding rural areas [27,28].In this context, cities that experienced a dramatic urbanization process are natural laboratories to investigate the responses of vegetation phenology to climate warming [3,18,[29][30][31].A number of recent studies have shown that vegetation GUDs and land surface temperature (LST) are negatively correlated [16,27], and vegetation green-up comes earlier in urban areas than that in rural areas [24,32].In addition, the study conducted by Zhang et al. [31] in the eastern North America reported that the differences in vegetation phenology in urban areas and surrounding regions are significantly correlated to temperature differences.In general, temperature rising would likely advance the green-up date and delay the end of season, leading to an extended growing season [2][3][4]27].However, vegetation phonological events differed among species and areas, and even between the same plant cover in different areas [3,33], and the urbanization-induced temperature changes also vary across different cities, making responses of vegetation phenology to changes in temperature across different areas diverse [1].Accordingly, case studies about what extent the impact of temperature warming on vegetation GUDs according to different urban-rural areas and species of vegetation is needed to better understand the responses of vegetation to temperature changes.
Remote sensing-based phenology has been significantly improved with the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard Terra and Aqua satellites [34,35].MODIS observations are now used routinely for building MODIS global products, such as the vegetation phenology product and land surface temperature (LST) product, providing useful data for regional-to-global phenology studies [24,34,36].Studies on urban-rural vegetation phenology using satellite-derived phenology and climate or environmental factors are gaining increasing interests [31,37,38].For example, Li et al. [3] characterized spatiotemporal dynamics in the phenology of urban ecosystems based on Landsat data; Ji et al. [22] quantified the urban heat island-induced variations in spring phenology using MODIS phenology and temperature products; Zhang et al. [2] explored the responses of four vegetation phenology indicators to urbanization amplitudes using the impervious surface area as an indicator.However, these studies were carried out either focusing on approach comparison [3], the changing trend of different vegetation phenology indicators over a long period [22,39], or on the relationships between vegetation phenology indicators and amplitude of urbanization [1,40].Vegetation species or temperature differences along urban-rural gradients were not considered in these studies.The species of vegetation is also an influencing factor affecting the phenological response to climate change [30].A significant difference in phenological shifts has been found between crops and other natural vegetation, such as grass and mixed forest [31].It is necessary to consider vegetation types in the investigation of vegetation responses to microclimate changes along urban-rural gradient.Furthermore, more studies are required to fully understand how temperature, in particularly remote sensing-based observations (e.g., MODIS Land Surface Temperature product (LST)), can be incorporated to better depict urban vegetation phenology.
In this study, we examined the effect of temperature changes on vegetation phenology along the urban-rural area of Beijing (China) using MODIS temperature and phenology products for the period 2001-2019.To quantify the dynamic responses of vegetation phenology to temperature changes, a series of buffer zones, i.e., concentric rings from the urban-dominated region to the surrounding rural areas, were created.The nighttime temperature and green-up dates of three main plant types, including forests, grasslands, and croplands, in different buffer zones, were computed.The main objectives of this study are to explore the relationships between green-up dates and temperature as well as the response rate of green-up dates to temperature changes over the past 19 years (2001-2019) along the urban-rural gradient, and to make some attempts to detail the spatial pattern of the effect of temperature warming on green-up dates to shed a little light on the understanding of future GUDs changes when temperature is supposed to be increasing.

Study Area
The study area is located in the suburban area around Beijing built-up lands.Beijing is in the north of China, which is latitudinally located between 39 • and 41 • .Beijing has a warm, temperate, semi-humid, and semi-arid monsoon climate with distinguished seasonal temperature differences, cold dry winter, with short spring and autumn.The average annual precipitation and temperature are around 500 mm and 11-12 • C, respectively [41].Over the past two decades, Beijing has experienced an accelerating urbanization process and noticeable expansion of its urban area [42,43].The huge built-up areas have dramatically affected LST in urban areas and its surrounding areas, featuring higher temperatures in central built-up areas than the surrounding suburban areas [44].Therefore, it is a good area to carry out the study of the response of vegetation phenology to temperature changes.The specific study area was set within an annular area extending 50 km outward from the perimeter of the built-up area of Beijing (Figure 1), including the central and southern areas of Beijing, the northwestern areas of Tianjin, and a small part of the Hebei province, which borders Beijing.
Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 18 phenology to temperature changes, a series of buffer zones, i.e., concentric rings from the urban-dominated region to the surrounding rural areas, were created.The nighttime temperature and green-up dates of three main plant types, including forests, grasslands, and croplands, in different buffer zones, were computed.The main objectives of this study are to explore the relationships between green-up dates and temperature as well as the response rate of green-up dates to temperature changes over the past 19 years (2001-2019) along the urban-rural gradient, and to make some attempts to detail the spatial pattern of the effect of temperature warming on green-up dates to shed a little light on the understanding of future GUDs changes when temperature is supposed to be increasing.

Study Area
The study area is located in the suburban area around Beijing built-up lands.Beijing is in the north of China, which is latitudinally located between 39° and 41°.Beijing has a warm, temperate, semi-humid, and semi-arid monsoon climate with distinguished seasonal temperature differences, cold dry winter, with short spring and autumn.The average annual precipitation and temperature are around 500 mm and 11-12 °C , respectively [41].Over the past two decades, Beijing has experienced an accelerating urbanization process and noticeable expansion of its urban area [42,43].The huge built-up areas have dramatically affected LST in urban areas and its surrounding areas, featuring higher temperatures in central built-up areas than the surrounding suburban areas [44].Therefore, it is a good area to carry out the study of the response of vegetation phenology to temperature changes.The specific study area was set within an annular area extending 50 km outward from the perimeter of the built-up area of Beijing (Figure 1), including the central and southern areas of Beijing, the northwestern areas of Tianjin, and a small part of the Hebei province, which borders Beijing.

Collection and Preprocessing of Remote Sensing Data
In this study, three remotely sensed observations were used, including vegetation green-up dates, land surface temperature, and land cover data.These data were derived from three remote sensing data products, i.e., MCD12Q2 (global vegetation phenology product), MCD12Q1 (the land cover product), and MOD11A2 (land surface temperature product), which are all produced by the Moderate Resolution Imaging Spectroradiometer

Collection and Preprocessing of Remote Sensing Data
In this study, three remotely sensed observations were used, including vegetation green-up dates, land surface temperature, and land cover data.These data were derived from three remote sensing data products, i.e., MCD12Q2 (global vegetation phenology product), MCD12Q1 (the land cover product), and MOD11A2 (land surface temperature product), which are all produced by the Moderate Resolution Imaging Spectroradiometer (MODIS) data (https://modis.gsfc.nasa.gov/,accessed on 1 January 2001).All used products are in version-6 and cover a period of 2001-2019.MCD12Q1 and MCD12Q2 are annual products, with a spatial resolution of 500 m, and MOD11A2 data are 8-day composites with a spatial resolution of 1 km.In order to be consistent with the spatial resolution of MCD12Q1 and MCD12Q2, MYD11A2 LST datasets were resampled to 500 m.
Resampling method is the nearest neighbor method that comes with the mask of update Mask function on GEE (Google Earth Engine) platform.That is, temperature products with a resolution of 1 km are masked with land use data of 500 m, and temperature data of 500 m can be obtained directly.As MCD12Q2 raw data are in the Julian day form, with 1 January 1970 as the start pointing, all MCD12Q2 data were first converted into the day of year (DOY).Since MCD12Q1 can be processed by five classification methods, the International Geosphere-Biosphere Programme (IGBP) classification was adopted and the obtained classification results were aggregated into four types, i.e., forests, grasslands, croplands, and urban.The data processing was mainly conducted through the GEE platform.Land use data in the study area for every year were obtained through screening of "QC" quality control bands.Land use data masks of corresponding years were used to obtain data of forests, grasslands, and croplands for this year for analysis.MOD11A2 LST data from the Terra satellite monitored at 22:30 (nighttime) local solar time, was used in this study.

Extraction of Vegetation Green-Up Dates (GUDs)
Phonological transitions like green-up dates (GUDs) can be identified from season dynamics via spectral vegetation indices [45] such as the enhanced vegetation index (EVI) and the normalized difference vegetation index (NDVI), which can be derived from various sensors of ground/satellite platforms, including the moderate-resolution imaging spectroradiometer (MODIS), the Landsat/TM data, and so on.Numerous studies have proved the reliable quality of vegetation phenology indicators derived from the MODIS phenology product [1,24,34,36,[46][47][48].In the present study, to identify the timing of vegetation greenup dates (GUDs) from 2001 to 2019, the Collection 6 MODIS Land Cover Dynamics product (MCD12Q2) (V6.0) was used.This product uses time series of the two band Enhanced Vegetation Index (EVI2) to define the timing of the total number of vegetation phenophase transition dates during each growing season in 500 m resolution.Compared with Collection 4 and Collection 5 products, the methodological approach in Collection 6 increased the reliability of retrieved phenometrics in tropical, arid, and semi-arid eco-systems, and can better capture phenometrics in systems with multiple vegetation cycles per year [49].Previous studies have reported that MCD12Q2 provides a reliable measure of vegetation phenology [36,46].MCD12Q2 (V6.0) contains several key phenology metrics, including the onset of greenness, green-up midpoint, maturity, peak greenness, senescence, green-down midpoint, dormancy, EVI2 minimum, EVI2 amplitude, integrated EVI2 over a vegetation cycle, as well as overall and phenology metric-specific quality information.For our analysis, the MCD12Q2 'green-up' metric, in which GUDs (Green-up Dates) is defined as the day of year (day of year, DOY) during the green-up phase in spring when the EVI2 time series at each pixel crosses 15% of its seasonal amplitude [49], was adopted for GUDs extraction.The annual vegetation green-up dates for forests, grasslands, and croplands during 2001-2019 were derived with each ring as the basic unit.

Methods of Correlation and Response Rate Analysis
To identify the spatial pattern of changes in GUDs and LST (Land Surface Temperature) along the urban-rural gradient, multiple concentric rings seeding upon the perimeter of an urban area with a step of 2 km and a total of 25 ring zones within 50 km beyond the urban area, were generated (Figure 2).Each buffer ring was used as a basic unit within which the following were computed: (1) the average values of GUDs and nighttime LST for the period 2001-2019; (2) the differences in GUDs and LST over the past 19 years, i.e., GUDs and LST in 2019 minus these in 2001, which is expressed as ∆GUDs and ∆LST, respectively; and (3) the absolute value of response rate of ∆GUDs to ∆LST, i.e., days advanced in GUDs with 1 K increase in LST (unit: days K −1 ).
for the period 2001-2019; (2) the differences in GUDs and LST over the past 19 years, i.e.GUDs and LST in 2019 minus these in 2001, which is expressed as ΔGUDs and ΔLST respectively; and (3) the absolute value of response rate of ΔGUDs to ΔLST, i.e., days advanced in GUDs with 1 K increase in LST (unit: days K −1 ).Steps for correlation and response rate analysis were shown in Figure 3.It is necessary to emphasize that in order to eliminate the effect of changes in land cover on calculations of vegetation GUDs and LST in different years, the MODIS products, including land cover product, land surface temperature product, and phenology product, were obtained for every year from 2001 to 2019, which means they are individual from each other and unaffected by land cover variations among different years.Since the GUDs are severely affected by temperature in spring, the springtime (March to April) observations of LST were adopted.The spatial pattern of the relationship between LST and vegetation GUDs, the relationships between temperature changes (ΔLST) and GUDs changes (ΔGUDs), and the response rate of ΔGUDs to ΔLST, were examined by general linear regression models, with the linear correlation coefficient (r) or determination coefficient (R 2 ) used as indicators.Steps for correlation and response rate analysis were shown in Figure 3.It is necessary to emphasize that in order to eliminate the effect of changes in land cover on calculations of vegetation GUDs and LST in different years, the MODIS products, including land cover product, land surface temperature product, and phenology product, were obtained for every year from 2001 to 2019, which means they are individual from each other, and unaffected by land cover variations among different years.Since the GUDs are severely affected by temperature in spring, the springtime (March to April) observations of LST were adopted.The spatial pattern of the relationship between LST and vegetation GUDs, the relationships between temperature changes (∆LST) and GUDs changes (∆GUDs), and the response rate of ∆GUDs to ∆LST, were examined by general linear regression models, with the linear correlation coefficient (r) or determination coefficient (R 2 ) used as indicators.

Spatial Patterns of GUDs and LST
The value of GUDs and LST of every year from 2001 to 2019 were computed and the spatial patterns were generated by 19-year average values of GUDs and LST (Figure 4).From the spatial patterns of the 19-year average GUDs and LST across the whole study area (Figure 3), it can be seen that vegetation in the southeast of the study area showed relatively earlier GUDs when compared with these in the northwest (Figure 4a), with the GUDs of most vegetation coming before DOY (Day of Year) 100, even before DOY 90 at some points.However, the corresponding average LST in this area were not significantly high (Figure 4c).Taking into account the vegetation types (Figure 1) in this area and the features of GUDs of the three investigated vegetation covers (Figure 4b), it can be inferred that sites that have GUDs earlier than DOY 90 were croplands, which have GUDs earlier than DOY100 (8.86% of croplands have GUDs earlier than DOY90 and 28.79% earlier than DOY100), in contrast to forests and grasslands.Meanwhile, it can be seen that GUDs in the northwest of the study area were mostly later than DOY110, with some sites later than DOY120 and a small amount of sites later than DOY130 (Figure 4a), corresponding to large areas of forests, some grasslands, and small areas of croplands (Figure 1).GUDs of forests were concentrated between DOY100 and DOY120 (about 97%), while those of grasslands and croplands were more dispersed, especially croplands, in which GUDs earlier than DOY 90 and later than DOY130 occupied about 8.86% and 9.35% respectively (Figure 4b).The phenomenon that the croplands in the south and southeast of the study area have GUDs smaller than DOY90 was mainly caused by winter wheat in these areas, and the greater GUDs (>DOY120) in the croplands in the northwest were mostly due the relatively lower temperature.On average, GUDs in croplands were earlier than forests, and forests earlier than grasslands.The earlier GUDs of croplands was in agreement with the results reported by Zhang et al. [1], who found that croplands have the earliest green-up date, followed by forests and grasslands.The 19-year average LST showed a decay pattern with increasing distances (Figure 4c), with the urban domains demonstrating an obvious higher temperature than the surrounding areas.The highest nighttime temperature happened in the central city with an LST of 282.7 K, and the lowest LST of 270.2 K was in the northwest, about 50 km away from the city central, with a temperature difference of 12.5 k, which fully displayed an urbanizationinduced warming.Besides, contrary to the dispersed distribution of GUDs in croplands, the smallest variations of LST in croplands were observed, with 89.5% croplands having an LST of 276-278 K (Figure 4d).LST in forests and grasslands showed a relative wide range, with 25% forests and 21% grasslands having LST lower than 275 K.
From spatial patterns of average GUDs and LST (Figure 4) as well as land cover types in the study area (Figure 1), it can be concluded that vegetation GUDs first and foremost were determined by vegetation types, highlighting the necessity to consider the vegetation types in the investigation of the responses of GUDs to temperature warming.This view was also supported by Li et al. [3] and Wang [33], who found distinct different phonological characteristics among different vegetation types and emphasized that the strength of climate modulations on phenology varies among vegetation types.
followed by forests and grasslands.The 19-year average LST showed a decay pattern with increasing distances (Figure 4c), with the urban domains demonstrating an obvious higher temperature than the surrounding areas.The highest nighttime temperature happened in the central city with an LST of 282.7 K, and the lowest LST of 270.2 K was in the northwest, about 50 km away from the city central, with a temperature difference of 12.5 k, which fully displayed an urbanization-induced warming.Besides, contrary to the dispersed distribution of GUDs in croplands, the smallest variations of LST in croplands were observed, with 89.5% croplands having an LST of 276-278 K (Figure 4d).LST in forests and grasslands showed a relative wide range, with 25% forests and 21% grasslands having LST lower than 275 K. From spatial patterns of average GUDs and LST (Figure 4) as well as land cover types in the study area (Figure 1), it can be concluded that vegetation GUDs first and foremost were determined by vegetation types, highlighting the necessity to consider the vegetation types in the investigation of the responses of GUDs to temperature warming.This view was also supported by Li et al. [3] and Wang [33], who found distinct different phonological characteristics among different vegetation types and emphasized that the strength of climate modulations on phenology varies among vegetation types.

Correlations between GUDs and LST
The quantitative relationships between average vegetation GUDs and LST over 2001-2019 (Figure 5) showed that vegetation GUDs were closely and negatively correlated with LST, either for forestlands, grasslands, or croplands, indicating that the rising in temperature tends to advance the vegetation GUDs, in agreement with results reported by previous studies [1,16,27,33,41].The coefficient of determination (R 2 ) of linear fitted models between GUDs and LST for forests, grasslands, and croplands were similar, around 0.74, 0.72, and 0.71 with an average value of 0.72, which meant 72% of the variations in GUDs can be explained by LST, indicating a significant correlation between GUDs and LST.
Besides, a rough idea of the range of GUDs and LST can be seen in LST, either for forestlands, grasslands, or croplands, indicating that the rising in temperature tends to advance the vegetation GUDs, in agreement with results reported by previous studies [1,16,27,33,41].The coefficient of determination (R 2 ) of linear fitted models between GUDs and LST for forests, grasslands, and croplands were similar, around 0.74, 0.72, and 0.71 with an average value of 0.72, which meant 72% of the variations in GUDs can be explained by LST, indicating a significant correlation between GUDs and LST.Besides, a rough idea of the range of GUDs and LST can be seen in Figure 4.As seen, LST over 2001-2019 varied between 272.4 k and 279.9 k, 272.7 k and 279.9 k, and 273.3 k and 279.9 k for forests, grasslands, and croplands, respectively, with the range of LST almost showing no difference.The scope of vegetation GUDs ranged from DOY 128 to DOY 101, DOY 136 to DOY 103, and DOY 118 to DOY 97 for forests, grasslands, and croplands, respectively, with the largest variations happening to grasslands, the earliest GUDs of DOY 97 given by croplands, and the latest DOY136 going to grasslands.

Spatial Patterns of the Correlation between GUDs and LST along the Rings
The vegetation GUDs and LST of forests, grasslands, and croplands were computed with each ring as a basic unit for every year from 2001-2019.The spatial pattern of the correlation coefficients between GUDs and LST from 2001 to 2019 were generated with each ring as a basic unit (Figure 6).The spatial pattern of the negative correlation between GUDs and LST along the urban-rural gradient followed a distance-decay pattern (Figure 6) and the extended distances of this declining trend were varied among vegetation types, about 32 km, 22 km, and 34 km beyond the perimeter of urban forests (Figure 6a,b), grasslands (Figure 6c,d), and croplands (Figure 6e,f), respectively.The correlation coefficient between GUDs and LST was slightly less negative in the ring near the urban area compared with those in the rural areas.A similar pattern was also reported by Meng et al.

Spatial Patterns of the Correlation between GUDs and LST along the Rings
The vegetation GUDs and LST of forests, grasslands, and croplands were computed with each ring as a basic unit for every year from 2001-2019.The spatial pattern of the correlation coefficients between GUDs and LST from 2001 to 2019 were generated with each ring as a basic unit (Figure 6).The spatial pattern of the negative correlation between GUDs and LST along the urban-rural gradient followed a distance-decay pattern (Figure 6) and the extended distances of this declining trend were varied among vegetation types, about 32 km, 22 km, and 34 km beyond the perimeter of urban forests (Figure 6a,b), grasslands (Figure 6c,d), and croplands (Figure 6e,f), respectively.The correlation coefficient between GUDs and LST was slightly less negative in the ring near the urban area compared with those in the rural areas.A similar pattern was also reported by Meng et al. [50], who found that the correlation coefficients between SOS (start of season) and air temperature were significantly weaker than rural ones in 43 cities in the United States.
Compared with those of grasslands (Figure 6c,d) and croplands (Figure 6e,f), the changing trend of correlation between GUDs and LST of forests (Figure 6a,b) was more regular, characterized by an increasingly significant correlation as distances were increasing.In the rings from 4 km to 50 km, the correlation coefficient (r) kept slightly fluctuating, with r ranging between 0.83 and 0.93, and the strongest correlation (r = −0.93)happened about 30-34 km away from the urban area.The GUDs of croplands also had a more significant correlation with LST in this area (32-34 km) than in other rings (Figure 6e,f), but those of the grasslands did the opposite (Figure 6c,d).From land cover (Figure 1) and the distribution of average GUDs and LST (Figure 4), it can be learned that there was a large area of croplands in the northwest of the rings 30-38 km surrounded with forests, and a lower LST and the later average GUDs of croplands in this area than those in the southeast of these rings were identified.Grasslands in these rings were mainly distributed near small built-up areas, where the temperature was slightly higher than other parts in these rings.So, it was no surprise that the correlation between GUDs and LST of grasslands in 32-40 km was similar to those in ring 2 km nearby the urban perimeter.
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 18 [50], who found that the correlation coefficients between SOS (start of season) and air temperature were significantly weaker than rural ones in 43 cities in the United States.Compared with those of grasslands (Figure 6c,d) and croplands (Figure 6e,f), the changing trend of correlation between GUDs and LST of forests (Figure 6a,b) was more regular, characterized by an increasingly significant correlation as distances were increasing.In the rings from 4 km to 50 km, the correlation coefficient (r) kept slightly fluctuating, with r ranging between 0.83 and 0.93, and the strongest correlation (r = −0.93)happened about 30-34 km away from the urban area.The GUDs of croplands also had a more significant correlation with LST in this area (32-34 km) than in other rings (Figure 6e,f), but those of the grasslands did the opposite (Figure 6c,d).From land cover (Figure 1) and the distribution of average GUDs and LST (Figure 4), it can be learned that there was a large area of croplands in the northwest of the rings 30-38 km surrounded with forests, and a lower LST and the later average GUDs of croplands in this area than those in the southeast of these rings were identified.Grasslands in these rings were mainly distributed near  As seen, the 2019-2001 differences in GUDs (∆GUDs) varied from 4.9 days to 20.1 days (forests), 9.8 days to 24.5 days (grasslands) and 2.7 days to 15.6 days (croplands), and the corresponding ∆LST were from 2.08 K to 2.94 K (forests), 2.27 K to 3.27 K (grasslands), and 2.20 K to 3.19 K (croplands), respectively.∆GUDs and ∆LST both showed the greatest values in 2 km ring nearest the urban perimeter for all three investigated vegetations, and then gradually decreased in rings further away (Figure 7), but the continuity of this declining trend varied between the vegetation types.
The ∆GUDs and ∆LST in forests almost keep the same pattern along the urban-rural gradient (Figure 7a), which were also demonstrated by the high coefficient of determination (R 2 = 0.7376) (Figure 7b).The similar changing trends in ∆GUDs and ∆LST were also observed in grasslands and croplands, but, differing from forests, these trends remained unchanged only within the range of 22 km away from the urban area for grasslands (Figure 7c) and croplands (Figure 7g) respectively.Within these ranges, ∆GUDs and ∆LST showed significantly positive correlation, with the biggest R 2 of 0.9059 (Figure 7d) given by grasslands, followed by R 2 of 0.7285 in croplands (Figure 7h).The positive correlation between ∆GUDs and ∆LST can be interpreted as the greater rising in LST being likely to lead to a larger advance of vegetation GUDs, i.e., a greater GUD change is associated with a higher temperature difference.Similar results were also reported by Zhou et al. [18] in China's major cities and by Zhang et al. [31] in eastern North American.
It is worth noting that the relationships between ∆GUDs and ∆LST in grasslands and croplands in rings 22-50 km, (Figure 7e,f,i,j), were either insignificant or opposite to the relationships that happened near the urban area.Specifically, ∆GUDs of grasslands 22 km away from the urban area were almost not related to ∆LST, and ∆GUDs of croplands 22 km away from the urban area showed a negative correlation with ∆LST.Further observations showed that the small ∆LST corresponded to big ∆GUDs in croplands 24 km to 36 km (Figure 7i), and that the large ∆LST was accompanied by small ∆GUDs in croplands 38 km to 50 km away from the urban area.The complicated spatial patterns in ∆GUDs and ∆LST of croplands may be due to being highly influenced by human activities, crop types, and human management practices [31,51].

Spatial Patterns of the Response Rate of ∆GUDs to ∆LST
To depict the response rate of ∆GUDs to ∆LST, the ratios of absolute ∆GUDs to ∆LST were computed, which were expressed as the number of days advanced in GUDs per 1 K increase in LST (unit: days K −1 ).The ratios were obtained for forests, grasslands, and croplands along the urban-rural gradient with each ring used as a unit.
Except for croplands, the most sensitive response of ∆GUDs to ∆LST in forests and grasslands both happened at the ring nearest the urban perimeter, in agreement with the greatest ∆GUDs and ∆LST in this ring, with the maximum rate of ∆GUDs to ∆LST being about 6.8 days K −1 and 7.5 days K −1 in forests and grasslands (Figure 8), respectively.The result is similar to the study conducted by Zhou et al. [18], who observed a strong phenological sensitivity to temperature with 9-11 days GUD advance per 1 • C increase of LST in spring in China's major cities.Some researchers have also reported that vegetation green-up dates advanced by 1-14 days per 1 • C increase of LST in urban centers compared to the corresponding rural surrounding regions in North America and China [3,18,52].
The ratios of ∆GUDs to ∆LST in forests showed a gradual decreasing trend from about 6.8 days K −1 in the ring nearest the urban area to the minimum 3 days K −1 in the ring 30 km away from the built-up region, which was then kept at a stable value of around 2.4-2.8days K −1 in the remaining areas.However, the ratios of ∆GUDs to ∆LST in grasslands and croplands decreased first within about 22 km beyond the urban perimeter and then climbed up to about 7.2 days K −1 (grasslands) and 6.7 days K −1 (croplands) at rings 46−50 km, respectively.The same turning points, distinguishing the upward and downward trend around 22 km further away from the urban perimeter, were identified in a spatial pattern of the response rate (Figure 8d,f) and in the changing trends of ∆GUDs and ∆LST in grasslands and croplands (Figure 7c,g), implying the necessity of paying more attention to the regions 22 km away from urban areas when exploring the effective distances of influences of urbanization on LST and vegetation GUDs.The ratios of ΔGUDs to ΔLST in forests showed a gradual decreasing trend about 6.8 days K −1 in the ring nearest the urban area to the minimum 3 days K −1 in th 30 km away from the built-up region, which was then kept at a stable value of aroun 2.8 days K −1 in the remaining areas.However, the ratios of ΔGUDs to ΔLST in grass and croplands decreased first within about 22 km beyond the urban perimeter and climbed up to about 7.2 days K −1 (grasslands) and 6.7 days K −1 (croplands) at rings km, respectively.The same turning points, distinguishing the upward and down trend around 22 km further away from the urban perimeter, were identified in a s

Correlationships of GUDs to Nighttime LST and GUDs to Daytime LST
As it is still up for debate about whether nighttime temperature or daytime temperature has a greater impact on GUDs, the correlations GUDs to nighttime LST and GUDs to daytime LST were compared (Figure 9) before the nighttime LST was finally adopted in the analysis.The comparison showed that both daytime LST and nighttime LST have a negative correlation with GUDs for all three vegetation types, with nighttime LST being found to be better correlated with vegetation GUD than daytime LST, which was especially true for forests and grasslands.For croplands, the relationships between nighttime LST, daytime LST, and GUDs were similar, with R 2 of the fitted model being lower than those of forests and grasslands.This may because croplands were highly influenced by human activities such as crop species selection.When vegetation types were not considered, vegetation GUDs had a more significant response to nighttime LST.In our opinion, there are two reasons for the better relationship of vegetation GUDs to nighttime LST: (1) vegetation GUD is mostly dominated by air temperature.Daytime LST is significantly higher than the corresponding air temperatures, while similar temperatures were found between nighttime LST and corresponding air temperatures [53]; (2) the nighttime LST is closer to the minimum temperature than daytime LST.There is evidence from other studies that (a) vegetation GUD was a low-temperature constraint and was mainly affected by the minimum air temperature [54], (b) nighttime temperature is a stronger modulator of the SOS (start of growing season) of deciduous broadleaved forest in northeast China than daytime temperature [33].The nighttime temperature tends to be more below the threshold temperature, and thus prevails over daytime temperature in affecting the fulfillment of the thermal requirement of vegetation in spring [55].Additionally, temperature records and future predictions all demonstrate that the temperature at nighttime presents a stronger warming than at daytime [33,56], which might contribute to the greater impact of nighttime temperature on vegetation GUDs than daytime temperature.

Different Spatial Pattern of Correlation between GUDs and LST along the Rings
On average, the results in this study showed that the correlation between GUDs and LST was the strongest in forest and was followed by croplands and finally grasslands.Wang et al. [33] also reported similar results in the study of relationships between vegetation phenology and climate factors in East China during the years 2001-2013.The phenomena that the GUDs of croplands were better correlated with LST in colder rings compared to warmer rings, and the correlation coefficients in the ring (0-2 km) nearby the urbandominate area were smaller than the surrounding areas, corresponds to the views that the vegetation green-up date is a low-temperature constraint [54] and that the higher sensitivity of vegetation phenology to temperature change is observed in colder regions than in warmer zones [18,57,58].
Additionally, some researchers have pointed out that apart from urban, topography also has large effects on solar radiation and land surface temperature [59,60].On the one hand, topography is likely to affect the retrieval accuracy of land surface temperature, and, on the other hand, it would complicate the interaction of urban effects and phenology.It is suggested that topography should be incorporated in future analysis to better describe the correlation between GUDs and LST along the gradient.

Different Response Rates of ∆GUDs to ∆LST along the Rings
The response rate of ∆GUDs to ∆LST varied along the gradient rings and among vegetation types, featured by a gradual decreasing trend in forests, and a downward then an upward trend in grasslands and croplands (Figure 8b,d,f).The correlation between ∆GUDs and ∆LST in grasslands were both negative, with the correlation of ∆GUDs and ∆LST being significant for rings 0-22 km (Figure 8d) but insignificant for rings 24-50 km (Figure 8f).It can be seen from Figure 7c that both ∆GUDs and ∆LST in grasslands showed a regular distance-decay pattern, with ∆GUDs gradually decreasing from 24.5 days to 11.8 days, and from 3.27 K to 2.50 K for ∆LST as the distance away from the urban area increasing from 0 to 22 km, indicating that the urban heat island-induced changes in GUDs and LST were weakening when further away from the urban domain.Meanwhile, in rings 24-50 km, ∆GUDs fluctuated around 14 days to 18 days with ∆LST being kept at around 2.4 K, both almost making no response to the increase in distances.This might mean that the extent of effect of urbanization on temperature warming and GUDs advancing in grasslands were most likely to be around 22 km in our case, and LST was most likely not the main influencing factor for changes in GUDs when the increase in LST was lower than 2.5 K.More work is needed to validate these views.The response rate in croplands 22-50 km beyond the urban perimeter showed an increasing trend with distances (Figure 8).The distribution of average LST (Figure 4) showed that the temperature in these regions was lower than near the urban areas.This meant that areas above 20 km beyond the urban perimeter were generally colder than areas near the urban areas.The upward trends of the response rate in croplands in a colder region demonstrated that vegetation ∆GUDs were more sensitive to ∆LST in colder regions when compared to warmer areas.This phenomenon was in agreement with previous studies [18,57,58].Prevey et al. [61] and Meng et al. [50] attributed the higher temperature sensitivity of GUDs in colder areas than in warmer areas in part to the small absolute changes in temperature constituting greater relative changes in the thermal balance under colder conditions.
negative correlation with GUDs for all three vegetation types, with nighttime LST being found to be better correlated with vegetation GUD than daytime LST, which was especially true for forests and grasslands.For croplands, the relationships between nighttime LST, daytime LST, and GUDs were similar, with R 2 of the fitted model being lower than those of forests and grasslands.This may because croplands were highly influenced by human activities such as crop species selection.When vegetation types were not considered, vegetation GUDs had a more significant response to nighttime LST.In our opinion, there are two reasons for the better relationship of vegetation GUDs to nighttime LST: (1) vegetation GUD is mostly dominated by air temperature.Daytime LST is significantly higher than the corresponding air temperatures, while similar temperatures were found between nighttime LST and corresponding air temperatures [53]; (2) the nighttime LST is closer to the minimum temperature than daytime LST.There is evidence from other studies that (a) vegetation GUD was a low-temperature constraint and was mainly affected by the minimum air temperature [54], (b) nighttime temperature is a stronger modulator of the SOS (start of growing season) of deciduous broadleaved forest in northeast China than daytime temperature [33].The nighttime temperature tends to be more below the threshold temperature, and thus prevails over daytime temperature in affecting the fulfillment of the thermal requirement of vegetation in spring [55].Additionally, temperature records and future predictions all demonstrate that the temperature at nighttime presents a stronger warming than at daytime [33,56], which might contribute to the greater impact of nighttime temperature on vegetation GUDs than daytime temperature.

Different Spatial Pattern of Correlation between GUDs and LST along the Rings
On average, the results in this study showed that the correlation between GUDs and LST was the strongest in forest and was followed by croplands and finally grasslands.

Conclusions
In this study, the influences of LST on vegetation GUDs were analyzed using long-term (2001-2019) MODIS products from two perspectives, i.e., the spatial patterns of correlation between GUDs and LST and the response rate of ∆GUDs to ∆LST.The highlight of the study is the consideration of vegetation types and distances away from the urban domain in the analysis of the relationships.The main conclusions are as follows: (1) Vegetation GUDs show divergent correlations with LST according to vegetation types and distances away from the built-up areas, featuring a negative and relative weaker correlation near the urban area (0-10 km) than the surrounding area.GUDs are more sensitively correlated to LST in colder rings than in warmer rings; (2) The magnitude of ∆GUDs to ∆LST are positively correlated, with the advanced days in GUDs and the rising in LST both being larger near the urban domains when compared to these in rings 22 km further away from the urban perimeter, proving a solid influence of urban warming on the changes in vegetation GUDs; (3) The spatial pattern of the response rate of ∆GUDs to ∆LST demonstrated a distinct trend in different vegetations and distances, characterized by a distance-decay trend in forests and a downward trend first and then an upward trend in both grasslands and croplands.Higher sensitivity of ∆GUDs to ∆LST is identified in colder sites than in warmer rings, especially for croplands.
Although a detailed analysis about the relationships between vegetation GUDs and LST of forests, grasslands, and croplands were carried out at different distances away from the urban domain area, and some quantitative results were obtained, more efforts are needed to validate and consolidate the conclusions above by incorporating more vegetation types, climate factors, and other remote sensing data, such as the new MODIS LST product.

Figure 1 .
Figure 1.The delineation of study area.

Figure 1 .
Figure 1.The delineation of study area.

Figure 2 .
Figure 2. The delineation of 25 ring zones along the built-up area of Beijing.

Figure 2 .Figure 3 .Figure 3 .
Figure 2. The delineation of 25 ring zones along the built-up area of Beijing.2, 14, x FOR PEER REVIEW 6 of 18

Figure 4 .
Figure 4. Spatial patterns of 19-year average (a) green-up dates and (c) land surface temperature over the period 2001-2019 and the area percentage charts for different GUDs (b) and LST (d).

Figure 4 .
Figure 4. Spatial patterns of 19-year average (a) green-up dates and (c) land surface temperature over the period 2001-2019 and the area percentage charts for different GUDs (b) and LST (d).

Figure 4 .
As seen, LST over 2001-2019 varied between 272.4 k and 279.9 k, 272.7 k and 279.9 k, and 273.3 k and 279.9 k for forests, grasslands, and croplands, respectively, with the range of LST almost showing no difference.The scope of vegetation GUDs ranged from DOY 128 to DOY 101, DOY 136 to DOY 103, and DOY 118 to DOY 97 for forests, grasslands, and croplands, respectively, with the largest variations happening to grasslands, the earliest GUDs of DOY 97 given by croplands, and the latest DOY136 going to grasslands.

Figure 6 .
Figure 6.Spatial patterns of correlation between vegetation GUDs and LST in forests (a,b), grasslands (c,d), and croplands (e,f) along the urban-rural gradient for the period 2001-2019.

Figure 6 .
Figure 6.Spatial patterns of correlation between vegetation GUDs and LST in forests (a,b), grasslands (c,d), and croplands (e,f) along the urban-rural gradient for the period 2001-2019.

Figure 7 .
Figure 7. Changing trends of the 2019-2001 differences in vegetation absolute GUDs (ΔGUDs) and in LST (ΔLST) along the urban-rural gradient and their relationships for forests (a,b), grasslands (cf), and croplands (g-j).

Figure 7 .
Figure 7. Changing trends of the 2019-2001 differences in vegetation absolute GUDs (∆GUDs) and in LST (∆LST) along the urban-rural gradient and their relationships for forests (a,b), grasslands (c-f), and croplands (g-j).

Figure 8 .
Figure 8. Spatial pattern of the response rate of absolute ΔGUDs to ΔLST in forests (a,b), gras (c,d), and croplands (e,f) along the urban-rural gradient.

Figure 8 .
Figure 8. Spatial pattern of the response rate of absolute ∆GUDs to ∆LST in forests (a,b), grasslands (c,d), and croplands (e,f) along the urban-rural gradient.

18 Figure 9 .
Figure 9.The relationship of vegetation green-up dates (GUDs) to daytime and nighttime land surface temperatures (LST) along the urban-rural gradient.
Wang et al. [33] also reported similar results in the study of relationships between vegetation phenology and climate factors in East China during the years 2001-2013.The phe-

Figure 9 .
Figure 9.The relationship of vegetation green-up dates (GUDs) to daytime and nighttime land surface temperatures (LST) along the urban-rural gradient.