Use of local climate zones to assess the spatiotemporal variations of urban vegetation phenology in Austin, Texas, USA

ABSTRACT Phenological changes caused by urbanization may provide evidence of how vegetation responds to global warming. However, the phenological characteristics of vegetation in metropolitan areas have been poorly studied, especially in terms of spatiotemporal variations. In this study, we explored the applicability of Local Climate Zones (LCZs) to investigate the spatiotemporal variations in urban vegetation phenology by linking MODIS-derived phenology metrics with LCZs in the Austin metropolitan area in Texas, USA. We extracted three vegetation phenology metrics from MODIS data between 2008 and 2018, including the start of growing season, end of growing season, and length of the growing season (i.e. SOS, EOS, and LOS, respectively). The results showed that during the study period, the EOS and SOS gradually advanced, while LOS showed no obvious change. Statistical analysis was conducted to examine the spatiotemporal variations of the phenology metrics among different LCZs and along “Urban-Rural Gradients” (URGs). There were 37.5%, 75.0%, and 74.3% pairs of LCZs, indicating statistically significant phenological differences in terms of SOS, EOS, and LOS in 2012, respectively. In contrast, most pairs of URGs showed almost no differences in phenological metrics, especially in EOS. Geographically, SOS showed a fluctuating change with an advancing tendency, whereas the EOS decreased very slowly with distance from the city center (i.e. along the URGs). LCZs can be used to help identify distinctive phenology metrics with statistically significant differences, especially in EOS and LOS. Compared to URGs, LCZs offer a unique analytical framework for studying urban ecosystem patterns, functions, and dynamics. Lastly, LCZs can enable the identification of sensitive areas for ecological protection in support of sustainable urban development and environmental stewardship.


