How urban growth changes the heat island effect and human thermal sensations over the last 100 years and towards the future in a European city?

Urban development exacerbates urban heat island effects of cities and thermal discomfort of the residents worldwide. However, such urban effects vary with the geographical and climatic conditions of the respective cities. This study is the first to quantitatively estimate the impact of past (1878), present (2012), and future (2050) urbanization on surface air temperatures and human thermal comfort in Sofia, the largest city of Bulgaria. We used the state‐of‐the‐art Weather Research and Forecasting (WRF) climate model, with a 1‐km horizontal resolution for three cases: LU1878, LU2012, and LU2050. We first validated the results of the control simulation (case LU2012) against actual observations. The model satisfactorily reproduced the diurnal variation of the observed monthly mean surface air temperatures for July (2011–2013) with a mean bias of ±0.6°C. We then estimated the potential impacts of past and future urbanization on surface temperature. Our results suggest that because of urbanization of Sofia, the surface air temperature at 0600 Local Standard Time (LST) has increased approximately 4.0°C in the past (1878–2012) and approximately 0.4°C in future (2012–2050). We also evaluated the impact of the urban thermal environment on human comfort levels. Considering the LU2050 case in the central part of Sofia, the temperature–humidity index during mid‐afternoon (1400–1700 LST) is estimated to be 24–26, indicating that approximately 50% of the residents could feel perceptibly uncomfortable.


| INTRODUCTION
In recent times, many countries have faced challenges presented by their growing urban populations, including health and human comfort, energy and infrastructure, and increasing temperature and urban heat islands (UHIs).
Historical and future impacts of urbanization on UHIs have been investigated for many cities across the world e.g., Tokyo, Japan (Kusaka et al., 2000Adachi et al., 2014), Sendai, Japan (Vitanova et al., 2019), Beijing, China (Miao et al., 2009), and Phoenix, United States (Georgescu et al., 2009a(Georgescu et al., , 2009b. For the Tokyo metropolitan area, Kusaka et al. (2000) numerically estimated the UHI intensity on a clear summer day to be 3.0-4.0 C due to landuse changes over the past 85 years. Doan et al. (2019) estimated that in the urban core area of Hanoi, Vietnam, the monthly mean temperature for July increased by 0.35 C as a result of urbanization during 1990-2010.
This study focuses on the largest Bulgarian city, Sofia. In recent decades, very few studies have been conducted on the thermal environment and UHIs in Sofia (e.g., Dimitrova et al., 2019;Kolev et al., 2000;Pongr acz et al., 2009;Syrakova & Zaharieva, 1998;Vitanova & Kusaka, 2018). However, there have been no research studies vis-a-vis Sofia on how urbanization alters the UHI effect and human thermal sensitivity. In fact, we are not aware of any such study being undertaken for any Eastern European city.
The primary goals of this study were to: a) numerically examine the potential impact of past (1878), present (2012), and future (2050) urbanization on surface air temperature of Sofia, and b) evaluate the impact of the urban thermal environment on human comfort.
Results from this study will provide valuable information for further research on UHIs, human health, energy-saving, adaptation and mitigation strategies, future projections, and urban planning. Additionally, the results of this study could be potentially applicable to other cities that exhibit similar social and demographic characteristics.
2 | DATA AND METHODOLOGY

| Study area
Sofia became the capital of Bulgaria in 1878, with a population of approximately 11,000; however, with a history of approximately 7000 years, it ranks among the oldest cities in Europe. The city's population crossed 500,000 in 1946, and presently Sofia is the largest city in Bulgaria with a population of more than 1.2 million.
Located in western Bulgaria, Sofia is surrounded by the Vitosha and Balkan Mountains in the south and north, respectively, with a mean elevation of 550 m (Figure 1b). The climate is humid continental, as defined by Köppen-Geiger (Kottek et al., 2006). The hottest months are July and August, with monthly mean temperatures of 21.4 and 21.8 C, respectively, recorded during the period 2016-2020 (Stringmeteo, 2020).

| Model configurations and design of numerical experiments
The Weather Research and Forecasting (WRF) model version 3.5.1 (Skamarock et al., 2008) with a 1-km resolution was used to conduct the LU1878, LU2012, and LU2050 cases (Table 1). For the numerical experiments, we used four nested domains with spatial resolutions of 27, 9, 3, and 1 km, and 30 vertical layers that increased with the altitude. The time period considered for the analysis was from 1 July 0000 to 30 July 2300 Local Standard Time (LST) between 2011 and 2013, i.e. 3 years ( Table 2).

