Modelling of the heat and the cold risks in Sofia and Varna – preliminary results

According to future climate projections, the expected thermal stress environmental conditions will get worse if the authorities do not apply appropriate measures for mitigation and adaptation to these changes. The health issues concerning air pollution and extreme temperatures have assumed great importance in recent years. The objective of this study is to estimate the thermal comfort in two of the biggest cities in Bulgaria - Sofia and Varna and their surroundings for the year 2017 by biometeorological indexes. We use computer simulations of the atmospheric parameters that define the thermal comfort indexes, by Advanced Research Weather Forecast Model WRF ARW version 3.9. We performed the simulations on four domains for 2017 with an output frequency of 1 hour. The outermost domain has a horizontal resolution of 9 km and encompasses the Balkan Peninsula. It uses initial and boundary conditions from the 0.25-degree NCEP Final Operational Model Global Tropospheric Analyses datasets with a time-frequency of 6 hours. The estimation of the thermal comfort conditions is performed with characteristics called indexes. The differences in the number of cases between the indexes are due to the specific defi - nitions and the meteorological factors that each of them takes into account. Some of these characteristics have applications depending on the specific thermal conditions.


Introduction
Thermal comfort and air pollution are among the most important factors for the human quality of life.According to future climate projections, the expected thermal environmental conditions will worsen, and the authorities should apply appropriate measures for mitigation and adaptation to climate change (Arias et al. 2021).The health issues concerning air pollution and extreme temperatures have assumed great importance in recent years (Monteiro et al. 2022), particularly in Bulgaria (Gadzhev and Ganev 2021;Georgieva 2021;Gadzhev et al. 2014).The studies of human thermal comfort in Bulgaria (Ivanov and Evtimov 2014a;Ivanov and Evtimov 2014b) and the Balkan Peninsula (Ivanov et al. 2020) suggest that the region is subjected to heat-related and cold-related illnesses.Other results based on different models for cold and hot thermal stress GeoStudies 1: 43-57 (2024), DOI: 10.3897/geostudies.1.e113477 Vladimir Ivanov et al.: Modelling of the heat and the cold risks at Sofia and Varna in big cities and their vicinities support similar findings (Ivanov and Dimitrova 2021;Gadzhev and Ganev 2021).They reveal that the thermal discomfort in Bulgaria could be a significant issue from the human point of view.
These studies use two of the so-called indexes for estimating the thermal comfort -Heat Index (Steadman 1979;Rothfusz 1990) and Wind Chill Index (Osczevski and Bluestein 2005).They have been used for warnings in some countries for many years (Environment and Climate Change Canada 2023; National Weather Service 2023).There are other widely used biometeorological indexes (Nastos and Matzarakis 2011;Di Napoli et al. 2018;Potchter et al. 2018), as the Predicted Mean Vote (Fanger 1972), the Physiologically Equivalent Temperature (Höppe 1999) and the Universal Thermal Climate Index (Błażejczyk et al. 2013), based on more comprehensive heat balance models of the human body in comparison to the previous ones.These indexes estimate both the heat and cold stresses within broad limits.Their differences are in the model definitions, parameters, and assumptions.We are going to use them together with the former ones to estimate the thermal stress conditions for human beings in order to get a wider perspective of thermal discomfort conditions.
The objective of this paper is to estimate the thermal comfort in the two of the biggest cities in Bulgaria -Sofia and Varna, and their surroundings for the year 2017, using several biometeorological indices.

