Determining the Influence of Long Term Urban Growth on Surface Urban Heat Islands Using Local Climate Zones and Intensity Analysis Techniques

Urban growth, typified by conversion from natural to built-up impervious surfaces, is known to cause warming and associated adverse impacts. Local climate zones present a standardized technique for evaluating the implications of urban land use and surface changes on temperatures of the overlying atmosphere. In this study, long term changes in local climate zones of the Bulawayo metropolitan city were used to assess the influence of the city’s growth on its thermal characteristics. The zones were mapped using the World Urban Database and Access Portal Tool (WUDAPT) procedure while Landsat data were used to determine temporal changes. Data were divided into 1990 to 2005 and 2005 to 2020 temporal splits and intensity analysis used to characterize transformation patterns at each interval. Results indicated that growth of the built local climate zones (LCZ) in Bulawayo was faster in the 1990 to 2005 interval than the 2005 to 2020. Transition level intensity analysis showed that growth of built local climate zones was more prevalent in areas with water, low plants and dense forest LCZ in both intervals. There was a westward growth of light weight low rise built LCZ category than eastern direction, which could be attributed to high land value in the latter. Low plants land cover type experienced a large expansion of light weight low rise buildings than the compact low rise, water, and open low-rise areas. The reduction of dense forest was mainly linked to active expansion of low plants in the 2005 to 2020 interval, symbolizing increased deforestation and vegetation clearance. In Bulawayo’s growth, areas where built-up LCZs invade vegetation and wetlands have increased anthropogenic warming (i.e., Surface Urban Heat Island intensities) in the city. This study demonstrates the value of LCZs in among others creating a global urban land use land cover database and assessing the influence of urban growth pattern on urban thermal characteristics.