Introduction
Vegetation phenology refers to plant cycle events triggered by seasonal changes (Richardson et al. 2018;Neil, Landrum, and Wu 2010). It is an important indicator of vegetation growth (White et al. 2009) and climate change (Zhou et al. 2016). Land cover changes resulting from urbanization, such as impervious surfaces and high-density buildings, are among the factors that cause local climate change (Roetzer et al. 2000). Urbanized regions are ideal places for studying the global warming mechanisms that influence vegetation phenology because of their relatively high temperatures (Luo et al. 2006;Mimet et al. 2009;Richardson et al. 2013). Rapid urbanization also intensifies the urban heat Island (UHI) effect, defined as an urban area experiencing higher air temperature than the surrounding rural area (Oke 2002). Studying the spatiotemporal variations of urban vegetation phenology can help to develop mitigation strategies to reduce the UHI thermal effect.
Urbanization occurs at the expense of natural land cover and strongly disturbs the plant growth environment, leading to large impacts on vegetation phenology in urban and surrounding areas. Urbanization reorganizes surface runoff networks and watershed hydrology by changing impervious surface coverage and subsurface drainage pathways (Price 2011). The impervious surface disturbs rainfall infiltration and may further reduce groundwater recharge (White and Greer 2006), leading to detrimental effects on plant growth. In addition, urban soil can be more compact and hydrophobic than rural soil and heavy metals accompanied by solid waste and water pollution can heighten the environmental stresses experienced by urban vegetation (Pouyat, McDonnell, and Pickett 1995). Urban soil is more heterogeneous than rural soil (Tang et al. 2011). While the relatively high concentration of CO 2 and the UHI effect tend to enhance vegetation growth, the concentration of gaseous and photochemical pollutants and the lower surface solar radiation obtained could inhibit vegetation production (Alpert et al. 2005;Gregg, Jones, and Dawson 2003). Monitoring urban vegetation phenology is important in understanding how climate change mechanisms influence vegetation cover and how plants respond to urban environmental changes.
Ecological field work and satellite remote sensing are the two main approaches for monitoring vegetation conditions and temporal variations (Li et al. 2017;White et al. 2002). Satellite imagery has increasingly been used to extract vegetation properties and monitor vegetation phenology because of its advantages in obtaining largescale and spatially explicit information. Results from satellite imagery indicate that the phenological changes caused by urbanization can promote the growth of vegetation, causing an earlier start of growing season (SOS) and significantly longer length of the growing season (LOS) in urbanized areas. For example, in USA (Neil, Landrum, and Wu 2010;White et al. 2002;Ziska et al. 2003), Europe (Rahman, Armson, and Ennos 2014;Roetzer et al. 2000), and China (Lu et al. 2006;Zhou et al. 2016), the SOS of plant species in urban areas takes place significantly earlier than in rural areas by several days to weeks. However, most previous urban phenology studies concentrate on the "urban-rural" binary classes, which is inadequate for systematically acquiring heterogeneous phenological characteristics across the entire metropolitan area.
Consequently, some recent studies have followed the "Urban-Rural Gradients" (URGs) zoning method to investigate the localized patterns in vegetation phenology in metropolitan areas. For instance, distinct urban-rural phenological variations over time and space have been detected in coastal Dalian, China (Yang et al. 2020). The effects of urbanization on phenology can vary based on geographic location (Diamond et al. 2014), urbanization intensity (Li et al. 2017), vegetation species, and changes in the urban environment (Jia et al. 2021). The URG method is straightforward and efficient for studying metropolitan areas at large scales. However, there is a demand for exploring and integrating different zoning approaches to comprehensively capture the spatiotemporal variation of urban vegetation conditions.
In contrast to the limited studies on spatiotemporal variations of urban phenology, the micro-scale effects and the spatial arrangement of urban canopy components on the surface UHI variation have been widely studied (Connors, Galletti, and Chow 2013;Weng 2009). Stewart et al. (2012) proposed a Local Climate Zones (LCZ) classification scheme to describe urban and rural landscapes from the perspective of thermal climate characteristics. LCZs are defined as "regions of uniform land cover, surface structure, construction material, and human activity that span hundreds of meters to several kilometers on a horizontal scale" (Stewart et al. 2012). LCZs are a robust, objective, and universal approach for characterizing urban forms and natural land cover for use in climatological research based on the zoning practices of urban spaces. Therefore, the LCZ method makes it possible to conduct a standardized and comprehensive spatiotemporal analysis of vegetation phenology in urban areas.
Recent studies have investigated the impact of UHI intensities (i.e. the quantitative thermal difference along the complete urban-rural gradient spectrum) on vegetation phenology when using LCZ mapping. These studies have provided a theoretical basis for using LCZs to study spatial differences in vegetation phenology (e.g. Kabano, Lindley, and Harris 2021). However, the impacts of urbanization and UHI on vegetation phenology at the neighborhood/local scale remain under-studied. Additionally, the relationship between vegetation phenology and the threedimensional structure of metropolitan areas has rarely been studied due to the complexities in urban morphology and lack of data for effective characterization. In this study, we tested the applicability of the LCZ classification system to investigate urban vegetation phenology. First, a framework was built to link remote sensing-derived vegetation phenology indicators with LCZ categories. The term LCZ mapping was introduced as the technique to implement the LCZ classification scheme in this study. The spatial variation of vegetation phenology in the Austin metropolitan area, Texas, USA, was then examined and compared using LCZ and URG methods. This study aims to develop a method for using LCZs to understand the spatiotemporal variations in vegetation phenology in rapidly urbanizing metropolitan areas.

Study area
Our study area is the Austin metropolitan area in central Texas,. The study area is characterized by its unique and narrow transitional climatic zone, intersecting with dry deserts in the southwest and humid regions in the southeast of the USA. Highway I-35 physically separates the study area from the Blackland Prairies in the east to the Edwards Plateau in the west. Correspondingly, the topography and vegetation distribution of the study area demonstrated contrasting characteristics on the east and west sides of the city of Austin (Figure 1). With a humid subtropical climate, Austin has hot summers, short and mild winters, and a pleasantly warm spring. During the summer (June to mid-September), the temperature is 32°C or warmer. The annual rainfall is 889 mm, which is close to the average value for USA. The number of sunny days per year in Austin reach 228 days. The highly variable humidity and prolonged summer heat lead to diverse and dense vegetation cover. Live oak, shinnery oak, mesquite, and juniper dominate the woody vegetation throughout the Edwards Plateau, with relatively high slope variation. Most fertile areas in the Blackland Prairies have been cultivated, and only small acreages of grassland remain in their original forms.
The urbanized area accounts for approximately 771 km 2 and the urban population was approximately 1 million people in 2019. Since the mid-2000s, Austin has adopted a series of ordinances and regulations for compacted urban design (Austin City Council 2012), which is still less compact than its comparable cities, such as Houston and San Antonio in Texas, in terms of population density and building volume (Zhao, Weng, and Hersperger 2020). These characteristics of climate, vegetation cover, and urban forms make Austin a representative area for studying the spatiotemporal distribution of urban phenology.