| Land-use/cover data
The land-use data for the LU1878 case were created based on the 3rd military mapping survey of Austria-Hungary, 1868-1880 (Department of Map Science and Geoinformatics, Eötvös Lor and University, n.d.) and CORINE Land Cover (Copernicus Land Monitoring Service, 2012). The LU1878 case consisted of one urban category ( Figure 2a). The average urban fraction was estimated using the city plans of Sofia from 1878 (not shown). Here, the cropland and forest areas are the same as those in the LU2012 case (Figure 2b), except for the artificial water surfaces.
The land-use data for the LU2012 case were created from the CORINE Land dataset 2012 (Copernicus Land Monitoring Service, 2012). Detailed urban classification of Sofia was obtained from Urban Atlas 2012. Information for the AH release of the LU2012 case is provided in the next section.
The land-use data of the LU2050 case were created based on the CORINE Land dataset 2012. We consider that land use in the 2050s is almost same as in 2012. The reasons for this assumption are as follows: (1) during 2012 and 2018, the urban area of the entire city has expanded very little (CORINE Land Cover 2012 and 2018), indicating that Sofia is already a well-established city and is unlikely to grow any further. (2) According to 'Sofiaplan, 2019' and 'Vision Sofia 2050', the existing houses and buildings are sufficient to accommodate the citizens already living in the Municipality area. Thus, it is widely considered that further urban expansion is not required in this stage for Sofia.
For the LU2012 and LU2050 cases, we considered three urban categories: low, medium, and high residential areas ( Figure 2b and Table 3), i.e. where the urban fraction of a grid cell was 0.05-0.55, 0.55-0.85, and >0.85, respectively. A change in urban fraction would indicate a change in spatial distribution of urban land categories.

| Anthropogenic heat fluxed data
The total anthropogenic heat (AH) flux release Q f (W m À2 ) was calculated using a top-down approach based on a simple calculation (Sailor & Lu, 2004): where Q b = Q be + Q bh . Q f was estimated from the heat produced by vehicles, Q v , the heat from buildings, Q b , and human metabolic heat, Q m . Here, Q b was calculated based on the heat produced by electricity, Q be , and heating fuels, Q bh . The average released Q m was assumed to be 100 W per person (Sailor & Lu, 2004). (Or details regarding the total AH release and the diurnal profile of the AH, see Vitanova & Kusaka, 2018).
We considered that the AH release in 1878 was 0-5 W m À2 , even though there was no electricity or automobiles. AH could be non-zero because of commercial activities during that period. The assumed AH of 5 W m À2 in 1878 is comparable with the rural areas of Bulgaria. According to the photographic archive provided by Stara Sofia, the buildings in 1878 were mostly small singlestorey houses interspersed with narrow, unpaved streets.
In order to clarify the uncertainty of AH in 1878, we conducted two simulations for the LU1878 cases with 0 and 5 W m À2 . The results of both simulations are almost same. Thus, we will show the results from the LU1878 case with AH of 0 W m À2 in this study (Table 3). Based on the model assumptions for decarbonization of electricity generation by 2050 approved by the European Commission, we considered that electricity consumption for Sofia would rise approximately 0.6% according to the PRIMES (EU reference scenario, 2016). Also, the population of Sofia is predicted to be approximately 1.4 million in 2050 based on a realistic consideration of the demographic and socio-economic development of the city projected by the National Statistical Institute (NSI), Bulgaria (2020). Based on the projected population and assuming a minor increase in the per capita energy consumption, we estimated for the LU2050 case the AH would rise by 15% over the LU2012 level (Table 3). Table 3 shows the building and ground heat capacity and thermal conductivity, which are important parameters for the calculation of urban geometry in the WRF model for each urban category in this study. 2.5 | Temperature-humidity index The temperature-humidity index (THI), derived from temperature and humidity, is a measure of the degree of discomfort experienced by an individual during warm weather. In this study, we used the THI to assess the impact of urban thermal environment on human health (Equation 2) where t is the air temperature ( C) and RH is the relative humidity (%) at 2 m above ground level. THI values 21-24, 24-26, and >26, respectively, indicate that 0%, 50%, and 100% of humans feel uncomfortably hot (Niewolt, 1977).   Shortwave radiation Dudhia Shortwave scheme (Dudhia, 1989) Longwave radiation RRTM scheme (Mlawer et al., 1997) Cumulus Kain-Fritsch scheme (for d01 only) (Kain, 2004) 3 | RESULTS

| Land-use changes
We compared the land-use maps of Sofia for 1878 and 2012. Compared with 1878 (Figure 2a), the residential area in 2012 ( Figure 2b) shows a significant increase across the entire city due to increased population and urbanization. Compared with 2012, by 2050, there is expected to be a minor change in land use of Sofia due to its geographical location and its proximity to mountains (not shown). The concept of 'Vision Sofia 2050' for a compact and concentrated city by 2050 was also considered.

