Evidence of urban heat island impacts on the vegetation growing season length in a tropical city Landscape and Urban Planning

Knowledge about the impacts of urban heat islands (UHI) and associated thermal gradients on vegetation sea- sonality (i.e. phenology) is vital for understanding spatial patterns in vegetation ecosystem functions. However, in contrast to temperate cites, there is little evidence to show how UHI influences landscape phenological pro- cesses in the tropics. In this study, we examined vegetation phenological responses to urban form, distance from the city centre and surface temperatures, in the tropical city of Kampala, Uganda. Estimates of vegetation growing season length and land surface temperature were derived from MODIS satellite imagery for multiple years (2013 – 2015) and urban form was characterised using the Local Climate Zone (LCZ) classification. We showed that growing season length increased along the urban – rural gradient (p < 0.001) and was longest in the least built-up LCZ class (p < 0.001). Growing season length was significantly reduced as land surface temperature increased (p < 0.001). These findings contrast with results reported for temperate cities, where higher temperatures are often associated with longer vegetation growing seasons. Our findings suggest that enhanced surface temperatures associated with UHI are a limiting factor to season length in the urban tropics. Urban planners in tropical cities should therefore account for vegetation sensitivity to UHI when developing targeted management strategies aiming to optimise the benefits accrued from vegetation.


Introduction
By 2050, 68% of the world's population will reside in urban areas (United Nations. (2019), 2019). Up to 90% of the projected growth of the global urban population will occur in Africa and Asia (United Nations. (2019), 2019), where many cities are already vulnerable to climate change and hazards caused by urbanisation (du Toit et al., 2018). Vegetation is posited to significantly improve urban resilience through the provision of ecosystem functions, such as thermal regulation, yet many cities in the global south are undergoing rapid vegetation loss (Lindley, Pauleit, Yeshitela, Cilliers, & Shackleton, 2018;Yao, Cao, Wang, Zhang, & Wu, 2019). It is therefore imperative to understand the impacts of urbanisation on vegetation to safeguard human and ecosystem health.
One of the most obvious effects of urbanisation is the phenomenon known as the Urban Heat Island effect (UHI), i.e. where urban areas typically experience higher temperatures relative to their surrounding landscapes (Landsberg, 1981;Taha, 1997;Voogt & Oke, 2003). UHIs are caused by high proportions of impervious land cover types that promote heat storage and restrict natural radiative cooling through evapotranspiration. UHI intensities are influenced by the size of cities and their characteristics in terms of urban morphology, biome, regional climate, seasonal changes in meteorological conditions and vegetation cover (Giridharan & Emmanuel, 2018;Imhoff, Zhang, Wolfe, & Bounoua, 2010;Roth, 2007). Urban vegetation mitigates UHI through increased thermal regulation via evapotranspiration (Norton et al., 2015) and the provision of shade (Li, Ratti, & Seiferling, 2018). In recent years, several studies have noted the importance of vegetation for UHI mitigation in tropical cities (e.g. Cavan et al., 2014;Feyisa, Dons, & Meilby, 2014;du Toit et al., 2018;Lindley et al., 2018). However, the specific impacts of urbanisation and UHI on vegetation phenology in the tropics remain less well studied.
Vegetation phenology, defined as the seasonal timing of vegetation growth and reproduction (Fenner, 1998), is essential for primary production and for sustaining many important ecosystem benefits for human populations (Denny et al., 2014). Vegetation responses to urbanisation have been studied extensively for temperate cities where UHIs result in earlier growing season start dates and longer growing season durations compared to surrounding rural areas (Jochner & Menzel, 2015). However, phenological patterns and processes within tropical cities are less well studied. The few urban tropical phenological studies that exist highlight differences in the timing of start of season tree budding compared to temperate cities (Gazal et al., 2008;Jochner, Alves-Eigenheer, Menzel, & Morellato, 2013). The evidence suggests that tree budding is more sensitive to UHI in temperate cities due to the sensitivity of temperate vegetation to springtime temperature increases (and photoperiod) after wintertime dormancy (Zhang, Friedl, & Schaaf, 2006). However, there remains a limited understanding of the effect of UHI on the timing of the end of the growing season and growing season duration in tropical cities.
In tropical natural habitats leaf growth usually occurs in the wet season whereas leaf fall occurs mainly in the dry season (de Camargo, de Carvalho, Alberton, Reys, & Morellato, 2018;Williams, Myers, Muller, Duff, & Eamus, 1997). Given that UHI intensities are most significant during the dry season in tropical cities (Giridharan & Emmanuel, 2018;Roth, 2007), UHIs might have a greater influence on leaf fall than on leaf growth. Moreover, the formation of an Urban Dryness Island (UDI) effect and increased plant water requirements due to UHI-induced potential evapotranspiration (Hao, Huang, Qin, Liu, Li, & Sun, 2018;Luo & Lau, 2019;Wang, Hutyra, Li, & Friedl, 2017;Zipper, Schatz, Kucharik, & Loheide, 2017), may be more pronounced during the dry season. Consequently, the impact of UHI on leaf growth and development may vary depending on the season. To investigate these impacts there is a need for phenological observations that span the entire vegetation growing season.
Satellite remote sensing can provide landscape scale phenological information on the start of the growing season (SOS), the end of the growing season (EOS) and the growing season length or duration (GSL) (Melaas, Wang, Miller, & Friedl, 2016;Zhang et al., 2003). Landscape phenology is inherently different from the phenology of individual species (Badeck et al., 2004), although information about both is useful for understanding phenology in urban environments (Jochner & Menzel, 2015). In temperate cities, vegetation growing in heavily built-up urban areas often has a long GSL due to elevated urban temperatures (Melaas et al., 2016;Zhang, Friedl, Schaaf, Strahler, & Schneider, 2004;Zhou, Zhao, Zhang, & Liu, 2016;Zipper et al., 2016). In contrast, the duration of the vegetation growing season declines along the urban-rural gradient as the degree of urbanisation diminishes (Zhang et al., 2004;Zhou et al., 2016). An equivalent understanding of landscape scale phenological processes in tropical cities is lacking. Whereas temperate regions experience extreme seasonal changes in temperature as the main trigger for vegetation phenology (along with photoperiod), changes in temperature are less drastic in the tropics. Tropical phenology is mainly controlled by rainfall (Clinton, Yu, Fu, He, & Gong, 2014;Zhang et al., 2006) but excess urban heat might act as a limiting (stress) factor in tropical urban contexts.
In this study, we examined the impact of UHI intensities on vegetation phenology in the tropical city of Kampala, Uganda. Our objectives were to: (i) determine the spatial variability in landscape phenology in respect to degree of urbanisation (i.e. urban form and distance from the city centre); (ii) determine the combined effect of urban form and distance from the city centre on Land Surface Temperature (LST); (iii) establish how spatial patterns of LST vary across years; and (iv) assess the effect of variations in LST on phenology.