Data and methods
The spatiotemporal investigation of vegetation phenology in this study was based on linking remote sensing-derived phenology indicators with zoning methods based on LCZs and URGs. Figure 2 shows the data processing steps for LCZ mapping, vegetation phenology calculation, and spatial statistical analysis based on the LCZs and URGs.

LCZ mapping
The three-dimensional urban and natural landscape structure in the Austin metropolitan area was investigated using the LCZ types devised by Zhao et al. (2019). Lidar data were used to extract threedimensional information for the surface roughness features, especially for buildings and vegetation. This information was combined with information from the National Land Cover Database (NLCD) to prepare three major surface cover properties (height of roughness, building surface fraction, and impervious surface fraction) as the main geometric and surface cover properties (Table A1 in the Appendix). Other geometric and surface cover properties were also collected and calculated based on remote sensing imagery, NLCD, and the city's urban planning datasets. The LCZ types were delineated using the top-down decision-making algorithm (Zhao et al. 2019). Considering the dominant role of building structures in the local climate, the LCZ mapping scale was determined to be 270 m. This was determined by investigating the ranges of the variograms for building heights. The sky view factor, aspect ratio, and terrain roughness were calculated to quantify the accuracy of the LCZ maps. The boxplot of the sky view factor and aspect ratio was consistent with results from a prior GISbased LCZ mapping study (Zheng et al. 2018). In addition, LCZ mapping results were linked with ERDAS and overlaid with the historical imagery provided by Google Earth, and LCZs matched with the examples provided in Table 2 on p. 1885 of Stewart et al. (2012). Figure 3 shows the LCZ distribution and four dominant LCZ examples illustrated on Google Earth.
Austin's LCZ "built-up" types (i.e. LCZ 1-10) are mostly concentrated in the center of the study area and are surrounded by "natural" types (i.e. LCZ A-G). The distribution of LCZ "natural" types indicates that the west is mainly LCZ A (dense trees), the east is mainly LCZ D (low plants), and the south is mainly LCZ C (bush, shrub). Table A2 and Figure A1 show the statistical summary of the individual LCZ patches and the proportion of each LCZ in the study area. LCZ D (low plants) cluster pattern was the most predominant feature in the study region, accounting for approximately 29.56% of the Austin metropolitan area and was the highest proportion of all LCZs.

Derivation of vegetation phenometrics
Remotely sensed vegetation phenology was characterized by a time series of spectral vegetation indices. Both the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) can capture seasonal changes in vegetation properties and productivity, while EVI is more suitable for monitoring the phenology of high biomass crops. These vegetation indices make it possible to monitor vegetation phenology events from local to large scales. The Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 (C6) Land Cover Dynamics Product (MCD12Q2) provides information related to global phenology at a 500 m resolution and an annual time step. This was applied as the data source for our study. The product was derived from MODIS observed 2-band EVI time series data, which was calculated from MODIS nadir bidirectional reflectance distribution function (BRDF) adjusted surface reflectance (NBAR-EVI2) (Gray, Sulla-Menashe, and Friedl 2019). In contrast to the C5 product, algorithms were designed to better capture phenometrics with multiple vegetation cycles in the C6 product; up to two detected vegetation cycles can be recorded to detect multi-cropping.
First, the 2-band EVI time series data were preprocessed to reduce the effects of clouds, aerosols, snow, and observation angles. A search sliding window (EVI values for five consecutive phases) was designed to determine the continuous increase and decrease intervals to fill the dormant period values. When the maximum value and amplitude of the interval met the given threshold conditions, the increase and decrease interval was determined to be a growth cycle by fitting a quality assurance (QA)/quality control (QC)weighted penalized cubic smoothing spline (Gray, Sulla-Menashe, and Friedl 2019). The extreme points of curvature changes were identified to calculate greenup, midgreenup, peak maturity, midgreendown, senescence, and dormancy (Table A3, Appendix). The QA scores were also calculated for each phenometric and for the entire greenup/ down segment in the MCD12Q2 product. The MODIS Reprojection Tool (MRT) software was used to splice and project the MODIS imagery. The phenology metrics were originally provided as days since January 1 st , 1970, and were converted to ordinal dates (i.e. day of year). The greenup and dormancy properties from MCD12Q2 were extracted to represent the SOS and the end of the season (EOS) ( Table A3, Appendix). The QA scores of greenup and dormancy properties were checked for the study area, and the proportion of pixels with best QA scores ranged from 99.98% to 100% during the study period. The LOS index was calculated as the difference between the EOS and SOS dates.

