Changes in urban heat island effect with the development of Newtowns

7 Newtown is a planned city built over a short time period. It is suitable for climate and thermal research, 8 particularly formulating urban planning strategies to analyse problems such as urban heat islands 9 (UHIs). Herein, a comprehensive approach was demonstrated for determining changes in UHI 10 distribution during 1989–2048 in two Newtowns with different urban planning. A significant increase in 11 built-up areas was observed from 1989 (< 5%) to 2018 (> 40%) in both Newtowns. However, this increase 12 significantly varied (approximately 12.25%) with urban planning in the areas where UHIs occurred 13 before and after development. Moreover, without effective mitigation, the built-up area in each Newtown 14 is estimated to increase to approximately 60%, and the surface UHI intensity in most areas to increase by 15 4 °C in 2048. Thus, these results combined with architectural assessment models can improve the 16 understanding of thermal environmental impacts of urbanisation and help mitigate heat island hazards. 17 Global population growth and urban expansion primarily cause land use and land cover (LULC) changes and 18 increases in built-up area. In 2018, approximately 55.3% of the world’s population resided in cities, among 19 which 60% will reside in cities with approximately 0.5 million inhabitants by 2030 1 . Rapidly increasing 20 economic development accelerates these changes, particularly in fast-growing urban areas, hindering sustainable 21 development 2 . LULC changes induced by human activities lead to different local climates than in surrounding 22 areas. This effect, termed as urban heat island (UHI), occurs worldwide 3,4 . UHIs primarily occur due to 23 increased solar radiation absorption and trapping in new surface materials of various infrastructure 5,6 . The 24 magnitude and