Methods
We use computer simulations of the atmospheric parameters that define the thermal comfort by Advanced Research Weather Forecast Model WRF ARW version 3.9 (Skamarock et al. 2008).The simulations were performed on the HPC System "Avitohol" at the Institute of Information and Communication Technologies at the Bulgarian Academy of Sciences (Atanassov et al. 2017;Karaivanova et al. 2022) and at the NESTUM cluster.They were performed on four nested domains (Fig. 1) for 2017 with an output frequency of 1 hour.The outermost domain (D1), with horizontal resolution 9 km, encompasses the Balkan Peninsula and the most western parts of Anatolia Peninsula.It uses initial and boundary conditions from the 0.25-degree NCEP Final Operational Model Global Tropospheric Analyses datasets (NCEP 2000) with a frequency of 6 hours.The domain D2 with horizontal resolution 3 km is nested in the D1 and encompasses the territory of Bulgaria and some parts of neighboring countries.The third (D3) and the fourth (D5) ones with spatial resolution 1 km are nested in D2 and cover the Sofia Valley and the Varna municipality with some adjacent land and marine areas, respectively.The model configuration uses 50 pressure-based terrain-following vertical levels from the surface to 50 hPa.The data assimilation is switched on for the D1 for all vertical levels and for the first ten levels above the ground for D2.The model physics subgrid processes cannot be resolved explicitly because of the high resolution.
For that reason, we use parameterization schemes for simulation of the subgrid processes -surface physics with the Noah land surface model (Tewari 2004), surface layer processes with MM5 surface layer scheme (Jimenez et al. 2012), planetary boundary layer with Shin-Hong scheme (Shin and Hong 2015), NSSL 2-moment scheme (Mansell 2010) for the microphysical pro-cesses, Kain-Fritsch scheme for the cumulus convection (Kain 2004) and the new version of Radiative Transfer Model for shortwave and longwave radiation processes (Iacono et al. 2008).The cumulus parameterization scheme is switched off for D3 and D5 domains, because the convective processes at their scales are simulated directly.The reliability of the model results is supported by its validation in a previous study (Ivanov and Dimitrova 2021).The choice of the year 2017 for the study is dictated by the data availability, and because our intention is to show thermal comfort conditions for some typical subperiod.It turns out that 2017 fulfills the last condition, according to the course of the yearly maximum and yearly minimum temperatures for the simulated period (Fig. 2).The thermal comfort conditions are studied by a numerical characteristic of deviations from the thermal comfort called index.The current research deals with the following indexes.The Heat Index (HI) is used in high air temperature and humidity with four categories, shown in the Table 1 (Environment and Climate Change Canada 2023; National Weather Service 2023; National Oceanic and Atmospheric Administration 2023).It is the apparent temperature in the reference environment, corresponding to the feeling of heat in the current air temperature and humidity (Steadman 1979).The human body model assumptions for that index are constant walking and wind speeds with some clothing equivalent to long trousers and short-sleeved shirt and no solar radiation.The index is based on heat balance model of human body consisting of core and skin nodes, giving an account of convection, radiation, evaporation and heat loss from exhaled air (Steadman 1979;Rothfusz 1990; The NCAR Command Language 2019).The Wind Chill Index (WCI) (Table 2) considering the processes of convection and radiation from the human face (Tikuisis and Osczevski 2003).It is the feeling temperature in calm winds, corresponding to the heat loss from unclothed body parts (cheek), as in the current air temperature, wind speed, and no solar radiation (Osczevski and Bluestein 2005), and is usable in low air temperatures (Table 2).
The Predicted Mean Vote (PMV) is based on the thermal sensation vote from a large group of people (Fanger 1972).It is the index of the human comfort feeling (Table 3), accounting for the air temperature, wind speed, relative humidity, solar radiation, clothing and physical activity.The Physiologically Equivalent Temperature (PET) is the indoor air temperature, which is a result of the energy consumption of the human body equal to the one in the current outdoor conditions, keeping the clothing and physical activity constant (Höppe 1999).Its thermal stress categories follow the thermal sensation vote for the PMV are shown on the Table 4 (Matzarakis and Mayer 1996).The PMV and PET represent the human body model as two nodes -body core and body surface.The Universal Thermal Climate Index (UTCI) is the apparent temperature of feeling (Table 4) in reference conditions, which correspond to the same thermo-physiological reactions of the clothed human body in the current conditions with constant moving speed, physical activity and clothing (Błażejczyk et al. 2013).The heat balance model for that index uses adapting clothing model and divides the human body in 12 parts and multiple tissues (Fiala et al. 2012).
Each one of the last three indexes uses a human heat balance model with different and generally more comprehensive approach for calculating of their components and the body and clothing parameters, in comparison with HI and WCI.All indexes, except the PMV, are in the temperature dimension.The naming convention of their values categorizes the thermal discomfort by means of the feeling of the environment temperature.These indexes are based on different assumptions and models for heat balance of the human body, which implies a difference in their values.The calculation of the PMV, PET and UTCI is performed by the RayMan -software for environmental radiation and biometeorological modelling (Matzarakis et al. 2007;Matzarakis et al. 2010).The input variables for the calculations are the air temperature at 2 meters, the relative humidity at 2 meters, the wind speed at 10 meters, the surface temperature, cloud cover, sur- face albedo, along with time and location data.The clothing insulation is variable only in the calculation of the PMV with three different values for the winter, the summer, and the transition seasons -spring and autumn, respectively.