Statistical analysis of phenology indicators based on LCZs and URGs
The division of URGs was based on the geometric center point of LCZ 1 (compact high-rise), which was assumed to be the urban center of Austin and the starting point of the URGs. To ensure that the URGs were significant, 2000 m was used as the radius gradient to divide the Austin metropolis into 11 gradients. The significance was p < 0.05, where the effective pixel proportion in terms of the three phenology indicators for different years was greater than 95%. The gradients include gradient 2 K, gradient 4 K, . . ., gradient 20 K, and gradient 20 K + (i.e. the area outside the circular area with a radius of 20,000 m) (Figure 3). The calculated phenology indicators were superimposed onto the LCZs and URGs. The "zonal statistics" tool in ArcGIS software was used to acquire the pixel values of vegetation phenology indicators for the two respective zoning methods, which linked the indicators with the LCZ and URGs zoning methods.
The R statistical package was used to analyze the LOS in different LCZs. The average values in 17 different LCZs and distance gradients were visualized. We performed a one-way analysis of variance (ANOVA) to determine the LOS variation among different LCZs and URGs. Next, Tukey's HSD (Honest Significant Difference) test (Tukey 1953) was used to evaluate the significance of the LOS difference between certain pairs of LCZs and URGs. The formula to calculate the HSD is as follows: where q is related to the significance level α, the number of averages in the experiment (k), and the degrees of freedom of error n À k. The difference between the pair of averages is significant if the difference between any pair of averages exceeds the HSD.

Spatiotemporal characteristics of vegetation phenology
The results showed that the SOS occurred in March in the west of Austin and in February in the east, from the 30 th to the 60 th day. The SOS gradually advanced from west to east throughout the study period (Figure 4). The spatial distributions of the EOS in 2012 and 2018 were relatively similar. The EOS in the western area was concentrated in November, while the EOS in the east was concentrated in July and August. The EOS generally advanced from 2010 to 2018. Vegetation in the middle of the study area had the longest LOS (more than 300 days) and this was shorter for vegetation in the west and east (approximately 180-270 days). Overall, the EOS and SOS gradually advanced while LOS showed no obvious trend during the study period.