Study area
The study focussed on the equatorial city of Kampala in East Africa located at 00 • 18 ′ 49 ′′ N 32 • 34 ′ 52 ′′ E. Kampala has rapidly urbanised in recent years and the urban extent of the Kampala Greater Metropolitan Area (KGMA) now covers more than 800 km 2 (Vermeiren, Van Rompaey, Loopmans, Serwajja, & Mukwaya, 2012). Our region of interest (ROI) covers an area extending 20 km from the city centre (approximately 1402 km 2 ) and includes the KGMA. Kampala has a population of over 1.5 million inhabitants, and this is expected to reach 5.5 million by 2030 (United Nations. (2019), 2019). The city has a tropical rainforest equatorial climate (Af) according to the Köppen climate classification. There are two wet seasons (March-May and September-November) and the city has a mean annual precipitation of about 1200 mm. Torrential rains are often observed from March to May and July is normally the driest month.

Data
Urban form, vegetation abundance, vegetation phenology and LST were characterised using remotely sensed satellite imagery for 2013-2015. We constrained our selection of satellite imagery to three years to minimise the effect of rapid changes in urban form across Kampala. The Local Climate Zone (LCZ) classification scheme (Stewart & Oke, 2012) was used to represent urban form, and LCZs were characterised using imagery from the US Geological Survey Earth Explorer Landsat 8 Operational Land Imager (OLI).
The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor was used to estimate spatial and temporal patterns of LST, vegetation abundance and phenology. We obtained vegetation abundance and phenology data using Vegetation Indices as proxy measures of canopy 'greenness' (Huete, Didan, van Leeuwen, Miura, & Glenn, 2011;Senanayake, Welivitiya, & Nadeeka, 2013;Yao et al., 2019;Yuan & Bauer, 2007). The subsections below provide further information on data sources and processing steps, and the overall methodology is summarised in Fig. 1.