Results and discussion
The simulations for the HI in Sofia (Fig. 3A) show that there are cases for all categories.It determines the sultriness conditions mainly by the Caution and Extreme Caution, and the last two index categories are presented by one case in the rural areas.The most cases with the Caution and Extreme Caution are simulated in urban ones.That pattern for the Caution is typical also for the most populated areas outside Sofia.However, there are more simulated cases for the Extreme Caution also in the rural and sparsely populated areas in the domain.These spatial features suggest an influence of the effects of heat island and closed terrain form, leading to higher temperatures in the urban areas.
The most frequent heat stress categories for the Varna domain (Fig. 3B) are Caution и Extreme Caution.The numbers of cases decrease, approaching the Black Sea coast, more prominent for the Extreme Caution.The decreasing pattern of the last category holds even in the urban areas, in contrast to the Caution conditions.There are Danger and Extreme Danger cases in the non-urban south areas of the Varna domain, despite its proximity to the coast.The heat island effect could be the main factor for the urban areas.In contrast to the Sofia domain, the Black Sea proximity with the usual breeze circulation makes the heat discomfort milder, especially near to the coast.
The cold stress conditions simulations for the Sofia domain show results only for the first four categories of the WCI (Fig. 4A).The cold risk is determined mainly from the Low Risk and Moderate Risk categories.High Risk is simulated outside the urban area of the Sofia, specifically in the southern mountain areas.Very High Risk is presented in the Vitosha Mountain.The counts of High Risk and Very High Risk cases are almost negligible in comparison to the Low Risk and Moderate Risk.Taking into account the definition of the WCI, we suggest that the spatial features of the categories are due to the heat island effect (Macintyre et al. 2021) during the winter and the topography.The first factor could be responsible for the decreasing of the Moderate and High Risk conditions in the Sofia city, compared to the other areas.The second and third ones one are the closed terrain form, which suggest a lower sunshine duration and varying wind speeds due to changing between urban and rural terrains as well as the building structures.
The simulations for the Varna region (Fig. 4B) show results only for the first three WCI categories.There are no simulated cases with Very High Risk conditions  in contrast to the Sofia domain.The cold risk conditions are presented mainly by the Low Risk и Moderate Risk.The difference between these two categories is the much lower number of Moderate Risk cases, and their increase above 300 in the northern areas, probably due to the combination of the effects of the heat island and the buildings on the inner city wind.There are High Risk cases in the Varna domain only in the northern areas, but a small amount compared to the Low Risk and Moderate Risk categories.The simulated WCI patterns in the northern parts suggest a considerable influence of the north and northeasterly winds.
The simulations of the PMV in the Sofia domain (Fig. 5A) show a non-zero number of cases of each one of the index categories, with more homogeneous spatial distributions for the Cold and Cool categories.The urban areas are characterized by increased number of cases for the Very Cold, Slightly Cool, Neutral and Hot.The Vitosha Mountain is not presented with Warm, Hot, and Very Hot conditions.
The domain is characterized by many Very Cold cases and much smaller frequency of Warm, Hot and Very Hot.
The simulation results of PMV for the Varna domain (Fig. 5B) reveal that as for Sofia, the hottest conditions are much smaller in number than the other ones.Very Cold and Cold have smaller numbers of cases in the urban and the adjacent areas of the Varna domain.The Cool, Slightly Cool Neutral and Slightly Warm PMV show the opposite behavior.The differences in the PMV between the land and marine/coast areas are smallest for the Slightly Warm conditions, and increase from Warm to the Very Hot, and from Neutral to Very Cold PMV.The urban topography and heat island effect could have a prominent influence on these spatial features.We suggest that the typical for the coast breeze circulation could be а more important factor for the Cool, Slightly Cool Neutral and Slightly Warm, Hot and Very Hot during the summer and partly the transition seasons.It is evident in the sharper land-sea change in the cases numbers for the hotter categories.The bigger number of coldest conditions than the hottest ones suggests a prevalence of the cold weather or that the last one gives more weight to the PMV.
The PET simulation results for Sofia domain (Fig. 6A) show cases for all categories.The spatial distribution of the number of cases is relatively more  There are no substantial differences between the urban and other areas of the Varna domain (Fig. 6B), except for the No Thermal Stress.There are a smaller number of cases of Moderate Heat Stress and higher temperature discomfort for the coast and the northern areas.However, there are more cases of Extreme Cold Stress in the last ones.Probably, the main effect on the simulated results for the Varna domain is the breeze circulation for the hot discomfort conditions.Some influence could emerge from the north winds, evident by the higher numbers of Extreme Cold Stress in the northern parts of the domain and lower number of high temperature categories.
The simulated differences between the PMV and PET, despite their equal number of categories, are due to their specific determination and assumptions because the PET is applicable for warmer conditions.
The simulations for the UTCI in the Sofia domain (Fig. 7A) show that all categories are presented, but the Extreme Heat Stress is presented by just one case in several locations in the west.The spatial distribution of the number of cases is relatively more inhomogeneous for the coldest categories.There  The results for the Varna domain (Fig. 7B) show that almost all categories are presented, with isolated cases of Extreme Heat Stress in the south and inhomogeneous spatial distribution for the colder categories.The south urban parts of the Varna are notable for having more No Thermal Stress cases.There are a few Extreme Cold Stress cases as for the Sofia domain.It turns out for both Sofia and Varna domains that the cold categories of UTCI have generally higher frequencies than the warm ones.As we noted the same pattern in PMV, we assume that in 2007 there is cold weather prevalence concerning PMV and UTCI.The factors responsible for that distribution of cold and hot UTCI categories are the same as for the other indexes.However, the UTCI shows a larger spatial inhomogeneity of the thermal stress.In addition, the fact that there is no or negligible recurrence for the most extreme UTCI suggest that it covers a wider range of meteorological conditions than the other indexes.
According to the definitions of HI and WCI, the first one is applicable in warm/hot weather, and the second one in cold conditions.On the other hand, they use only air temperature and humidity or wind speed, which limits their use in more specific weather conditions.The PMV and UTCI are applicable for a broader range of meteorological conditions, which distinguishes them from the other three indexes.The PET is somewhat in the middle because, by definition, its categories are discernible generally for positive temperatures range.Therefore, the explanation for the bigger number of cases of Extreme Cold Stress PET in comparison to those of the Very Cold PMV is that the first one probably comprises part of the cases of the cold categories of PMV.Another feature is the difference between Strong Cold Stress PET and Moderate Cold Stress PET.There are more numbers of cases of Moderate Heat Stress and Strong Heat Stress PET than warm and hot PMV, but it is opposite for the hottest categories of the two indexes.These two patterns for cold and warm categories could give a warmer PMV-based stress or colder PET-based stress under the same conditions.The possibility for changing of the clothing and activity in the calculation of the PMV, but not in PET, would reinforce these differences.However, the more advanced heat balance model of the PET introduces more uncertainty in that conclusion.That suggests that the PMV and UTCI are applicable for broader meteorological conditions than the other indexes.The UTCI is the newest one and with the most sophisticated model of the human body heat balance.Probably that is the reason for its ability to categorize a wider range of meteorological conditions than PMV and PET, which is evident in two features.Firstly, the UTCI spatial distribution inhomogeneity is larger than that of the PMV and PET.Secondly, the UTCI has more locations with a zero number of cases for the extreme categories.Therefore, it is applicable in a wider range of outdoor thermal conditions than PMV and PET, and is probably better from the physical point of view.On the other hand, the HI and WCI have larger frequency than the UTCI in extreme conditions.However, probably UTCI outperform HI and WCI concerning the capability for modelling the human thermal comfort.