Vegetation phenology by LCZ types and along URGs
Considering the desirable quality of data for our study area according to phenometric-specific QA/QC information provided by the MCD12Q2 product, we assume that the vegetation phenology pixels in each LCZ can sufficiently reflect the vegetation phenology characteristics for statistical tests. Figures 5 and 6 display box plots of LOS values for different LCZs and URGs, respectively, from 2008 to 2018. The LOS was observed to vary according to LCZ and year, as shown in Figure 5. The longest LOS was generally observed in LCZ 1-4 "high-rise" types, especially in 2009, 2011, 2012, and 2018. Most of the LOS in 2009 was significantly lower than in other years, which may be due to the severe drought in Texas in the aforementioned year. For LCZ 1 (compact highrise), LCZ 4 (open high-rise), LCZ 9 (sparsely built), and LCZ A (dense trees), the box plots were leftskewed (i.e. mean values close to the lower quartile). Most of the LCZ 1 (compact high-rise) type and a considerable amount of LCZ 4 and LCZ 9 were distributed inside the city, and the adverse effects of drought could be alleviated with human intervention, which led to some higher values of LOS. LCZ A (dense trees) was considerably distributed in the Edwards Plateau in the western part of the study area. Besides the exposed limestone, the soil is dominated by calcareous rubble and clay and the water storage capacity is relatively strong (Texas 2021). Additionally, the dominant woody vegetation reduces water loss due to transpiration by defoliation. The plants utilize residual moisture to rehydrate the stem tissue (Borchert 1994); thus, the LOS can be extended for some vegetation types classified as LCZ A (dense trees). In 2012, the LOS variation among LCZs was the largest and most significant, which may be related to the 2012 Austin LCZ map used in this study. Furthermore, LCZ E (bare rock or paved) and LCZ F (bare soil or sand) are more susceptible to high temperatures because of the low thermal heat capacity of water relative to other properties; thus, the LOS values tend to be lower with high variations compared to those in other years.
The LOS differences among URGs were less apparent than those based on LCZ types, and the LOS range within each URG was smaller and included fewer outliers ( Figure 6). The shortest LOS was mostly found for Gradient 20 K+; however, there was less variation among different URGs. LOS also varied annually. Nevertheless, the phenological characteristics of LCZs and URGs showed apparent annual variations ( Figures  5 and 6). Besides the LOS variations for some URGs for the years 2009 and 2012, there were considerable variations for Gradient 20 K + . The zoning of the Gradient 20 K + includes all the outside circular area with a radius of 20,000 m, where the diversity of plant cover is relatively high, and the LOS can be different for different vegetation types. The LOS differences in 2015 were the smallest, perhaps because of the low availability of vegetation phenology pixels extracted in that year.
The p values produced by ANOVA from 2008 to 2018 indicate variations in the LOS among the LCZs and URGs. The HSD test was used to compare the 11-year variance analysis results and visualize the LOS differences among the LCZs. Compared with the Kruskal-Wallis and post hoc Conover's test, the parametric one-way ANOVA and Tukey HSD analysis indicated the same findings, and the latter was used to show variations. Most pairs of LCZs indicate differences in terms of LOS (i.e. boxes in yellow, green, and light blue). The differences in 2015 were smaller, perhaps because of low availability of vegetation phenology pixels extracted in 2015 due to cloud coverage. Figure 7 visualizes the significance of phenological differences between all pairs of LCZs and URGs for 2012. Basically, the differences between pairs of LCZs and URGs were more noticeable for EOS than for SOS. Among the 136 LCZ pairs, there were 51, 102, and 101 pairs of LCZs, indicating statistically significant phenological differences in terms of SOS, EOS, and LOS, respectively (i.e. 37.5%, 75.0%, and 74.3%, respectively). There were nearly 25 pairs of built-up LCZs, indicating statistically significant LOS differences (Figure 7). In contrast, for the area within 10, 000 meters (i.e. URG 2 K to URG 12 K) and area in the ring between 14, 000 meters and 20, 000 meters (i.e. URG 14 K to URGs 20 K) from the city center, URG pairs did not indicate statistical phenological differences.

The distribution of vegetation phenology and LCZs along URGs
The spatial characteristics of the 2012 phenology data were further analyzed to compare phenology distributions among the LCZ types and distance gradients. The average SOS value gradually increased from the center to the periphery, indicating that the LCZ "builtup" types had an earlier SOS, while the LCZ "natural" types had a later SOS (Figure 8a). LCZ D had the latest SOS, followed by dense trees, bushes, and scrubs. The central region had the largest EOS values, which gradually became smaller when moving away from the center, indicating that the EOS in LCZ "built-up" types occurred later, while EOS was earlier in the "natural" type ( Figure 8b). Additionally, the EOS occurred most frequently in dense trees, followed by bush, scrubs, and low plants. The LOS map for the year 2012 in Figure 4 shows that the LOS values for the central region were significantly larger than those of the surrounding regions, and gradually became smaller when moving away from the center. In other words, the LOS was longer in the LCZ "built-up" types, and the LCZ "natural" types had a shorter LOS (Figure 8c). The LOS for LCZ A was significantly larger than that for LCZ C, whereas LCZ D had the smallest LOS. Figure 9 shows the statistical distribution of LOS in 2012 for different LCZs and URGs. Most of the LOS in the LCZ "built-up" types were greater than 270 days, while the LOS in the "natural" types could be less than 220 days. The average SOS of "built-up" LCZs fell on the 52 nd day, EOS was on the 296 th day, and the LOS was on the 265 th day. In terms of vegetation phenology in the LCZ "natural" types, the mean SOS and EOS were the 88 th and 224 th days ( Figure A2). The latest SOS occurred for LCZ D, followed by LCZ A, and LCZ C. The EOS occurred earliest in LCZ D, followed by LCZ C, and finally occurred in LCZ A. The LOS for LCZ A was significantly longer than that for LCZ C, and the LOS for LCZ D was the shortest. The LOS for the LCZ "built-up" types   was significantly longer than that of the "natural" type. Hence, urbanization promoted the extension of the growing season, which generally increased toward the city center.
As shown in Figure 10, there was a fluctuating change in SOS with an increasing tendency based on the distance from the city (i.e. URGs). EOS decreased very slowly along the URGs, but it clearly advanced when the distance from the city center was greater than 20,000 m. Generally, it was observed a farther LOS from the city center indicated a shorter LOS. However, the URG zoning methods can reflect changing trends in vegetation phenology in terms of SOS but not consider other phenological characteristics. Combined with the proportion of different LCZs along URGs, it was found that LCZs "built-up" types (i.e. LCZ 1,4,7,9) were distributed at various distance gradients and the proportional change tendency was not in line with the distance from the centers (i.e. URGs). For instance, LCZ 1 (compact high-rise) was mostly distributed in the downtown area, but was also distributed in strips in Austin (Figure 2). On the one hand, the LCZ "built-up" types mostly covered the areas within 6,000 m from the city center with the highest variations of the proportional destitutions among different LCZ "built-up" types, which led to considerable fluctuations in phenological indicators. On the other hand, the proportion of the LCZ "built-up" types showed a small change from the URG 16 K to URG 20 K, which leads to tiny changes in phenological indicators. These findings were consistent with the Tukey's HSD test, as shown in Figure 7.