Urban form and vegetation abundance
The LCZ framework is a robust, objective and universal approach for characterising urban form for use in climatological research (Kotharkar & Bagade, 2018;Mushore et al., 2019;Stewart & Oke, 2012). There are seventeen classes contained in the LCZ framework, each representing distinct characteristics of surface cover and structure. We used the World Urban Database and Access Portal Tools (WUDAPT) LCZ classification method (Bechtel et al., 2015) for mapping Kampala's LCZs. The WUDAPT method uses supervised machine learning to generate a citywide LCZ map from multispectral Landsat 8 OLI imagery in a 3-step process; namely: (i) acquisition and pre-processing of cloud-free images, (ii) creation of training areas in Google Earth and (iii) implementation of the classification in the open-source GIS software, System for Automated Geoscientific Analyses (SAGA-GIS).
It is recommended that more than one satellite image is used in LCZ classification to minimise the impact of spectral changes in vegetation over time (Bechtel et al., 2015). We selected two Landsat 8 OLI scenes (LC81710602015074 and LC81710602015058), both with low cloud cover. The images were radiometrically calibrated to Top of the Atmosphere reflectance and resampled from 30 m to 100 m. Resampling allows the spectral signals of multiple urban features to be combined and thus facilitates attribution of local climate to zones (Bechtel et al., 2015;Stewart & Oke, 2012).
Training areas with a minimum width of 200 m were digitised within Google Earth for the most common LCZ classes in Kampala: Compact midrise (LCZ2), Compact low-rise (LCZ3), Open midrise (LCZ5), Open low-rise (LCZ6), Large low-rise (LCZ8), Sparsely built (LCZ9), Dense trees (LCZA), Scattered trees (LCZB), Bare soil or sand (LCZF) and Water (LCZG). Subclass LCZ3_F was added to the selection to indicate compact low-rise neighbourhoods (LCZ3) with mostly bare soil surfaces (LCZF), as was done for the city of Nagpur in India (Kotharkar & Bagade, 2018). To minimise the effect of mixed pixels, we digitised polygons in locations with broadly homogenous land cover composition and avoided boundaries between any two given LCZ types. The number of training sample polygons per LCZ class ranged between 3 and 20, and depended on the area covered by each LCZ class and how well areas could be digitised.
The LCZ classification was implemented in SAGA-GIS using its Random Forest classifier, which has high accuracy and computational performance for LCZ classification (Bechtel et al., 2015;Breiman, 2001).
New training areas were iteratively selected, and LCZs reclassified to obtain an LCZ map that compared favourably with existing maps of Kampala's urban morphology. Accuracy assessment of the final LCZ map was performed using the Semi-Automatic Classification QGIS plugin and a total of 764 test polygons (100 m cell size) selected through stratified random sampling (Congedo, 2016). The overall accuracy of the final LCZ map was 73.2% with a Kappa coefficient of 0.674. This is within the expected accuracy range, for instance, within the 60-89% range reported for LCZ maps produced for 20 cities in China (Ren et al., 2019). Four of the eleven LCZ classes (i.e. Compact low-rise and bare soil (LCZ3_F), Open low-rise (LCZ6), Sparsely built (LCZ9) and Scattered trees (LCZB)) covered 87% of the ROI, and all LCZs had producer and user accuracies above 90% and 85% respectively. Given their high accuracy statistics, distinctiveness and large spatial extent, these four LCZ classes (Fig. 2) were selected for further analysis. To validate our choice of LCZ classes, we ascertained that they were different in terms of surface cover using a Kruskal-Wallis H-test performed on Enhanced Vegetation Index (EVI) data. We chose EVI for characterising surface cover (i. e. vegetation abundance) because it minimises canopy background variations and maintains sensitivity over dense vegetation in urban environments (Huete et al., 2011). EVI is an optimised combination of Blue, Red, and Near Infrared (NIR) bands and can be derived using the where ρ is reflectance; L is the canopy background adjustment factor; C 1 and C 2 are aerosol resistance weights. The coefficients are L = 1, C 1 = 6 and C 2 = 7.5 (Huete et al., 2011).
We used the MOD13Q1 EVI product (16-day composites and a spatial resolution of 250 m) acquired from the Oak Ridge National Laboratory (ORNL) Distributed Active Archive Center (DAAC) of the National Aeronautics and Space Agency (NASA). To characterise vegetation abundance across Kampala whilst minimising the effects of seasonal changes in vegetation, we calculated the average of 23 MOD13Q1 EVI images for 2013 using the raster calculator function in ArcGIS 10.0. The vegetation abundance and LCZ maps were converted to polygon vector format and an intersect function of the overlay toolset in ArcGIS 10.0 used to acquire the LCZ classification at each EVI pixel location.