| Model validation of surface air temperatures and relative humidity
This study used observation data provided by Stringmeteo, the National Institute of Meteorology and Hydrology, and TV Met. The results of the control simulation (case LU2012) were validated against actual observations. While the simulation data used for the LU2012 case are the same as those of the control case reported by Vitanova and Kusaka (2018), the validation approach differed slightly. The results show that the WRF model satisfactorily reproduced the hourly average surface air temperatures and relative humidity (RH) at four observation stations over the month of July consistently for the period 2011-2013 (Figure 3; see Figures 1b and 2b for locations). The mean model biases, (i.e., the arithmetic differences) for the monthly mean air temperatures and RH ranged from À0.6 to 0.6 C and À 3.96% to À0.18%, respectively. For both the measured and the simulated data, we also calculated the uncertainty or error range, i.e., between quantiles 0.25 and 0.75 of data at each hour (Figure 3). For a certain hour, it was based on 90 values (i.e. 3 years Â 30 days). The differences between the simulated values and observations were small, and in most cases fell within the error range of the observed data.

| Potential impact of past, present, and future urbanization on surface air temperatures
Urban expansion and population growth are known to be factors that affect surface air temperature increase in urban areas. The surface air temperature in central Sofia at 0600 LST for the LU2012 (present) and LU2050 (future) cases was approximately 4.0 C higher than that for the LU1878 (past) case owing to land-use change and AH release during periods 1878-2012 and 1878-2050 (Figure 4a,b). The simulated results also showed that the difference between the temperatures in the surroundings of the central area at 0600 LST for the LU2012 and LU2050 cases was approximately 0.4 C because of the impact of the future AH increase (Figure 4c). The surface air temperature around the central area at 1500 LST for the LU2012 and LU2050 cases was, respectively, 0.8 and 1.0 C higher than that for the LU1878 case (Figure 4d,e). Comparing the simulated results at 1500 LST between the LU2012 and LU2050 cases, we find that the difference between the temperatures in the urban areas of Sofia was 0.1 C (Figure 4f). We computed the minimum (0600 LST) and maximum (1500 LST) surface air temperatures at the Central and Boyana stations for all study cases ( Figure 5). We selected the Central station as a typical built-up, busy, and high-density urban area. Conversely, the Boyana station represented a typical low-urbanized area with low urbanization surrounded by forests (Figure 2b).
It is noteworthy that the air temperature at 0600 LST in the LU1878 case for the central station was 1.2 C lower than that for Boyana (Figure 5a), although the sensible heat flux (transferred from the surface to the atmosphere) at the Central station is higher than that at Boyana (Table 4). We speculate that the advection effect could explain this difference. Because the Central station is located at the valley bottom, Boyana is on the valley slope (about 200-m difference from the valley bottom). Mountain winds, which occur during night, blow cold air downslope, filling the valley with cold air. This explains why the temperature at the Central station was cooled down deeper than that at Boyana in LU1878. However, as the UHI effect is much stronger in the LU2012 than in LU1878, the cooling effect of mountain wind becomes not obviously seen.
Moreover, we found high UHI intensity (air temperature anomaly between the Central and Boyana stations) of 2.7 C at 0600 LST in both the LU2012 and LU2050 cases, though no difference in UHI intensity between the two cases is recognized. The high UHI intensity is due to the increased urban fraction and AH release in the highdensity residential area where the ground heat flux releases more heat back into the air compared with that at the Boyana station (Table 4). The higher UHI intensity in the early morning (0600 LST) versus that during the afternoon (1500 LST) is consistent with the results reported by several other studies (e.g., Memon et al., 2009;Yadav et al., 2017). At 1500 LST, air temperature at the Central and Boyana stations increased by 0.3 and 0.6 C, respectively, for the LU1878 and LU2012 cases. The lower increasing rate at the Central station is attributed to a higher increase in thermal inertia as more vegetation land areas were converted to urban structures,  (WRF)) results from the control simulation (case LU2012); black circles are observational (OBS) data; the red and grey colours signify the error ranges for model and observation data, respectively compared with that at Boyana. The unnoticeable change in UHI intensity between two cases LU2012 and LU2050 ( Figure 5) is interpreted as future urban change will not play a significant role on temperature climate in Sofia.