Conclusions
The HI in Sofia suggests risk mainly for fatigue and heat cramps, and to a lesser extent, from heat stroke.There are also several cases of increased risk from heat cramps, heat exhaustion and heat stroke.The WCI in Sofia is presented by its first three categories, which suggest that the low discomfort conditions are prevalent.The risks of hypothermia and frostbite are the second ones by recurrence and the risk of freezing is the third one.The hot and cold risks in the Varna domain decrease in comparison to that of Sofia.The spatial distribution of the PMV is relatively homogeneous in the Sofia domain, with the differences for the Vitosha Mountain because of the increasing altitudes.The coldest PMV cases numbers are much smaller in the Varna domain.The PET in Sofia domain is characterized by more cases of cold discomfort due to the index definition.The differences between urban and non-urban areas are notable for the cold, the neutral, and the first high-temperature discomfort categories.The Varna domain is characterized by fewer cases of PET heat discomfort.The UTCI has a more complex spatial distribution for Sofia and Varna domains.The hottest discomfort category is presented by several cases, meaning that the UTCI and the HI are appropriate for extreme heat risk warning.
The changing of the discomfort for the Sofia domain is due to the terrain form and proximity to the mountains, as well as the urban heat island effect.The Varna domain is influenced by the proximity of the Black Sea, with the breeze circulation in the warm season, the warming effect of the sea in the cold season and the geographic location and the topography in their vicinity.The differences of the simulated numbers of cases between the indexes are due to their specific definitions, number of categories, and the number of meteorological factors that each one of them takes into account.Another possible factor is the length of the simulation period.Further, we plan to use computer simulations for a longer period, which would reveal spatial and temporal features that cannot be observed yet.