Introduction
Urban areas continue to expand in population and built-up extent, with faster rates in developing countries [1][2][3][4][5][6]. Whereas urban growth is often considered as a sign of economic vitality [3], its adverse impacts that include increased air pollution, Surface Urban Heat Islands, dust and haze, significantly influence urban micro-and macro-climate and affect urban environmental quality and human health [7,8]. Globally, urbanization has caused climate modifications, most evident in higher temperatures in urbanized than the non-urbanized surroundings [9][10][11][12][13]. Such growth often exacerbates heat stress in the already warming (due to global climate change) urban areas, leading to deterioration of outdoor thermal comfort [14]. Urban growth and associated surface changes induce nearsurface warming, which increases energy and water demand due to search for indoor and outdoor thermal comfort [15][16][17][18]. Hence, urban growth assessment techniques that consider both surface and near surface characteristics and atmospheric gas/pollutant emissions are valuable in determining the influence of anthropogenic processes on the urban thermal emissions to inform sustainable urban growth.
Remote sensing offers a variety of data for analyzing spatial and temporal effects of urban land surface changes on the urban thermal environment. This has largely been facilitated by advances in sensor development that has improved the quality and availability of remotely sensed data. For instance, missions such as Landsat offer large archival data spanning as far back as 1972 at reasonable and improved radiometric, spectral and spatial resolution, valuable for assessing both large scale and localized temporal and multitemporal landscape and environmental patterns [19][20][21][22][23]. A large body of literature has investigated the impact of land use land cover (LULCs) changes on surface and near surface temperature, e.g., Kumar and Shekar [24] and Uddin et al. [25]. These studies have indicated that impervious, bare and built-up surfaces result in urban warming while vegetation and water bodies act as thermal sinks. However, whereas LULC-based techniques and surface characterization schemes account for the contributions of land surface changes to temperature dynamics, they often exclude other anthropogenic contributions such as emission, which are key drivers of temperature changes associated with urban growth. Hence, schemes that account for both land surface characteristics and gas emissions/pollutants as drivers of changes in the thermal environment are necessary to adequately explain climatic changes in such complex environments.
Due to the dependence on LULC schemes, studies on effects of urbanization on temperature have mostly defined Surface Urban Heat Island (SUHI) as the difference between "urban" and "rural" temperature [26][27][28][29][30][31]. In such studies, rural and urban are vaguely defined by differences in population and built-up extent in a manner that is not universal [32][33][34]. However, this separation is no longer always clear cut as traditional and non-traditional urban and rural land uses increasingly continue to coexist [33]. The traditional classification scheme is area specific, making it difficult to make global comparisons as LULC characteristics vary between cities of the same country and between countries. However, analysis based on Local Climate Zones (LCZ) provide understanding of urban structures and land uses in a globally standardized manner, useful for understanding the influence of urbanization on urban climate [33][34][35]. The LCZ scheme is local and climatic in nature considering surface cover, three dimensional surface structures (such as height and density of buildings and vegetation) as well as anthropogenic thermal emissions [8,[35][36][37][38][39]. LCZs provide useful information for assessing adherence of cities to the 2030 agenda for sustainable development Goal 11 [40] both in the form of LULC transitions and anthropogenic emissions effects on climate. Hence, LCZ provide a standard basis upon which the impacts of urban growth on the thermal environment can be monitored and assessed.
The LCZ scheme emphasizes the difference in temperature among the categories within and between cities, thus directly linked to climate of an area, while contributing to the creation of the global urban database [11,32]. Close association between LCZ and LST shows that LCZs are helpful for examination of evolution of SUHI over time [41]. LCZs are more conducive to analysis and less prone to confusion because they highlight common exposure characteristics and invite physically based explanations of SUHI magnitude [42][43][44][45][46][47][48][49]. Studies which used LCZs in SUHI analysis mostly focused on short temporal scales such as diurnal, seasonal and annual, e.g., United Nations General Assembly [40] and Ardiyansyah et al. [49]. Focus on long term interactions between LCZs and SUHI have remained understudied, especially in Africa. Furthermore, although Stewart and Oke [33] showed the effectiveness of LCZs in defining SUHI in cities, application of inter-LCZ temperature difference to define SUHI has remained minimal, especially in the analysis of long-term changes. Most of the studies that analyzed the relationship between LCZs and SUHI either used LST to directly represent SUHI [33,46,49], reclassified LST into different SUHI categories [45], converted LST into other forms such as Distribution Index [41] or normalized temperatures [47,50], or used urban to rural temperature gradient [29,48,51,52].
Of the few studies that used temperature difference between LCZ types to quantify SUHI, most of them used the Low plants LCZ as a reference against which LSTs of other LCZs were compared [50,[53][54][55]. A number of studies on LCZ have mainly focused on the current state and short-term variations as well as contributing to World Urban Database and Access Portal Tool (WUDAPT), e.g., Cai et al. [32], Cai et al. [36], Danylo et al. [37], Qiu et al. [39], and Demuzere et al. [56], with little effort towards understanding and explaining their long term changes that affect a city's SUHI intensity. In long term analysis, single date imageries are commonly used to develop static LCZ for each year, ignoring the value of combining multi-seasonal data for the same analysis. Long term changes provide a better understanding of the contribution of human activities to local climate change and assessment of adherence of city growth patterns to Agenda 2030 for Sustainable Development Goal 11 [39]. Hence, there is a need to use multi-season imageries to generate LCZ and LCZ-based SUHI to determine long term effects of urban growth on a city's thermal environment.
Cities in developing countries have been experiencing growth characterized by massive changes in land surface characteristics and intensification of activities, which have the potential to pollute the atmosphere and exacerbate changes in local thermal environments. Although the trends have been observed in different parts of the world, actual changes vary between and within countries, triggering the need for detailed and city-specific assessments. Hence, in-depth understanding of a city's specific influences on LCZ is important for ensuring that further development is climate smart and sustainable. However, available literature on long term LCZ changes [52,[57][58][59] uses the traditional "from to" change detection approach which lacks in depth analysis to provide detailed understanding of long term LCZ transitions and their potential impacts on local climate. Although not yet applied to understand long term LCZ transitions, in depth analysis of changes based on transition matrices of different periods is better done using intensity analysis than the traditional "from to" change detection approach. Intensity analysis is useful for effecting classifications of different time intervals to understand sizes and intensities of temporal changes among categories [60][61][62][63], as it provides details on whether transition from one category to another deviates from a uniform process [61]. It also identifies time intervals when rate of change was fast or slow, identifies whether category changes were active or dormant in a time interval and whether a category was targeted or avoided by changes during an interval [62,64,65]. Furthermore, it analyzes land changes relative to size of category to identify systematic transitions over time [62]. As such, it reveals information such as underlying processes associated with changes which ordinary change detection conceals [66]. For instance, Alo and Pontius [67] revealed that protected areas in Ghana experienced systematic transitions from closed forest to bare and bush fires, while Ekumah et al. [66] revealed that between 1985 and 2017, human induced LULCs grew at the expense of natural categories in the Densu Delta, Sukumo II and Muni Pomodze Ramssar sites in Ghana. Intensity analysis is thus valuable for obtaining an in-depth and detailed analysis and understanding of long term LCZ changes, especially for cities such as Bulawayo where spatial and temporal temperature patterns are not yet documented. Combining intensity analysis and SUHI retrieval in the context of LCZ will therefore provide a detailed understanding of the effect of urban growth patterns on the thermal environment of cities.
Hence, the aim of this study was to integrate intensity analysis and LCZ-based SUHI retrieval to provide detailed analysis of the impact of urban growth on the thermal environment in Bulawayo metropolitan city in Zimbabwe. Specifically, this study sought to determine long term effect of urban growth on the thermal environment using LCZ between 1990 and 2020 in Bulawayo city, Zimbabwe. The study also utilized intensity analysis for an in-depth assessment of the changes in LCZ in Bulawayo between 1990 and 2020. Additionally, and contrary to the broad literature that uses the between rural-urban difference in temperature, this study enhanced the use of inter-LCZ temperature difference approach to quantify SUHI intensity and their long term changes.

Description of the Study Area
Bulawayo is the second largest city in Zimbabwe ( Figure 1). It is located to the southeast of the country (Figure 1a) at an elevation of approximately 1358 m above sea level. The period between October and March is hot and wet with a minimum of 16 • C and maximum of 30 • C, with an average temperature of 25 • C, while the rest of the year is cool and dry, with a minimum of 10 • C, maximum of 25 • C and average temperature of 15 • C [68]. Generally, the area receives erratic rainfall, with annual average precipitation of 600 mm that ranges from 199.3 mm to 1258.8 mm, typical of a semi-arid climate. Bulawayo (Figure 1b) lies in the subtropical steppe (Bsh) according to Koppen climate classification. The period from December to February is the wettest. Most rain falls from December to February and the area is vulnerable to droughts due to proximity to the Kalahari Desert [69,70]. 2020. Additionally, and contrary to the broad literature that uses the between rural-urban difference in temperature, this study enhanced the use of inter-LCZ temperature difference approach to quantify SUHI intensity and their long term changes.

Description of the Study Area
Bulawayo is the second largest city in Zimbabwe ( Figure 1). It is located to the southeast of the country (Figure 1a) at an elevation of approximately 1358 m above sea level. The period between October and March is hot and wet with a minimum of 16 °C and maximum of 30 °C, with an average temperature of 25 °C, while the rest of the year is cool and dry, with a minimum of 10 °C, maximum of 25 °C and average temperature of 15 °C [68]. Generally, the area receives erratic rainfall, with annual average precipitation of 600 mm that ranges from 199.3 mm to 1258.8 mm, typical of a semi-arid climate. Bulawayo (Figure 1b) lies in the subtropical steppe (Bsh) according to Koppen climate classification. The period from December to February is the wettest. Most rain falls from December to February and the area is vulnerable to droughts due to proximity to the Kalahari Desert [69,70].

Field Observations of Local Climate Zones in Bulawayo
Since the WUDAPT places little emphasis on field data collection, a survey to identify and obtain ground truth samples of LCZ categories in the study area was important to guide digitizing of training polygons in Google Earth Field observations, which allowed for identification of inter-and intra-category variabilities that could not be adequately captured from Google Earth. The sample coordinates of each LCZ category (ground truth data) were obtained between 18 and 27 October in 2020. This experience also guided selection of training areas for the historical periods using Google Earth in the absence of field measurements for that period. Field observations increased the validity of the analysis instead of exclusive reliance on Google Earth retrievals. Generally, 8 LCZ categories were identified in the study area that fit into the description of LCZs provided by Stewart and Oke [33]. The categories were three land use-based LCZs, namely Compact low rise (LCZ3), Open low rise (LCZ6) and Light weight low rise (LCZ7), as well as three land cover-based LCZs, which were Dense forest (LCZA), Low plants (LCZD) and Water (LCZ

Field Observations of Local Climate Zones in Bulawayo
Since the WUDAPT places little emphasis on field data collection, a survey to identify and obtain ground truth samples of LCZ categories in the study area was important to guide digitizing of training polygons in Google Earth Field observations, which allowed for identification of inter-and intra-category variabilities that could not be adequately captured from Google Earth. The sample coordinates of each LCZ category (ground truth data) were obtained between 18 and 27 October in 2020. This experience also guided selection of training areas for the historical periods using Google Earth in the absence of field measurements for that period. Field observations increased the validity of the analysis instead of exclusive reliance on Google Earth retrievals. Generally, 8 LCZ categories were identified in the study area that fit into the description of LCZs provided by Stewart and Oke [33]. The categories were three land use-based LCZs, namely Compact low rise (LCZ3), Open low rise (LCZ6) and Light weight low rise (LCZ7), as well as three land cover-based LCZs, which were Dense forest (LCZA), Low plants (LCZD) and Water (LCZ G). The study used LCZs definitions and pictorials provided by Stewart and Oke [33] as reference to identify similar classes in the study area for global comparability.

Multi-Temporal Remotely Sensed Datasets
Multi-temporal Landsat 5, Landsat 7 and Landsat 8 Operational Land Imager datasets were downloaded from the United States Geological Survey's (USGS) earth explorer website for analysis. Cloud free imageries for the wet and dry vegetation periods were downloaded to minimize compromising effects of atmospheric noise on image radiometric values and LCZ mapping accuracy. Table 1 shows the imageries used for 1990, 2005 and 2020. In this study, Landsat thermal, panchromatic as well as bands for monitoring cirrus clouds and coastal aerosols were not used for analysis. For each year, dry and wet periods were selected in order to capture seasonal variations in LCZ, especially in areas with vegetation. The post rain period was chosen to represent the wet biomass period because during that period, trees and grasses are vibrant after a rainy season. The rainy season was avoided to attain both temporal and multi-temporal cloud free imagery. Amorim [71] indicated that Surface Urban Heat Island intensity is influenced by responses of vegetation to rainfall patterns before the imagery date. Amorim [71] observed that Heat Island Intensities increased during periods of high biomass, which reduce temperatures in vegetation areas. Therefore, precipitation patterns of up to 10 days prior to overpass (amounts and rainy days) are shown in Table 1. Generally, the number of rainy days prior to overpass was higher in the post-rain than other seasons, with the lowest number in the cool season. Cumulative rainfall amounts in 10 days before overpass were also low in all seasons (less than 20 mm).

Mapping of LCZ Using Dry and Wet Season Imagery
The advantage of the LCZ scheme is that their mapping follows an easy and standardized approach for mapping LCZ, which involves downloading of imagery, digitizing of training sites (for classification and accuracy assessment) on Google Earth and supervised classification using the random forest (RF) classifier [37,72,73]. Local Climate Zones for Bulawayo were thus mapped following the WUDAPT L0 procedure [8,36,56,73,74]. The procedure involves downloading suitable imagery of the study area, on-screen selection and digitizing of training areas on Google Earth, supervised image classification using the Random Forest (RF) Classifier and post classification accuracy assessment in SAGA GIS. The procedure was adopted due to its use of readily available and freely downloadable imageries as well as easy access to the SAGA GIS software for implementation. Additionally, the steps have been followed in different parts of the world, making it very easy to follow. Use of the RF classifier makes the procedure attractive as it can perform bootstrapping analysis which is used to assess quality of the LCZ database [56,75]. The RF model is a collection of decision trees and each tree is made up of a subset of training dataset for a subset of predictors [73,76,77]. It is termed RF because its subsets are randomly formed. In RF classification, the predicted value is the mode of the predictions from all trees. Main advantages of RF are the use of both categorical and numerical values, the evaluation of the precision of prediction, the robustness in the presence of outliers, noise, and overfitting. The RF model can quantify the contribution of each predictor to the total spatial variability of the target and assigns a variable importance score to each predictor. Random forest requires a small amount of training data, yet provides competitive results and can handle a large volume of input data without deletion, while still capable of identifying important variables for classification [78][79][80][81][82][83]. In addition, RF is not sensitive to over-training or noise and is desirable for multi-source remote sensing and geographical information systems data [75,81]. LCZ maps were produced for the years 1990, 2005 and 2020 using the same training areas. These were collected from locations whose LULC category did not change over time in order to eliminate the effect of differences in ground truth data on mapping accuracy. The use of the same training areas was made possible by availability of historical Google Earth imagery where the areas could be clearly identified at different periods. In order to adhere to the definition of LCZ, which requires that they cover at least a hundred meters to several kilometers to influence temperature [8,11,41], the maps were resampled using a 5 by 5 pixels window. A LCZ must be large enough to influence temperature of an area. In order to quantify the effect of seasonality on mapping accuracy, LCZ maps were also produced using data for the hot and dry season for comparison with analysis based on a combination of data from the post rain, cool dry and hot dry periods. Figure 2 provides a summary of the procedure followed. The broken arrow in the figure shows the approach taken by previous studies which largely skip the iterative step of qualitative accuracy assessment, further improving accuracy before change detection. [56,75]. The RF model is a collection of decision trees and each tree is made up of a subset of training dataset for a subset of predictors [73,76,77]. It is termed RF because its subsets are randomly formed. In RF classification, the predicted value is the mode of the predictions from all trees. Main advantages of RF are the use of both categorical and numerical values, the evaluation of the precision of prediction, the robustness in the presence of outliers, noise, and overfitting. The RF model can quantify the contribution of each predictor to the total spatial variability of the target and assigns a variable importance score to each predictor. Random forest requires a small amount of training data, yet provides competitive results and can handle a large volume of input data without deletion, while still capable of identifying important variables for classification [78][79][80][81][82][83]. In addition, RF is not sensitive to over-training or noise and is desirable for multi-source remote sensing and geographical information systems data [75,81]. LCZ maps were produced for the years 1990, 2005 and 2020 using the same training areas. These were collected from locations whose LULC category did not change over time in order to eliminate the effect of differences in ground truth data on mapping accuracy. The use of the same training areas was made possible by availability of historical Google Earth imagery where the areas could be clearly identified at different periods. In order to adhere to the definition of LCZ, which requires that they cover at least a hundred meters to several kilometers to influence temperature [8,11,41], the maps were resampled using a 5 by 5 pixels window. A LCZ must be large enough to influence temperature of an area. In order to quantify the effect of seasonality on mapping accuracy, LCZ maps were also produced using data for the hot and dry season for comparison with analysis based on a combination of data from the post rain, cool dry and hot dry periods. Figure 2 provides a summary of the procedure followed. The broken arrow in the figure shows the approach taken by previous studies which largely skip the iterative step of qualitative accuracy assessment, further improving accuracy before change detection.

Accuracy Assessment
The WUDAPT procedure automatically splits training areas into 50% for classification and 50% for accuracy assessment. A confusion matrix is formulated in a tabular form for the purposes of comparing reference class labels (ground truth) with labels for corresponding pixels on a derived classification of remote sensing imagery [83]. The diagonal values on the matrix indicate where categories assigned on the classified map correctly corresponded with ground truth. For example, the matrix shows the number of pixels which were assigned same LCZ value as observed on the ground as well as those that are misallocated following classification of remote sensing data. The confusion matrix was used to generate indicators of accuracy at class level (Producer Accuracy and User Accuracy) as well as at entire study area level (Overall Accuracy-OA-and Kappa-K). The use of OA and K with the inclusion of ground data is the most common and reliable way of assessing accuracy [84]. Accuracy assessment is important for the separation of real changes from changes due to errors for LCZ maps of different periods. Error analysis was done for all the study years. Accuracy was also assessed qualitatively by overlaying the produced LCZ maps with a corresponding Google Earth imagery in combination with expert judgement. The number of training areas was objectively and iteratively increased in areas where marked mismatches were observed to capture for intra-and inter-class variabilities using Google Earth.

Detection of Long-Term Changes in LCZ in Bulawayo
A 30-year period (1990 to 2020) was selected, as the study aimed at using LCZ dynamics as a proxy for temperature changes in Bulawayo. The World Meteorology Organization (WMO) recommends a minimum of 30 years for a representative climate change analysis. The period was further split into two 15 year periods (i.e., 1990 to 2005 and 2005 to 2020) for understanding of rapid changes that occur at local scale and for intensity analysis purposes. Additionally, Coppin and Bauer [85] recommended a period of at least 3 years for change detection involving forests and other vegetation types. In an urban setting, there is a mix of rapid and slow LCZ making a period of at least 15 years enough to detect effect of all change trajectories on the climate of an area. A post classification change detection approach was used. The approach produces a change matrix/table which shows the number of pixels which were converted to other LCZ types or remained in the same categories in the considered interval. For instance, the table indicated the number of pixels which were in LCZ1 at the beginning and remained in the same as well as those that were changed to other LCZ categories during the same time interval. Although the change matrix is useful in depicting the quantities and directions of change, it does not adequately explain the changes [60,85,86], hence the need for intensity analysis.

Intensity Analysis for In-Depth Characterization of LCZ Changes
Intensity analysis was used to obtain a better understanding of LCZ transitions between 1990 and 2020 in Bulawayo. It was used to assess locations and intensities of temporal changes among categories. The analysis was done at the interval, category and transition levels [87][88][89][90] using freely available software on Pontius Clarke University web page (https://www2.clarku.edu/faculty/rpontius/ (accessed on 15 April 2021)). The site provides an easy to use Excel sheet where the change matrix for a given time interval is entered. Varga et al. [91] provide descriptions and defining equations used in the analysis that were adopted in this study. S t is the uniform intensity in interval t. An interval change is fast if it exceeds uniform intensity and slow if otherwise.

Category Level Analysis
Changes in gross loss or gain in intensity among different categories was described by category level intensity analysis during time interval t (where t separately represents the two interval periods-1990 to 2005 and 2005 to 2020). Category level loss (L ti ) and gain (G tj ) during the interval t are obtained using Equations (2) and (3) as According to [64], calculation of category persistence (P ti ) is done using Equation (4): S t = L ti = G tj for all categories i and j if changes are uniformly distributed across spatial extent. Uniform transition assumes category i uniformly changes to other categories during the time interval. If L it > S t then the loss of category i is active in the interval t while loss is dormant if L it < S t . Dormant implies that the loss of category i slowed down or stopped within the interval t. Similarly, gain of category j in the interval t is active if G tj > S t and dormant if G tj > S t . When two intervals are considered and the status of a change as dormant or active is same in both intervals, then the category's loss or gain is said to be stationary.

Transition Level Analysis
The analysis describes the variation in intensity with which the gain of a particular category transitions from other categories within a time interval [60]. Transition level intensity analysis was used to determine whether change of a category avoids or targets other categories. If the intensity of the change from category i to category j exceeds uniform intensity, then category i targets j otherwise it avoids.

Retrieval of Changes in SUHI in Response to Long Term LCZ Dynamics
Thermal data of Landsat 5, 7 and 8 were used to compute land surface temperature (LST) for 1990, 2005 and 2020. Initially, the data were corrected of differences in solar zenith angles. While Landsat 8 has two thermal infrared bands, a single channel technique was applied for all the periods to minimize effects of differences in computation algorithms on LST variations between time periods. Digital numbers of thermal data were converted to radiances, which were then used to determine brightness temperature (T b ) and surface temperature (T s ) using Equations (4) and (6), respectively [91][92][93][94]. correction was applied on brightness temperature to obtain actual land surface temperature using Equation (6) [96].
where λ is the central wavelength of emitted thermal radiance (11.5 µm for Landsat 5 and Landsat 7 and 10.9 µm for band 10 of Landsat 8) and ρ is equal to 1.438 × 10 −2 mK. The procedure above was used to retrieve LST on two different dates so that independent sets were used for training and accuracy assessment of the developed estimation algorithm. The spatial structure of LST intensities were used for visual and quantitative analysis of changes which occurred between 1990 and 2020. Stewart and Oke [33] defined SUHI as the difference in LST between LCZs. In this study, we adopted an approach by Dimitrov [97], which defined SUHI as the maximum temperature difference between LCZs. For each year, the Surface Urban Heat Island was computed as the difference between the average T s of LCZ category and the mean surface temperature of the water LCZ (LCZ G), which was identified as giving the highest LST difference with other LCZs consistently in all years using Equation (7).
SUHI LCZ is the SUHI derived from LST difference between other LCZs (X) and the Water LCZ (Y). Although studies such as Stewart and Oke [33] used LCZ D as reference, it was not applicable in this study due to the varying and opposing effects of the LCZ in different places and seasons for daytime analysis. For instance, agricultural areas had thermal values that vary in space and time and between seasons, rendering them as heat sinks in some instances and heat sources in others. As such, during the growing season, they acted as heat sinks while in the dry season, their heat mitigation value was reduced. In other areas, they were completely removed, as they turned to dry biomass or bare soil areas. Similarly, grasslands (also in LCZ D) vary in heat mitigation value depending on season and maintenance, making them another example of inconsistency of LCZ D. The LCZ G was chosen as a reference since it was the coolest and more stable than the vegetation based categories, which have broad seasonal and long term variations in characteristics. Mean LST per LCZ category was obtained using the Zonal Statistics overlay function in ArcGIS version 10.2 for each year. Changes in SUHI per LCZ strata were monitored and linked with observed changes in LCZ from 1990 to 2020 in 15-year intervals.

LCZ Maps Based on Multi-Seasonal Image Analysis
The use of imagery for the wet and dry vegetation periods reduced the confusion between light weight low rise and low plants in the western areas ( Figure 3). The overall classification accuracies were 98%, 98.2% and 95%, for 1990, 2005 and 2020, respectively. Visual inspection shows that between 1990 and 2020, light weight low rise LCZ was spread westwards into areas formerly occupied by low plants. All built local climate zones increased in spatial coverage while low plants and dense forest LCZ decreased in coverage. The water LCZ also decreased in coverage during the study period.     Table 3 shows that Compact low rise increased by 1.5% between 1990 and 2005 and a further 2.2% between 2005 and 2020, giving a 30-year expansion of 3.7%. A significant decrease in coverage was observed in the dense forest LCZ, which experienced a net reduction of 42% between 1990 and 2020, with greater change in the 1990 to 2005 than 2005 to 2020 periods. Generally, all built-up LCZ increased in coverage between 1990 and 2020, with larger expansion in the light weight low rise than other built-up LCZs. Sustained contraction was recorded in dense forest and water LCZs. The low plants LCZ, which in   Table 3 shows that Compact low rise increased by 1.5% between 1990 and 2005 and a further 2.2% between 2005 and 2020, giving a 30-year expansion of 3.7%. A significant decrease in coverage was observed in the dense forest LCZ, which experienced a net reduction of 42% between 1990 and 2020, with greater change in the 1990 to 2005 than 2005 to 2020 periods. Generally, all built-up LCZ increased in coverage between 1990 and 2020, with larger expansion in the light weight low rise than other built-up LCZs. Sustained contraction was recorded in dense forest and water LCZs. The low plants LCZ, which in this study included croplands, grasslands and parks increased by 4.3% over the entire period, except a 0.9% decrease recorded between 2005 and 2020. Low plants were active losers in the 1990 to 2005 interval and became dormant losers (the loss of the class was slow or stopped along the interval) in the 2005 to 2020. The water LCZ was an active loser in both intervals. The light weight low rise LCZ was a dormant gainer in the initial interval, but turned into an active gainer in the interval 2005 to 2020.      The gain of the dense forest LCZ in the 1990 to 2005 interval targeted low plants (Figure 7a). This could imply growth of trees in grasslands such as parks in addition to other tree planting efforts. The expansion of low plants LCZ targeted dense forest, open low rise and light weight low rise, while it avoided water and compact low rise in the 2005 to 2020 interval (Figure 7b). The study also noticed slight spectral confusion between light weight low rise and low plants, especially in the western parts of the study area.  The gain of the dense forest LCZ in the 1990 to 2005 interval targeted low plants (Figure 7a). This could imply growth of trees in grasslands such as parks in addition to other tree planting efforts. The expansion of low plants LCZ targeted dense forest, open low rise and light weight low rise, while it avoided water and compact low rise in the 2005 to 2020 interval (Figure 7b). The study also noticed slight spectral confusion between light weight low rise and low plants, especially in the western parts of the study area.  The gain of open low rise LCZ in the 1990 to 2005 interval targeted low plants (Figure 8a), while in the interval 2005 to 2020 it targeted low plants, water, dense forest and compact low rise (Figure 8b). The gain targeted dense forest more than low plants LCZ, which could be due to the nature of the LCZ that consists of a few and well-spaced buildings surrounded by trees and grass. Expansion into water reveals an adverse environmental impact, with growth intrusion into wetland areas. The gain of open low rise LCZ in the 1990 to 2005 interval targeted low plants ( Figure  8a), while in the interval 2005 to 2020 it targeted low plants, water, dense forest and compact low rise (Figure 8b). The gain targeted dense forest more than low plants LCZ, which could be due to the nature of the LCZ that consists of a few and well-spaced buildings surrounded by trees and grass. Expansion into water reveals an adverse environmental impact, with growth intrusion into wetland areas.    In each year, SUHI intensities between the Compact low rise and Low plants LCZs were comparable, although in 2020 the Low plants were slightly warmer ( Figure 10). SUHI intensities increased with built-up density as well as density of tall buildings evidenced by largest intensity in Compact low rise (e.g., 11.7 °C in 2020) and lowest in Open low rise In each year, SUHI intensities between the Compact low rise and Low plants LCZs were comparable, although in 2020 the Low plants were slightly warmer ( Figure 10). SUHI intensities increased with built-up density as well as density of tall buildings evidenced by largest intensity in Compact low rise (e.g., 11.7 • C in 2020) and lowest in Open low rise (10 • C). Dense vegetation LCZ areas were cooler than built LCZ in all the periods. In each year, SUHI intensities between the Compact low rise and Low plants LCZs were comparable, although in 2020 the Low plants were slightly warmer ( Figure 10). SUHI intensities increased with built-up density as well as density of tall buildings evidenced by largest intensity in Compact low rise (e.g., 11.7 °C in 2020) and lowest in Open low rise (10 °C). Dense vegetation LCZ areas were cooler than built LCZ in all the periods.

Discussion
Mapping accuracy was very high, exceeding 95% for 1990, 2005 and 2020 LCZ maps. Accuracy of at least 95% in LCZ mapping were also recorded by Danylo et al. [37] in Kyiv and Lviv in Ukraine for city specific analysis, which was higher than the below 75% accuracy they achieved using training data for multiple city mapping. The differences in accuracy between city level and large scale LCZ mapping approaches demonstrate that better LCZ maps are generated with city specific efforts than when training areas of another city are used for LCZ mapping [55]. The high accuracy achieved stresses the high performance of random forest classifier in comparison to other classifiers such as support vector machine [75,80,81,83]. RF uses a set of classifiers that make it superior to individual classifier based approaches [81]. The use of multi-seasonal data also enhanced discriminability of LCZs in this study, which resulted in high yearly accuracies. LCZ maps did not show most of the linear features such as roads and rivers due to the large filter used. As noted by Kotharkar and Bagade [8] and Gal et al. [72], an LCZ must be large enough to affect an area's temperature, while use of a large filter removes linear features and expands urban area. According to Kotharkar and Bagade [8] and Gal et al. [72], shifting to a coarser scale results in loss of a number of LCZs.
Consistent with global trends, there was a general increase in built LCZ in Bulawayo between 1990 and 2020. This is a characteristic of city growth globally typified by conversion from natural land covers to urban fabric. For instance, in Kigali, Rwanda, Akinyemi et al. [60]  This may signify growth of the industrial and CBD to serve the surrounding residential areas, which were also expanding. The intensity of transition of water to compact low rise was greater in the interval 1990 to 2005 than 2005 to 2020. In both intervals, the gain of compact low rise also avoided low plants and dense forest LCZ areas. This may signify adherence to the Environmental Management Act, which was signed into law in 2013 [98]. The Act has increased protection of wetlands and the general environment with non-compliance, especially of business enterprises attracting heavy fines. Avoidance of water areas by expansion of compact low rise may also be a result of costs associated with construction in wetlands.
The The gain of the dense forest LCZ in the 1990 to 2005 interval which targeted low plants could imply growth of trees in grasslands such as parks in addition to other tree planting efforts. Over a period of 15 years, the planted and naturally growing trees can increase in canopy size, density and leaf area, enough to be separable from low plants. The expansion of low plants LCZ targeted dense forest, open low rise and light weight low rise, while it avoided water and compact low rise in the 2005 to 2020 interval. The spread into dense forests indicates deforestation, which turns formerly dense forest LCZ into low plants areas such as grasslands and croplands. As the light weight low rise occupied formerly low plant areas, demand for areas such as peri-urban agriculture increased. This may cause communities to spread activities into unused areas by clearing some of the dense forests. Low plants targeting light weight low rise areas could be associated with growth of vegetation within the built-up LCZ. The vegetation includes edible vegetation in small gardens as well as lawns around households. The study also noticed slight spectral confusion between light weight low rise and low plants, especially in the western parts of the study area. However, most of the confusion was eliminated through the use of multi-season data, which increased inter-class discriminability in supervised classification.
The gain of open low rise LCZ in the 1990 to 2005 interval targeted low plants while in the 2005 to 2020 interval it targeted low plants, water, dense forest and compact low rise. The gain targeted dense forest more than low plants LCZ, which could be due to the nature of the LCZ that consists of a few well-spaced buildings surrounded by trees and grass. Expansion into water demonstrates adverse environmental impact with growth intrusion into wetlands. The expansion of open low rise into low plant areas could indicate development of spacious settlements into formerly grassland, bare and agricultural areas. Furthermore, the expansion into other LCZs could indicate the advantage of wealth, as this LCZ is mostly occupied by medium to high income strata which can afford land and develop in any area. The expansion of this spacious LCZ could also be a sign of economic emancipation which enables residence of the city to purchase tracts of large land. In Zimbabwe, this includes existing land owners in densely built-up low income areas but who prefer accommodation in low density built-up areas. Due to increased demand for such spacious settings, developments encroach into formerly protected LCZs such as wetlands and dense forests.
Intensity analysis provided details of LCZ transitions in Bulawayo between 1990 and 2020 beyond the usual "from to" change detection analysis. It revealed the tendency of built-up growth, which was at the expense of vegetation areas, especially low plants. The analysis also showed that the growth of light weight low rise targeted low plant areas and completely avoided compact low rise and open low rise areas. This could be associated with the cost of land in the compact low rise and open low rise, mainly found in the central business district and low-medium density residential (largely occupied by medium to high income strata), respectively. The findings of this study emphasize the argument by Niya et al. [63] that intensity analysis clarifies substantial causes and processes of land use changes. Additionally, in agreement with Huang et al. [66], intensity analysis can assess evidence of a particular change and help develop hypothesis concerning processes of change.
High temperature surfaces expanded while LST temperatures increased between 1990 and 2015. According to Nayak and Mandal [99], urbanization causes temperature change due to both alteration of land use land cover and greenhouse gas concentrations. Similarly, in this study, we attributed surface warming to both LCZ transitions and background cause by anthropogenic activities such as industrial emissions. Blake et al. [100] also reported that Harare was warming, despite cooling in the decade from 1900 to 2002. The expansion of high LST areas and increase in LST intensities was largely due to replacement of natural land covers with built-up LCZs. Buildings and impervious surfaces have high heat absorption capacities which cause elevation of LSTs, especially where vegetation fraction and surface wetness are low. High LSTs were recorded in compact low rise areas especially as their coverage increased over time. This is in agreement with the sentiment that large areas of densely packed buildings create homogeneous areas with high LST [49]. Therefore, the growth patterns of Bulawayo have caused warming due to massive replacement of natural surfaces with buildings and impervious surfaces.
Each year, SUHI intensities between the Compact low rise and Low plants LCZs were comparable, although in 2020, the Low plants were slightly warmer. This is because field observations showed that during the hot and dry seasons, low plant LCZ areas are characterized by dry vegetation and even bare or close to bare ground in areas that are cleared at the end of farming seasons. On the other hand, the link between heat stress and built-up extent is also complex, as buildings also provide shading while on the other hand reducing ventilation leading to opposite effects on thermal comfort [100,101]. Additionally, vegetation within built-up LCZs has a heat mitigation effect, which can reduce differences in SUHI intensities between built-up and natural (land cover based) LCZs.
Surface Urban Heat Island intensities increased with built-up density as well as density of tall buildings evidenced by largest intensity in Compact low rise. According to Stewart and Oke [33], built-up LCZs vary in air temperature depending on factors that include density and height of buildings as well as type and density of vegetation within built-up areas. Consistent with this study, Lelovics et al. [11] stressed that LCZ2 is warmer than LCZ3, which is warmer than LCZ6. They also concluded that LCZ maps can distinguish areas based on degree of LULC modification. Contrasts in temperature between classes with differences in geometry/cover can exceed 10 • C, while classes with few physical differences can be less than 2 • C [33]. Similarly, Lau et al. [14] recorded the highest temperature (38.9 • C) and lowest (29.9 • C) in land cover LCZs in Hong Kong. According to Qiu et al. [39], LCZ scheme considers three-dimensional surface structure and anthropogenic parameters such as heat from human activities that influence temperature. Based on this understanding, comparatively very high land surface temperatures and SUHI were observed in the compact low rise area of the city, followed by the densely packed light weight low rise while the open low rise LCZ was the coolest of the built LCZ.
Dense vegetation LCZ areas were cooler and had lower SUHI intensities than built LCZ in all the periods. Even as vegetation cover declined, their heat mitigation value remained remarkably high within artificial LCZs that cause SUHI intensification as the city grows. Low plants LCZ areas were warmer than Open low rise areas in all the periods. This could be due to the fact that the Open low rise areas constitute residents of the middle to high income strata who have resources to ensure that vegetation around their homes is well-maintained and healthy throughout the year. This is because vegetation in urban areas serves as temperature refugees in streets and parks providing cooling effect through evapotranspiration and shading [14]. According to Lu et al. [59], vegetation within buildings reduces patch sizes of built-up LCZs thus lowering their thermal effect on the surrounding LCZs. High SUHI in low plants LCZ contradicted with other studies such as Shi et al. [102] and Lu et al. [59] which had low plants as heat sinks. The disparity was because open low rise includes natural grassland areas, which in Bulawayo experienced drying of vegetation during the hot dry season. Low plants also include bare areas and croplands whose cover during the dry season could be bare or dry vegetation, reducing the surface cooling effect of latent heat transfer. Therefore, spatial and temporal variations in the thermal characteristics of low plants significantly reduced their heat mitigation value in Bulawayo. Higher UHI in 2005 and 2020 may partly be explained by higher average precipitation around satellite overpass dates in those years than in 1990. This is also in tandem with Amorim [71], that the heat mitigation value of urban greenery and SUHI varies with seasons and are increased during wet periods around overpass when the vegetation biomass is increased.

Conclusions
In order to understand the effect of urban growth on the thermal environment, the study used multi-temporal Landsat data to map LCZ with very high accuracy and retrieve SUHI for Bulawayo metropolitan city in Zimbabwe for 1990, 2005 and 2020. LST of the water LCZ were used as reference for quantifying SUHI intensities instead of the subjective traditional "rural-urban" LST difference. The high LCZ mapping accuracies were attributed to precise generation of training data and the robustness of the RF classifier, an ensemble based technique, compared to single classifier-based approaches. Furthermore, the use of dry and wet biomass periods significantly improved LCZ mapping accuracy in all years. In both intervals, i.e., 1990 to 2005 and 2005 to 2020, built LCZs monotonically expanded at the cost of vegetation-and water-based LCZs. Intensity analysis showed that the growth of lightweight low rise mainly targeted low plant areas. Deforestation in the city was expressed by the gain of low plants, which targeted dense forests. Intensity analysis also showed that the growth of compact low rise occupied mostly by low income strata avoided eastern direction where there are compact low rise and open low rise generally characterized by high cost per land unit. Due to expansion of built and polluting LCZs, the SUHI intensities rose monotonically during the study period. SUHI intensities varied between LCZs as they intensified with built-up proportion and density of tall building while decreasing with abundance of healthy vegetation. Based on the findings, the study concluded that human activities and growth induced LCZ changes have continued to trigger warming in Bulawayo. SUHI retrieval based on LCZ scheme proved effective in determining effects of urban growth on the thermal environment.

Data Availability Statement:
Remotely sensed data used in this study can be freely downloaded from United States Geological Survey (USGS) Earth Explorer website (www.earthexplorer.usgs.gov (accessed on 27 January 2022)). The training data used to map LCZ have been uploaded on the WUDAPT website (https://lcz-generator.rub.de/factsheets/3eb90c0ab4210886bdd0b23f1fc7dccae8 93c005/3eb90c0ab4210886bdd0b23f1fc7dccae893c005_factsheet.html) (accessed on 27 January 2022).