Vegetation phenology
Time series analysis and extraction of vegetation phenology seasonal parameters was performed in the TIMESAT 3.2 program (Eklundh & Jönsson, 2015;Jonsson & Eklundh, 2004) using MOD13Q1 EVI images spanning three years (2013)(2014)(2015) as shown in Fig. 3. The adaptive Savitsky-Golay algorithm which uses local polynomial fitting functions was used to smooth the EVI time series. This approach preserves important features of the dataset (i.e. maximum, minimum and width) whilst reducing the noise due to possible cloud cover (Yao et al., 2017). We set the adaptation strength to 3.0, seasonal parameter to 0.5 (2 seasons per year) and used a Savitzky-Golay window size of 3. We weighted the data for each pixel using the pixel reliability band (i.e. cloud-free assigned 1, marginal data assigned 0.5 and cloudy data assigned 0.1). We extracted the start (SOS), end (EOS) and length (GSL) of the vegetation growing season using the amplitude method. The amplitude method is better suited than the threshold method in urban settings because low vegetation cover in cities makes it hard to find a threshold for determining the start and end of the season (Zhou et al., 2016). The SOS and EOS were the dates that the fitted curve increased and declined to the proportion of the amplitude (set to 20%).

Land surface temperature
We derived LST for 2013-2015 from the MODIS MOD11A2 LST product. The MOD11A2 LST product has a high temporal resolution (8day composite) and is a useful proxy for estimating surface UHI intensity (Yao et al., 2018). We again used the Savitzky-Golay smoothing algorithm and TIMESAT 3.2 to identify the maximum LST for each of the three years as the basis for assessing impacts on vegetation development. As with the EVI data, the MODIS quality assessment images were used to de-noise the data.

Data analysis
Vegetation phenology, vegetation abundance, LST and LCZ data were assigned to each MOD13Q1 EVI pixel as the unit of analysis. Firstly, all data was converted from raster to vector format (conversion toolbox in ArcGIS 10.0) due to differences in the spatial resolution of the input data (vegetation abundance & phenology, LCZ and LST; Fig. 4). The dissolve function (Manage data toolset in ArcGIS 10.0) was applied to the LST and LCZ vector data separately to merge polygons that had similar data and shared boundaries. The intersect function (overlay toolset in ArcGIS 10.0) was then used to acquire LCZ and LST data for each MOD13Q1 EVI pixel. To account for the diminishing effect of degree of urbanisation along the urban-rural gradient, we obtained the distance between each pixel and the city centre using the Near Analysis tool in ArcGIS 10.0. Each data point therefore represented an EVI pixel location with its attributes (i.e. phenology, vegetation abundance, LCZ, distance from the city centre and LST) for each year. As with vegetation abundance, differences in LST and phenology among LCZs were computed using Kruskal-Wallis H-tests.
A series of linear mixed models were used to analyse the effect of urban form (LCZ) and distance from the city centre on SOS, EOS, GSL and LST. We fitted a model for each response variable with LCZ and distance as explanatory variables. Year was used as a fixed effect to control for differences across years in the phenology models. Interactions between year and LCZ and between year and distance from the city centre were included in the LST model to account for varying effects of urban form and distance from the urban core across the three years. Location was used as a random effect to allow for correlated error terms caused by repeated observations (each year) at the same location. The modelling was done using the "lmer" function of the "lme4" package in R (Bates, Mächler, Bolker, & Walker, 2015;R Core Team, 2018). We used a likelihood ratio test to establish the significance of the final full model against a null model comprising the intercept only. The importance of each explanatory variable was determined using a likelihood ratio test (R function "drop1") that compared the full model with reduced models. Variance inflation factors were computed from standard linear models to assess collinearity (Zuur, Ieno, & Elphick, 2010), and no issues were found. Normality and independence of the residuals were confirmed by inspecting QQ plots and plots of residuals against fitted values.
Linear mixed models were also used to assess the influence of LST on SOS, EOS and GSL using LST and Year as fixed effects. This allowed us to examine the influence of LST patterns while controlling for meteorological differences across years. Location was included as a random effect and model diagnostics (significance of final full model and normality of residuals) performed as described in previous models.