temperature 9 . SUHI is particularly evident in spatial variations of upwelling thermal radiance caused by LULC 28 changes and is commonly influenced by the surrounding sub-urban environment 8,9 . 29 Newtown, also called a planned city, is built in a short time period within a pre-determined boundary for 30 specific purposes. Since the mid-to-late twentieth century, Newtowns have been constructed worldwide, 31 contributing to population growth and inflation in large cities 10,11 . Newtowns facilitate climate and thermal 32 research through formulation of urban planning strategies to analyse problems, such as UHIs, and by providing 33 information on the urban temporal temperature variation mechanism 12 . Comparison of UHI changes in 34 Newtowns have not yet been conducted. Carrying out comparative studies on climate effects of urbanisation 35 under different urban planning conditions is particularly difficult because of different urban environments, 36 economic situations, and climates, as well as inconsistent data. 37 Since 1990, 14 Newtowns have been repopulated or built in sub-urban areas in South Korea to manage 38 population, transportation, and environmental concerns in several large cities. Urban planning in the first-39 generation Newtowns, providing indiscriminate housing, was not systematic and resulted in negative impacts, 40 such as unplanned urban expansion, environmental degradation, and low greenspace ratio in housing complexes. 41 The second-generation Newtowns were developed through systematic and environmentally friendly urban 42 planning, such as low-density urbanisation and expansion of green areas (Table 1). However, in both cases, an 43 increase in UHI is estimated because of a rapid infrastructural development and vegetation loss. Moreover, the 44 UHI phenomenon may intensify with further urban expansion. 45 Herein, expansion and intensification of UHI due to Newtown development was empirically analysed 46 using satellite data in two different-generation Newtowns in South Korea (Fig. 1). The SUHI intensity of each 47 newly constructed housing complex was found to be higher than in Pangyo Newtown, which led to increased 115 UHI intensity. 116 Predicted LULC for 2028, 2038, and 2048.
The cellular automata (CA)-Markov chain model (MCM) 117 analysis predicted that the proportion of built-up areas would increase by approximately 10% from 16.44 km 2 118 (49.16%) to 19.78 km 2 (59.12%) between 2018 and 2048 in Bundang Newtown (Fig 2a). Moreover, it predicted 119 decreases in forest areas from 35.61% to 29.9% and the grass cover from 12.76% to 10.69%. As Newtown 120 development in the past primarily occurred through transformation of agricultural areas to built-up areas, it was 121 not predicted that a significant urban expansion would occur through deforestation. In addition, most of the 122 buildings in the housing complex of Bundang Newtown were completed in 1990, over 25 years ago. Therefore, 123 renovations are planned for most of these old apartment complexes to improve the poor residential environment 124 and meet the latest urban housing requirements. Hence, most urban expansion was predicted to occur through 125 renovation within the existing built-up areas and partial transformation of the forest surrounding the Newtown. 126 In the case of Pangyo Newtown, the proportion of urban expansion between 2018 and 2048 was 127 predicted to be higher than that of Bundang Newtown. According to the CA-MCM prediction, built-up areas 128 would increase by approximately 18.42% from 40.81% to 59.23%, the forest areas would decrease from 40.84% 129 to 32.25%, and the grass cover including golf courses would decrease from 15.34% to 7.92% (Fig. 3a). The 130 primary trend observed in the predicted urban expansion was that non-urban areas, such as forest and grass, 131 surrounding the main road were transformed into built-up areas. In contrast with Bundang Newtown, Pangyo 132 Newtown is public-transportation-oriented. During the past Newtown development, the areas surrounding the 133 main road that existed outside the city were underdeveloped. However, if urban expansion occurs in the future, 134 it would be evident primarily in areas with good road proximity. In addition, urban expansion due to the 135 completion of development in the open spaces that were under development in 2018, and further development 136 within the city was also predicted. In terms of agricultural area and water, both Newtowns were predicted to 137 remain almost unchanged from 2018, with little fluctuation. 138 Predicted SUHI distribution for 2028, 2038, and 2048. CA-MCM predicted the increase in area and 139 intensity of the SUHI phenomenon in both Newtown and, unlike LULC prediction, a significant change was 140 predicted. In Bundang Newtown, the areas where the SUHI phenomenon occurs would increase by 141 approximately 5% between 2018 and 2048. For SUHI intensity distribution, the areas with SUHI ≤ 4 ℃ would 142 decrease from 17.12 km 2 (51.16%) to 11.44 km 2 (34.21%). Simultaneously, the areas with SUHI > 4 ℃ was estimated to increase from 4.25 km 2 (12.73%) to 10.68 km 2 (34.71%), affecting the lower SUHI intensity areas. 144 It is predicted that SUHI intensity would expand and increase from the existing residential area, which may 145 reflect the renovation trend partially occurring between 2000 and 2018. Therefore, development of sustainable 146 renovation guidelines is required such as thermal insulation, replacement of the insulation material, and 147 improving the air tightness of the building envelope through renovation using insulation materials 19 . In addition, 148 the areas with SUHI > 6 °C are predicted to increase from 0. 56 km 2 (1.7%) to 2.77 km 2 (8.28%). It has been 149 observed that the higher the LST, the higher the frequency of heat waves at regional scales 20 . In the future, 150 additional thermal environmental policies and energy policies are required for areas where SUHI intensity is 151 expected to increase significantly (Fig. 3a). 152 In the case of Pangyo Newtown, the areas where the SUHI phenomenon occurred were predicted to 153 increase by 20%. The affected areas are similar to those predicted to change from forests existing around the 154 main road to built-up areas. For SUHI intensity distribution, the area with SUHI ≤ 4 ℃ would decrease from 155 7.75 km 2 (43.97%) to 5.08 km 2 (28.83%). Moreover, the areas with SUHI > 4 ℃ would increase from 2.53 km 2 156 (14.34%) to 8.7 km 2 (49.36%), and most areas were in the range 2 ℃-4 ℃ (49%) (Fig. 3c). Therefore, it can be 157 predicted that urban features, such as structural characteristics, materials, and building disposition type would 158 change according to the housing complex newly built through Newtown development. 159

Discussion 160
This study is the first attempt to simulate and compare the pattern of UHI occurrence according to Newtown 161 development using remote sensing and GIS technology. This discussion focuses on the principal two 162 contributions of the proposed research in comparison with previous studies. Afterwards, the limitations are 163 discussed. 164 The main contribution of our study is that the different patterns of changes in land use land cover and 165 SUHI phenomenon depending on urban planning were visually and quantitatively shown for the study sites While the presented study provides useful method and information regarding the current and future status 202 of the UHI phenomenon, it is still faced some limitations. This study does not consider additional parameters 203 typically influencing the urban growth because of the specificity of the study area. As mentioned, Newtown is 204 the planned city where the physical and legal aspects of the site were reviewed through feasibility analysis 205 beforehand, the complication associated with urban expansion is relatively low for Newtown. However, the 206 factors for urbanisation are related to the complexity of the terrain, degree of socio-economic development, 207 urban regulations, etc 24 . Therefore, it is necessary to consider additional factors for urban expansion when 208 applying this methodology to a region other than Newtown in the future. In addition, a model that explains the 209 detailed behaviour of UHI using a combination of building renovation and structural characteristics is still 210 necessary. Future research studies should attempt to obtain structural and temporal data over the same period of 211 time and develop models able to explain the change of UHI based on structural characteristics changed by 212 building renovation. 213