Figure 3 .
Figure 3. Spatial distribution of the number of Heat index categories for A) Sofia domain and B) Varna domain.

Figure 4 .
Figure 4. Spatial distribution of the number of Wind Chill index categories for A) Sofia domain and B) Varna domain.

Figure 6 .
Figure 6.Spatial distribution of the number of Physiologically Equivalent Temperature categories for A) Sofia domain and B) Varna domain.

Figure 5 .
Figure 5. Spatial distribution of the number of Predicted Mean Vote index categories for A) Sofia domain and B) Varna domain.
is a higher number of No Thermal Stress cases in the urban part of the domain, which is opposite to the tendency for the Very Strong Cold Stress, Strong Cold Stress and Moderate Cold Stress.There are no or negligible number of cases of the Moderate Heat Stress and the hotter categories in the Vitosha Mountain.There are no simulated cases of Extreme Cold Stress in most of the urban areas.The

Figure 7 .
Figure 7. Spatial distribution of the number of Universal Thermal Comfort index categories for A) Sofia domain and B) Varna domain.

Table 1 .
Thermal stress categories of the Heat index.

Table 2 .
Thermal stress categories of the Wind chill index.

Table 3 .
Thermal sensation categories of the Predicted Mean Vote.

Table 4 .
Thermal stress categories of the Physiologically Equivalent Temperature and the Universal Thermal Climate Index.