Surface cover differences across LCZ types
Vegetation abundance (i.e. mean EVI) differed significantly among the LCZs and between all LCZ pairs (Table 1). EVI decreased from the least to most built-up LCZ type. LCZ9 (Sparsely built) and LCZ6 (Open low-rise) had EVI values that were 17% and 34% lower respectively than those in LCZB. LCZ3_F (Compact low-rise and bare soil) recorded the lowest overall EVI at 47% lower than LCZB.

Phenology
The combined effect of LCZ, distance from the city centre and year had a significant influence on SOS (likelihood ratio test: X 2 = 1097, df = 6, p < 0.001), EOS (likelihood ratio test: X 2 = 7963, df = 6, p < 0.001) and GSL (likelihood ratio test: X 2 = 284.2, df = 11, p < 0.001). LCZ class had a significant influence on SOS, EOS, and GSL (Table 1; Table 2). LCZB (Scattered trees) experienced the earliest SOS and the latest EOS dates resulting in longer growing seasons (Fig. 5). In comparison to LCZB, GSL was shorter in LCZ9 (Sparsely Built) (estimate = − 3.3 days, standard error = 0.5, p < 0.001), LCZ6 (estimate = − 8.6 days, standard error = 0.7, p < 0.001) and LCZ3_F (estimate = − 13.2 days, standard error = 0.7, p < 0.001). Distance from the city centre significantly (positively) influenced GSL but had negligible effects on SOS and EOS   Table 2 Mixed Models testing whether phenology start of season (SOS), end of season (EOS) and growing season length (GSL) vary across years, Local Climate Zone (LCZ) class, distance from city centre and Land Surface Temperature (LST). DenDF and NumDF represent the Numerator and Denominator degrees of freedom. Pr(>F) is the significant coefficient (P) for the F statistic.
LST declined with increasing distance from the city centre across all years, but the magnitude of LST change from urban to rural areas varied across years (interaction between year and distance: F 2, 43819 = 1247.7, P < 0.001). LST ranged between approximately 26-28 • C (rural-urban) in 2014 and 2015 and around 24-28 • C (rural-urban) in 2013 (Fig. 6). The Kruskal-Wallis H-test confirmed significant differences in LST across years (χ 2 = 8575.1, df = 2, P < 0.0001).