Conclusions 214
Although the research methods and measures face certain conceptual and practical challenges, this study 215 suggested a proximate causal relationship between urban expansion and SUHI phenomenon change according to 216 urban planning. It is easy to apply for practitioners and the necessary data for application are available without 217 complex acquisition procedures or unopened access datasets. Therefore, the proposed novel method may be 218 applied to both existing and newly-built cities to predict future UHI distribution according to urban planning. 219 Furthermore, the findings and methods constructed through this research can be useful to policy makers, urban 220 planners, researchers, and citizens to adopt sustainable thermal environment management practices including 221 adaptation and mitigation strategies for the city. atmosphere and the surface radiative properties that influence the emission and reflection of radiation within the 231 spectral wavelengths detected by the sensor 9 . Atmospheric correction using the dark object subtraction (DOS) 232 method and radiometric correction for pre-processing using the semi-automatic classification (SCP) plugin in 233 QGIS 3.14, were applied to the images. Atmospheric scattering and absorption caused the imaging system to 234 Assessment of classification accuracy is necessary to ensure that classification data can detect changes; 252 this was conducted on the resulting classified imagery through an error matrix and kappa index that enables 253 differentiation between ground-truth and predicted classification 24,27 . High-resolution Google Earth data and 254 aerial photographs provided by the National Geographic Information Institute (NGII) of South Korea were used 255 to establish ground-truth regions for the evaluation of classification accuracy (http://map.ngii.go.kr/). High-256 resolution data from Google Earth have been used as reference in many classification studies and national 257 standardised land cover maps; NGII provides high-resolution aerial photographs captured since 1945, and can 258 also be used for accuracy assessment 22,24,28 . The kappa coefficient was calculated using equation (1): 259 (1) where i is the class number; n is the total number of points; nii is the number of pixels belonging to the actual 260 data class i, which were classified as class i; Ci is the total number of classified pixels belonging to class i; and 261 Gi is the total number of actual data belonging to class i. Fifty sample points per class for each Newtown, except 262 water class, were selected automatically by QGIS 3.14. A minimum of 50 samples must be collected for each 263 land cover class in the error matrix to avoid the risk of a biased sample during accuracy assessment 29 . 264 LST estimation. LST estimation using ArcMap 10.5 includes transforming DNs to radiance (Lλ), measuring 265 radiance brightness temperatures (TB), and adjusting emissivity to extract surface temperature from brightness 266 maps 30 . The LST values were obtained using thermal bands from Landsat TM (B6) and Landsat OLI/TIRS 267 (B10) because of the USGS recommendation to avoid using TIRS band 11 because of its higher calibration 268

uncertainty. 269
Every object on the Earth emits thermal electromagnetic radiation when its temperature is above absolute 270 zero (K), and the signal received by the thermal sensors can be transformed to radiance (Lλ) using equation (2): 271 where Lλ is the spectral radiance in W/(m 2× sr×μm); ML is the radiance multiplicative scaling factor for the band; 272 AL is the radiance additive scaling factor for the band; and Qcal is the level 1 pixel value in DN, whose values are 273 obtained from the metadata of the Landsat images. After the DN value was converted to radiance, the radiance 274 values were converted to TB using equation (3): 275 where TB is the At-satellite brightness temperature and K1 and K2 represent the band-specific thermal conversion 276 constants from the metadata. To obtain the temperature in Celsius, the radiant temperature is revised 30 . The final 277 step in estimating the LST is to rectify the TB using land surface emissivity (LSE, ε) correction as shown in The obtained values of TB were referenced as a black body, whose properties are different from that of 283 real objects on the Earth's surface and would also be different from real LST 33 . The LST values across a city can 284 have a wide range, and it depends on LULC states constructed within the city. Furthermore, LSE, which is 285 essential for estimating the LST, has strong land use/land cover dependence 34,35 . 286 The LSE value is calculated conditionally using equation (5)  hybrid and robust algorithm in spatial and temporal dynamic modelling of LULC changes that includes the 303 deterministic modelling framework, spatially explicit approach with stochastically based temporal 304 framework 43,44 . In addition, CA-MCM analysis allows the user to add factors related to urban expansion into the 305 model to improve accuracy, and it can be a support tool for land use planners and policy makers to establish 306 future land use policies 45 . Furthermore, MCM is a tool used to evaluate adjustments in land use among cycles by 307 a sequence of values that depend on the present state 46 . MCM defines the present temporal LULC change to 308 predict future change, and equation (8) presents the calculation of land use change prediction 47 : 309 where S(t) is the system state at time t, S (t+1) is the system state at time t+1, and Pij is the transition probability 311 matrix in a state, which is calculated using equation (9). A contiguity filter of 5×5 pixels was used to define the effect of neighbouring pixels on the central pixel. 347 Mapping and prediction of the SUHI distribution. The UHI effect occurs due to the anthropogenic 348 modification of natural landscapes in the city boundary layer, and as the urban area increases, the UHI intensity 349 also increases 14 . In addition, LST and SUHI effects are particularly related to the surrounding sub-urban 350 environment 8,14 . To analyse this trend, the SUHI intensity of each Newtown was defined as the difference 351 between the temperatures of an urban area and its surrounding areas (LULC, excluding built-up area) within the 352 boundary 4,13,15 . Thus, the SUHI intensity distribution maps for each Newtown and each period were constructed 353 using two steps. (1) The SUHI intensity variation was calculated using equation (10): SUHI intensity distribution = Ts -(T mean + 0.5 × δ) surrounding area (10) where Ts is the LST (℃) distribution of Newtown, and T mean and δ are the mean and standard deviation of LST 355 in non-urban areas of Newtown. By subtracting the average temperature of non-urban areas from the 356 temperature of the entire city, it may be verified that the actual SUHI effect was due to urban expansion, rather 357 than the temporary LST value. In addition, the water bodies were excluded while calculating the SUHI intensity 358 because it can irregularly influence the surface temperature (Lee et al, 2020). (2) The SUHI intensity variation 359 was classified into six appropriate ranges: (ⅰ) value ≤ 0 ℃, (ⅱ) 0 ℃ < value ≤ 2 ℃, (ⅲ) 2 ℃ < value ≤ 4 ℃, (ⅳ) 360 4 ℃ < value ≤ 6 ℃, (ⅴ) 6 ℃ < value ≤ 8 ℃, (ⅵ) 8 ℃ < value. Thus, the difference in distribution and intensity 361 of the SUHI phenomenon can be compared according to the change in LULC for each Newtown at each time 362 period. In addition, classes are divided into value ranges, to facilitate future SUHI intensity distribution 363 prediction using CA-Markov analysis. The indices, which were positively and negatively correlated with LST, 364 were used to develop transition suitability maps for predicting the SUHI distribution. The normalised difference 365 built-up index (NDBI) was used as the index that highly correlated with LST 56 . NDBI is the most widely 366 accepted tool for the identification of built-up areas and has shown a high surface temperature correlation in 367 previous studies 13,22,57 . The NDBI value was calculated using equation (11): 368 (11) Built-up areas are sensitive under the 1.55-1.75 wavelength range in the short-wave infrared (SWIR) band; 369 however, they are less sensitive under the 0.79-0.90 wavelength range in the NIR band 58 . The NDBI values 370 range from -1 to +1, and values near +1 generally represent highly dense built-up areas. Furthermore, NDVI was 371 used as the index that weakly correlated with LST. NDVI is the most common index for vegetation extraction 372 and has shown a strong negative correlation with LST in previous studies 32,57,59 . Fuzzy membership functions 373 were also used to standardise the factor maps to 0-1, where 0 represents a low SUHI potential and 1 represents a 374 high SUHI potential. 375

Data availability 376
Satellite images from 1989 to 2018 used in this study are freely available at httl://earthexplorer.usgs.gov/. Other 377 datasets are available upon request from K. Lee (leedake@korea.ac.kr). 378