| Impact of human thermal comfort (THI)
The impact of the urban thermal environment on human comfort was also evaluated. For the LU1878 case, when the city area was less developed (Figure 6a), THI varied between 21 and 24, indicating that 100% of the citizens felt comfortable, in both the Central and Boyana stations. However, for LU2012 and LU2050, we found that in the already-developed central Sofia, during mid-afternoon (1400-1700 LST, the period of the usual daily maximum temperature), THI varied between 24 and 26. This indicates that approximately 50% of the residents could feel uncomfortable. Moreover, THI for LU2012 and LU2050 cases was similar (Figure 6b,c), suggesting that if Sofia follows the concept of a compact city and adopts a F I G U R E 4 Differences in the simulated July mean surface air temperatures between the cases (a) LU2012 and LU1878 at 0600 LST, (b) LU2050 and LU1878 at 0600 LST, (c) LU2050 and LU2012 at 0600 LST, (d) LU2012 and LU1878 at 1500 LST, (e) LU2050 and LU1878 at 1500 LST, and (f ) LU2050 and LU2012 cases at 1500 LST. Large dotted circles represent Sofia and its surrounding areas F I G U R E 5 Simulated means surface air temperature for July (2011)(2012)(2013) at the Central and Boyana stations for the LU1878, LU2012, and LU2050 cases at (a) 0600 LST, and (b) 1500 LST moderate energy consumption policy, the THI may not increase substantially over the next 30 years.
Additionally, the above results suggest that the central area of Sofia should be improved in terms of thermal comfort along the lines of Boyana. This could be achieved by increasing and improving the city's green cover, integration of green infrastructure, particularly in areas having fewer green spaces, developing green corridors, and implementing tree-planting activities.

| DISCUSSION
This study is among the first to address urbanization-related warming for a major city in Eastern Europe, enriching the current understanding of the issue, which has implications for large cities around the world. Using numerical modelling approach (which has been well verified against in situ observational data), we show that the future urban changes in Sofia would be unlikely to contribute substantially to urban warming. The possible change in air temperature for Sofia due to urbanization is approximately 0.4 C, much lower than the rise of 4.0 C during the last century. Our results are in line with those of previous studies on cities such as Paris, Brussels, and Tokyo (Adachi et al., 2012;Hamdi et al., 2015;Kusaka et al., 2016), which suggest that future urban change may not contribute much to urban warming, but global climate change will do. Researchers argue that since cities like Paris, Brussels, and Tokyo are already well-developed, there is little scope for their further expansion in the future. Sofia is similar to those cities in this sense, as the future urban expansion is not a likely scenario T A B L E 4 Simulated sensible heat (SH), latent heat (LH), and ground heat (GH) fluxes at the Central and Boyana stations for the LU1878, LU2012, and LU2050 cases On the other hand, our results contradict some studies that argue the imperative impacts of future urbanization on surface air temperature (i.e., Doan et al., 2019;Georgescu et al., 2013). Notably, such results emerge from studies on cities where topography allows for urban extension possible, and which are characterized by dispersed urban structures. For example, Doan et al. (2019) found that AH in Hanoi did not contribute significantly to the temperature increases during 1990-2010. However, during 2010-2030, the increase in AH is estimated to contribute 30-50% of the total UHI, comparable with that due to land-use change.
This study is limited to discussing the primary signal of urban change on the temperature regime in Sofia. We did not account for the impact of background climate changes induced by greenhouse gas emissions. Our discussion focused on urban impact, without investigating the interaction between the urban area and other mesoscale effects such as mountain-valley circulations. It is noteworthy that Sofia being a mountain city at a relatively high elevation, the impact of the topography on the surface air temperature could account for the high UHI at 0600 LST in this area (Figure 4a, b). Understanding the urban-mountain interaction and how it could change in the context of global warming will be the goal of our future studies.

| CONCLUSIONS
This is the first study to quantitatively evaluate the potential impact of past (1878), present (2012), and future (2050) urbanization on surface air temperatures and human thermal comfort in the city of Sofia, Bulgaria's capital. The conclusions of this study are as follows: 1. Control simulation results (case LU2012) were validated against observations. The Weather Research and Forecasting model could satisfactorily reproduce the diurnal variation in the observed surface air temperatures and relative humidity. 2. The surface air temperature in the central area of Sofia at 0600 Local Standard Time (LST) for LU2012 case was approximately 4.0 C higher than that for the LU1878 case due to the impact of land-use change and anthropogenic heat (AH) release during 1878-2012. 3. The difference between the temperatures in the surroundings of the central area of Sofia at 0600 LST for LU2012 and LU2050 was approximately 0.4 C because of the impact of future AH increase.
4. The results confirmed that the contribution of urbanization to surface air temperature was larger in the past than that in the future because Sofia is already a well-developed city. 5. Vis-a-vis the impact of the thermal environment on human comfort, in central Sofia, in the mid-afternoon, temperature-humidity index for LU2050 was estimated to vary between 24 and 26, meaning that approximately 50% of the residents could feel uncomfortable.