Local climate zones vs. urban-rural gradients
Overall, the LCZ method is more suitable for comprehensively interpreting the spatial variation in vegetation phenology. The results of the analysis are consistent with the fact that Austin extends from the central city to the suburban area from inside to outside, and the LCZs "built-up" type are distributed in strips in the city center, rather than in a simple concentric circle (Figures 2 and 10). The results of SOS and LOS in the vegetation phenology reflect the change in vegetation phenology with consideration of the distance from the city center. To a large extent, this is  related to the fact that the distribution of the "builtup" types in the study area are not concentrically distributed. There are also LCZ "built-up" types far away from the city, which weakens the trend in EOS change. Altitude, slope, aspect, and other influencing factors such as topography and soil conditions affect surface temperatures, which have a certain effect on the spatial differences in vegetation phenology. However, the method of URGs was only delineated from the horizontal distance perspective, without considering the vertical factors, such as altitude, slope. Thus, the results can be considered one-sided. Based on the three-dimensional (3-D) structure of the city, we comprehensively considered the 3-D information and divided the Austin metropolis into multiple LCZs according to different factors. The vegetation phenology of the different LCZs was then studied. In summary, the zonal statistics combined with the 3-D structure of the city are better than the URGs.
SOS showed a fluctuating change with an advancing tendency, whereas the EOS decreased very slowly based on the distance from the city (i.e. URG). This phenomenon was apparent for areas within 20,000 m from the city center. SOS was significantly delayed. LCZ "built-up" types promote the advancement of SOS and the delay of EOS with considerable fluctuations, especially for SOS. This may be largely caused by the higher coverage of impervious surfaces in the "built-up" type, higher air temperature than that of the "natural" type, and the significant UHI phenomenon. It may also largely offset the decline in ecosystem productivity caused by the replacement of construction land and sparse urban vegetation to a large extent. These findings confirm the dependence of vegetation phenology on the expansion of metropolitan urbanization. These results are consistent with the patterns found in Eastern USA (Meng et al. 2020), Northeast China (Yao et al. 2017), many urban areas in China (Zhou et al. 2016), and other northern mid-latitude large cities (Qiu et al. 2020). Urban canopy air temperature is affected by streets and building forms, building materials, surface permeability, and urban anthropogenic heat. Additionally, urban areas with higher impervious surface coverage have higher sensible heat fluxes which contribute to the UHI effect, causing the temperatures in urbanized areas to be higher than in surrounding areas. Because spring phenology is mainly regulated by temperature, an increase in temperature causes earlier growth in spring (Körner and Basler 2010;Menzel et al. 2006). In addition, the UHI also affects the autumn phenology (EOS) (Richardson et al. 2013;Walther et al. 2002). Autumn phenology determines the LOS and extends the LOS. This is consistent with the results of this study.
LCZs describe the shape, structure, and coverage of urban land surfaces and provide a comparable urban land use classification scheme (Stewart et al. 2012). Urban areas with higher impervious surface coverage will increase sensible heat fluxes, resulting in the UHI affecting LCZs, which is in contrast to natural surfaces that have higher pervious coverage (Avissar 1996). The LCZs consider 3-D urban information and offer a more comprehensive insight into the complexity of vegetation phenology. Statistical analysis revealed significant differences in vegetation phenology among different LCZs. Hence, LCZs can not only be used to study urban ecosystems and their response to UHI (e.g. Parece and Campbell 2018), but also be used to understand the feedback mechanisms that urban canopies may have on the environment and climate.