Discussion
Our results show that landscape phenology in the tropical city of Kampala is influenced by LST, distance from the city centre and degree of urbanisation (LCZ). During 2013-15, heavily built-up locations experienced high LST, early EOS and a short GSL in comparison to less built-up locations. LCZB (Scattered trees) had the highest vegetation cover and experienced the lowest LST, earliest SOS, latest EOS and longest GSL. LCZ9 (Sparsely built) and LCZ6 (Open low-rise) had GSLs that were 3 and 9 days shorter respectively than in LCZB. LCZ3_F (compact low-rise and bare soil) recorded the shortest GSL overall at 13 days shorter than in LCZB. The order of LCZs in respect to GSL was LCZB > LCZ9 > LCZ6 > LCZ3_F which mirrored the order of LCZs in respect to LST (LCZB < LCZ9 < LCZ6 < LCZ3_F). Furthermore, the decline in LST along the urban-rural gradient varied between years. Rural areas experienced temperatures that were 2 • C cooler in 2013 in comparison to 2014 and 2015.
The LCZs exhibited differences in surface cover characteristics. The relative order of vegetation abundance in each LCZ (i.e. EVI) was LCZB > LCZ9 > LCZ6 > LCZ3_F. This mirrored the order of LCZs in respect to LST (LCZB < LCZ9 < LCZ6 < LCZ3_F). Several studies have attributed lower temperatures within cities to higher vegetation abundance (e.g. Mushore et al., 2019;Senanayake et al., 2013;Yuan & Bauer, 2007) due to enhanced latent heat flux through evapotranspiration (Cavan et al., 2014;Feyisa et al., 2014). On the other hand, a high proportion of impervious land cover enhances thermal admittance and high heat storage resulting in higher temperatures (Landsberg, 1981). In the case of LST, a high amount of vegetation cover is particularly important in terms of its influence on surface albedo and shading (Taha, 1997).  Our findings on the effect of urbanisation on landscape phenology in Kampala contrast with the established phenological patterns in temperate cities, where vegetation experiences earlier SOS, later EOS and longer GSL due to UHI (Melaas et al., 2016;Zhang et al., 2004;Zhou et al., 2016;Zipper et al., 2016). The long growing season observed in highly vegetated locations in Kampala agree with the phenological patterns observed by Guan et al. (2014) in tropical natural habitats. However, the cause and mechanisms determining longer growing seasons were not well established in Guan et al. (2014).
The Urban Dryness Island effect is a phenomenon that describes lower moisture availability in cities compared to the surrounding landscape as a consequence of high proportions of impervious land cover types (Hao et al., 2018;Luo & Lau, 2019;Wang et al., 2017). The spatial patterns of UDI are similar to the thermal gradients associated with UHI (Hao et al., 2018;Luo & Lau, 2019;Wang et al., 2017) and high plant water requirements have been attributed to UHI-induced potential evapotranspiration (Zipper et al., 2017). The UDI in Kampala is therefore also expected to have diminished from the most to least built-up LCZ category (LCZ3_F (Compact low-rise and bare soil) > LCZ6 (Open low-rise) > LCZ9 (Sparsely built) > LCZB (Scattered trees)). Observed phenology patterns in the city are likely to be driven by factors relating to both UHI and UDI.
The spatial pattern of LST in 2013 differed from 2014 and 2015 due to the meteorological anomalies of the drought and El Niño effect in 2014 and 2015 respectively (Zhang, Dannenberg, Hwang, & Song, 2019). Regional temperature anomalies in 2014 and 2015 had strong effects on LST in rural areas, and the effect declined along the rural--urban gradient. These findings suggest that inter-annual differences in regional climate exacerbate the effect of UHI.

Conclusions
In this paper, we provide substantial new evidence about the role of temperature as a limiting factor for GSL in tropical cities. Our study results for Kampala demonstrated that GSLs in LCZ6 (Open low-rise) and LCZ3_F (compact low-rise and bare soil) were 8.6 and 13.2 days shorter respectively compared to LCZB (Scattered trees). Shorter vegetation growing seasons, in turn, limit the provision of beneficial ecosystem functions. For instance, a shorter growing season will influence some aspects of urban agricultural productivity and regulating ecosystem roles such as those associated with human thermal comfort. Urban planning in tropical cities could focus on strategies that aim to mitigate UHI and extend the GSL, for example, through the enhancement of green spaces in highly built-up zones. Knowledge about the environmental processes and intrinsic attributes (surface cover, meteorology and phenology) of LCZ classes could provide useful information to support urban planning in Kampala. For example, city planners could aim to increase vegetation cover in LCZ3_F (compact low-rise and bare soil). Additionally, city planners could restrict further expansion of LCZ3_F in favour of LCZ classes that have higher vegetation cover, lower temperature (and dryness) and therefore longer GSLs (e.g. LCZ6 and LCZ9). These strategies could also be applicable in other tropical cities that are faced with similar urbanisation challenges.
The combined effect of UHI and elevated regional LST might have a more substantial impact on phenology than UHI alone. As the current study focussed solely on local-scale UHI effects on phenology, future work could explore the combined effect of UHI and LST anomalies at a regional scale.