Uncertainties of LCZs for vegetation phenology analysis
The datasets and methods may introduce uncertainties and affect the interpretation of results. The MCD12Q2 product used in this study has a longer period and a resolution of 500 m. Compared with LCZs with a resolution of 270 m, errors will occur when using LCZs to explore vegetation phenology. In this study, the footprints and heights of individual buildings were extracted from Lidar and served as an essential input for LCZ mapping. The average building surface fraction was calculated as the proportion of building footprints in a 30 × 30 m area and was further aggregated to obtain the LCZ mapping. According to Stewart et al. (2021), LCZ 8 should is prevalent in North American cities owing to its large footprint. However, from our mapping scheme, there are not many areas where the building surface fraction ranges from 30 to 50 and while the height of roughness is between 3 and 30 m. Areas with a building surface fraction in this threshold can be classified as LCZ 5 (open mid-rise) or LCZ 6 (open low-rise). The mismatch between vegetation phenology data (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and LCZ mapping (2012) may also introduce uncertainties, especially considering the rapid urbanization process in the USA. Therefore, a higher spatial and temporal resolution is required to reduce uncertainty.
In addition to topography and landforms, other factors are also related to urbanization. For instance, species composition (Liu et al. 2018), hydrological conditions (Guan et al. 2014), photoperiod (Körner and Basler 2010;Zhang, Tarpley, and Sullivan 2007), management measures (Wang 2020), and atmospheric environment, including CO 2 concentration (Busetto et al. 2010) may also affect phenology. Urban vegetation often includes non-native species, which may exhibit phenological patterns differently compared to native plants. Urban vegetation is often strictly managed, and may also affect phenological patterns. In particular, cultivated land introduces great uncertainty to our results because it depends, to a large extent, on humans and the local climate (Lee, Kim, and Yoon 2019). Therefore, the driving forces of multiple factors will be further analyzed in future studies.
Our study quantified the impact of urban structures on vegetation phenology. However, it is challenging to explore the interactions between potential environmental driving factors. Future studies should focus on unraveling the potential mechanisms of vegetation on spatial differences in urban structure. This study provides evidence for LCZs as the basis for exploring the spatial differences in vegetation phenology in rapidly urbanizing metropolitan areas and provides a possible reference for studying the phenological response to continuous urban warming. Nevertheless, considering that the focus of this study is on the comparison of LCZ mapping and URG zoning in studying the phenological variation, we did not extend the study to other classification schemes. Neither did we apply the study to other cities.

Conclusions
This study explored the applicability of LCZs in investigating the spatiotemporal variations of urban vegetation phenology by linking MODIS-derived phenology indicators with LCZs in the Austin metropolitan area, Texas, USA. It confirmed the response mechanism of urbanization to vegetation phenology. The results showed that, in 2012, the SOS in the LCZ "built-up" types was about 35.73 days earlier than in LCZ "natural" types, and EOS was delayed by about 72.37 days in the LCZ "built-up" types. A considerable number of LCZ pairs showed statistically significant phenological differences, especially in SOS and LOS, which indicated that LCZs can help identify urbanized areas with distinctive and homogeneous phenological indicators. LCZ mapping provides rich 3-D information which can be used to standardize the description of the urban morphology of a city and is valuable for studying the spatial characteristics of vegetation phenology. Additionally, compared with previous distance gradient methods, this study found that zonal statistics based on LCZs are significantly more accurate than distance gradients.
With the tools used in this study, the study provides a novel approach to reveal the spatiotemporal dependence of vegetation phenology in a metropolitan area. Combining LCZs and vegetation phenology indicators contributes to the identification of ecological protection barriers to enhance ecosystem services and support the planning and management of land resources for sustainable urban development.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported the National Natural Science Foundation of China [No. 42101465].

Data availability statement
The data and codes in this manuscript are freely available in GitHub at https://github.com/wangsiyu99/Spatiotemporal_ investigation_of_the_vegetation_phenology.