Changes Detected in Five Bioclimatic Indices in Large Romanian Cities over the Period 1961–2016

: Bioclimatic indices are very important tools to evaluate the thermal stress of the human body. The aims of this study were to analyze the general bioclimatic conditions in ten big cities in Romania and to ﬁnd out if there has been any change in ﬁve bioclimatic indices over a 56-year period: 1961–2016. The indices considered were: equivalent temperature, e ﬀ ective temperature, cooling power, universal thermal climate index and temperature-humidity index. They were calculated based on the daily meteorological data of air temperature, relative humidity, and wind speed recorded in 10 weather stations in Romania: Bucharest-B ă neasa, Botos , ani, Cluj-Napoca, Constant , a, Craiova, Galat , i, Ias , i, Oradea, Sibiu and Timis , oara. The features investigated for trend detection consisted of the frequency and length of the occurrence period for each class and for each index. The test used for trend detection was Mann-Kendall and the magnitude of the trend (the slope) was calculated by employing Sen’s slope method. The main results are based on frequency analysis. Three indices showed comfort class as dominant whereas the other two indicated cold stress conditions as dominant in the area. There was a shift from the cold stress conditions to the warm and hot ones for all the indices. The most stressful conditions for hot extremes did not indicate signiﬁcant change. The climate in the big cities of Romania became milder during the cold season and hotter during the warm period of the year. The analysis of the length of each thermal class indicated mainly longer occurrence periods during the year for comfortable or warm stress classes.


Introduction
Much interest is shown nowadays, in the era of climate change, to information about the local impact of climate changes in every domain of our lives, from politics to public health and even science [1]. It is already widely known that one of the main consequences of global air temperature rises is the increasing frequency of intense, extreme events [2][3][4]. They lead, mostly during extreme seasons, to severe impact on different social and economic sectors, causing serious health problems to the population and disturbing agriculture, transportation, building industry and tourism activities. It is considered that the most affected climatic variables will be temperature, precipitation, humidity and wind, mostly because of their increase in extreme values, in terms of frequency, intensity and persistence [5]. Furthermore, international specialists in territorial planning state that climate change also affects, in a negative way, life quality in urban areas, especially through heat waves that occur in urban climate islands (ICU), to an increase in the energy and water consumption, as well as to an increase in polluting substances. Additionally, there is a higher risk of infrastructure damage and winter tourism. For all these, the mitigation and adaptation strategies propose both general and specific measures such as the continuous information and education of the population, preservation of the natural resources, encouraging infrastructure and geo-system resilience processes, as well as expanding green spaces and parks, building green infrastructure, innovation and diversification of winter tourism through finding new solutions independent of snow [6]. However, some of the effects of ICU under high temperature conditions (e.g., heat waves intensification) on human body are represented by the extreme thermal stress for the cardiovascular system and increase in the mortality rate, especially for children and elderly people [7,8]. Therefore, some of the proposed solutions by the same institution are developing pilot projects for acclimatization, infrastructure and green spaces, as well as the development of local, regional and national climatic change adaptation strategies.
In this context, common people show a high interest for knowledge, forecasts and predictions about weather conditions [5,9].
It has been believed since ancient times that weather plays an important role on the human body, but the latest scientific research has proved it. The sensitivity of the organs and the psycho-physiological reactions are increased by external atmospheric conditions. The adaptation ability to sudden changes in weather conditions differ from individual to individual and widely depends on genetic predisposition and specific characteristics [10]. Over the last 20 years, a multitude of studies have been conducted in order to explain the relationships between humans and their environment, trying to investigate how thermal comfort or thermal stress in outdoor environment influences human behavior during daily activities. Two comprehensive reviews have been developed [11,12], whereas some other research papers focused on specific issues. The great majority of these studies have focused on thermal sensations in urban areas [1,[13][14][15][16][17][18][19][20][21][22][23][24][25][26].
One of the most realistic and objective ways to assess thermal perception and stress for humans is based on using appropriate indices [1]. Even though inside the meteorological and biometeorological community there are different opinions about the indices that should be used or not, in general it is considered that the more indices used, the better and more reliable the image of the changes is created [5,27]. Indices, serving as tools for heat stress and thermal comfort analysis, have a major role in describing the combined effect of meteorological variables on humans in terms of thermal stress or comfort [28].
In Romania, the number of biometeorology studies is quite small, most of them being focused on bioclimatology and more likely presenting theoretical aspects [29][30][31][32][33][34][35]. In 2008, impressive work was done by Nicoleta Ionac and Sterie Ciulache that represented the Romanian Bioclimatic Atlas [36]. The number of studies that integrated biometeorological indices is quite low. One of them focused on the analysis of some biometeorological indices in the southern Dobrogea, which is one of the most hot and dry regions of the country [37]. Most recently, a bioclimatic analysis in the context of urban environment and tourism was developed, based on the temperature-humidity index, described by fractal Higuchi Dimension, and covering a period of 17 years (2001-2017) in a mid-sized city (Focs , ani). The study emphasized the increasing air temperature defined by this index [38].
The main objectives of this study are: (i). to assess the general bioclimatic conditions, based on five bioclimatic indices in 10 of the largest cities of Romania, and (ii). to find out if there has been any change in the bioclimatic indices over a 56 year period (1961-2016) in terms of duration of their occurrence period (DOP) and frequency of occurrence considering the number of days for each class (FO). With this study, we intend to present regional differences in changes in bioclimatic conditions in Romania.

Study Area
Romania is located in South-Eastern Europe, on the north-western shore of the Black Sea (extending on about 9 • on longitude and approximately on 5 • on latitude) and it covers more than 237,000 km 2 (Figure 1). With a general temperate continental climate, there is a large diversity of climate sub-types induced by several important influences: extreme continental in the eastern and southern regions, more humid conditions generated by the moist oceanic air masses originated over the North Atlantic in central and western regions and altitude influenced climate in the mountain region. In general, they lead to specific bioclimatic conditions: milder and more humid in the western regions and hotter in summer and colder in winter and dryer over the entire year in eastern and southern regions. The mean multiannual temperature varies from more than 11.0 • C in the southern regions and on the coastline to sub-zero values in the mountains. Extreme temperatures can reach more than 40 • C in the summertime, and below −20 • C in the wintertime. The largest daily and annual temperature ranges are specific to eastern and southern Romania. In terms of precipitation, the multiannual amount is above 500 mm in western and central regions, but between 300 to 700 mm in Eastern and Southern Romania [39]. In the mountains, the value rises above 1000 mm. Generally, the important difference in climatic conditions between western and central regions on the one hand and the eastern, southern, and southeastern ones on the other hand, are imposed by the presence of the Carpathians, which are considered to be a natural barrier for western moist air masses toward Eastern Europe [8]. In Romania, most studies conducted on climate change focused on air temperature and precipitation and they revealed significant increase in air temperature (mean and extreme temperature indices) and no significant change for precipitation [8,[39][40][41][42][43]. Sunshine hours significantly increased, especially in spring and summer, whereas wind speed significantly decreased in most of the locations considered [41]. No study on changes in the relative humidity or cloudiness has been conducted in Romania so far.
Atmosphere 2020, 11, x FOR PEER REVIEW 3 of 24 the North Atlantic in central and western regions and altitude influenced climate in the mountain region. In general, they lead to specific bioclimatic conditions: milder and more humid in the western regions and hotter in summer and colder in winter and dryer over the entire year in eastern and southern regions. The mean multiannual temperature varies from more than 11.0 °C in the southern regions and on the coastline to sub-zero values in the mountains. Extreme temperatures can reach more than 40 °C in the summertime, and below −20 °C in the wintertime. The largest daily and annual temperature ranges are specific to eastern and southern Romania. In terms of precipitation, the multiannual amount is above 500 mm in western and central regions, but between 300 to 700 mm in Eastern and Southern Romania [39]. In the mountains, the value rises above 1000 mm. Generally, the important difference in climatic conditions between western and central regions on the one hand and the eastern, southern, and southeastern ones on the other hand, are imposed by the presence of the Carpathians, which are considered to be a natural barrier for western moist air masses toward Eastern Europe [8]. In Romania, most studies conducted on climate change focused on air temperature and precipitation and they revealed significant increase in air temperature (mean and extreme temperature indices) and no significant change for precipitation [8,[39][40][41][42][43]. Sunshine hours significantly increased, especially in spring and summer, whereas wind speed significantly decreased in most of the locations considered [41]. No study on changes in the relative humidity or cloudiness has been conducted in Romania so far.

Data Used
The historical data derived from direct observations at 10 weather stations were used to calculate a set of five bioclimatic indices over a 56-year period (1961-2016). The meteorological parameters employed for indices calculation were daily average data for: air temperature (T) (°C), relative humidity (RH) (%), wind speed at 10 m (v10) (m/s), and cloudiness (N) (%). In this research, we also used the daily maximum air temperature (TX) and the daily minimum relative humidity (RHmin).
The weather stations used in the study provide a good spatial coverage for the whole country and include all the climatic regions of Romania. Their location is presented in Figure 1 and their

Data Used
The historical data derived from direct observations at 10 weather stations were used to calculate a set of five bioclimatic indices over a 56-year period (1961-2016). The meteorological parameters Atmosphere 2020, 11, 819 4 of 24 employed for indices calculation were daily average data for: air temperature (T) ( • C), relative humidity (RH) (%), wind speed at 10 m (v10) (m/s), and cloudiness (N) (%). In this research, we also used the daily maximum air temperature (TX) and the daily minimum relative humidity (RHmin).
The weather stations used in the study provide a good spatial coverage for the whole country and include all the climatic regions of Romania. Their location is presented in Figure 1 and their geographical coordinates and elevation are listed in Table 1. Since all the weather stations considered are inside the built area of the cities (except for Craiova weather stations, which is located 200 m away from the built area limit), we consider that they catch weather conditions quite well in the cities low-rise building areas. The climatic data for this study were derived from four main sources. Most of them were provided by the National Meteorological Administration (NMA). The missing values for two of the weather stations were supplemented with data from the existing online databases. The missing data from the Bucharest-Băneasa weather station for N for 2001 were completed with data available on www.meteomanz.com. For the Cluj-Napoca weather station, the T values were freely downloaded from the European Climate Assessment and Dataset project (ECA&D) database (non-blend data) [40], and from www.meteomanz.com [44]. For 2016, the values of the parameters of v10, RH and N were downloaded from the databases www.meteomanz.com and www.rp5.ru [45]. Moreover, for all the weather stations and the entire period, the TX values were extracted from ECA&D and www. meteomanz.com databasis [46,47]. While four databases were used, the data homogeneity was assured by the common source of the data: raw SYNOP messages issued by the weather stations considered. Random checking was performed for the common periods where possible.

Biometeorological Indices
In the present study, we analyzed the following bioclimatic indices: the equivalent temperature (TeK), the Effective Temperature (TE), the Cooling Power (H), the Universal Thermal Climate Index (UTCI), and the Temperature-Humidity Index (THI). They are largely used for the bioclimatic conditions assessment in different regions [18,20,23,[47][48][49]. For some cities in Romania, these indicators were calculated and presented in other short studies [50,51].
To calculate all the indices mentioned above, as well as their parameters, the freely available software BioKlima 2.6 (https://www.igipz.pan.pl/bioklima.html) [52] was used. For wind speed conversion from 10 m altitude (v10 m) to 1.2 m (v), we employed the formula for wind speed extrapolation (1), available on: https://websites.pmc.ucsc.edu/~{}jnoble/wind/extrap/ [53]. (1) where: -v is the velocity to be calculated at height z; -z is the height above the ground level for velocity v; -v ref is the known velocity at height z ref ; -z ref is the reference height where v ref is known; -z 0 represents the roughness length in the current wind direction. Below we presented briefly the main information related to the indicators we calculated (focusing on the calculation formula and bioclimatic comfort classes).

TeK Calculation
TeK evaluates the influence of air temperature and water vapor pressure (e) on the human body. This index was introduced by Dufton [54,55] and Bedford described its use [56,57]. For its calculation, Equation (2) was used: Vapor pressure (e) (hPa) was alos derived from air temperature and relative humidity by using the BioKlima software [52] as Equation (3): where T is the air temperature and RH is the relative humidity.

TE Calculation
TE evaluates the common influence of air temperature, wind speed, and the relative humidity of air. The index establishes a relationship between the identical state of the human body thermoregulatory capacity (warm and cold stress sensation perception) and the differing temperature and humidity of the surrounding environment [48]. It was calculated by using two different Equations (4) and (5), depending on the wind speed values [49].
-for v ≤ 0.2 m/s: -at v > 0.2 m/s: where T is the air temperature, v is the wind speed, and RH is the relative humidity.

H Calculation
The H (W/m 2 ) index was calculated according to Hill's empirical Equations (6) and (7), depending on the wind speed [49,[58][59][60]: -at v ≤ 1 m/s: -at v > 1 m/s: Atmosphere 2020, 11, 819 where T is the air temperature and v is the wind speed. For this index, thermal sensations are assessed according to the modified scale developed by Petrovič and Kacvinsky [49,58,59].

UTCI Calculation
UTCI is the most comprehensive index for calculating heat stress in outdoor environment [61,62]. Its values are calculated as a polynomial regression function up to the sixth order and the input data include meteorological (temperature, mean radiation temperature, the pressure of water vapor or relative humidity, and wind speed (at the elevation of 10 m), and non-meteorological (a metabolic rate of 135 W m −2 and a walking speed of 1.1 m s −1 , clothing thermal resistance) data [57,60,63].
The mean radiant temperature was calculated using Equation (8) [63]: where R is the solar radiation absorbed by the outer layer of clothing in a standing man. This was calculated using the statistical SolAlt model [63][64][65][66][67]. They were identified based on Equations (9) and (10) [63].
where T is air temperature, e is the vapour pressure, Tg is the temperature of the ground surface, as in Equations (11)-(13) [63]: for N < 80% and t ≥ 0 for N < 80% and t < 0 where T is the air temperature, and N is the cloudiness.

THI Calculation
THI is an index developed for warm and hot periods of the year and it is the only index used by the NMA in Romania in the national weather forecast and for releasing early warning messages during summertime for heat stress [35].
In this study, it was calculated based on Equation (14), considering daily maximum temperatures and the daily minimum relative humidity.
where TX is the daily maximum temperature and RHmin is the daily minimum relative humidity.

Thermal Comfort and Discomfort Classes
Each of the thermal indices has been described as having a specified number of classes, usually ranging from extreme cold to extreme heat discomfort. In most situations, the classes do not coincide among the indices and this situation makes comparison among indices very difficult. The thermal comfort classes of the indicators used for this study are summarized in Table 2. Table 2. Classes for thermal comfort and discomfort of the used indices.

Thermal Sensation/Index
TeK

Calculating Indices Parameters
Investigating features such as FO and DOP is of great importance to highlight the possible period of a certain thermal comfort/stress class during the year and to identify any change all over the analyzed period. Therefore, for all the indices, first we identified the feature (thermal comfort/discomfort class) of each day and year in the 56-year period. After that, the FO was calculated for each year. FO is the number of days included in each class (a certain comfort/discomfort condition), between the first and the last day of occurrence of the specified class and year.
For the DOP calculation, the first and the last days of occurrence of each class for every year were detected. The DOP was calculated as the total number of consecutive days between the first and the last day of occurrence for each class each year, no matter if all the days belonged to the same class or not.
As a calculation procedure, for warm and hot stress classes, a number from 1 to 366 was assigned to each day of the year, from the 1st of January to the 31st of December. However, for the cold stress classes, the numbering was made differently, so that the DOP, which was calculated as the difference between the last day and the first day of occurrence, could be homogenous and not present long missing periods. For instance, the first three classes of comfort for the H index are missing every year during the summer months, hence, the numbering for these classes started from the 1st of September to the 31st of July, so that the length of the DOP could be as exactly as possible.
All these operations were made using the Macro option in Microsoft Excel 2013 software.

Trend Detection and Mapping Methods
For trend detection in the two features of the considered indies, we used the Mann-Kendall test [70,71] and the magnitude of the trend (the slope) was calculated by employing the Sen's slope method [72]. The significance level was established at 0.05.
All the calculations were performed by employing XLSTAT ProPlus software (Addinsoft, Paris, France) and the spatial distribution of trend types was mapped using ArcMap10.2 software (ESRI, Bucharest, Romania). For most of the cases, the DOP exceeds 150 days/year, with the cool and slightly cool classes exceeding 200 days/year, or even 250 days/year in western Romania; thus, these conditions are mostly accomplished all year round, except for hot summer days, which are characterized by DOP varying from 100 to 150 days/year corresponding to slightly sultry sensations, respectively less than 50 days/year (apart from Constant , a) for sultry conditions. However, the DOP values are considerable higher for the extra-Carpathians regions, especially for those located in plain and low hilly or tableland areas. The lowest values characterize the center of the country (Cluj-Napoca and Sibiu weather stations) ( Figure 2).

Changes Detected in TeK Index Parameters
Increasing trends were dominant for both parameters considered. Most of them were found to be statistically insignificant, when all the data sets were considered: 53% of the FO series and 72% of the DOP (Figure 3a,c).
The analysis revealed that the FO by comfort/discomfort classes increased for half of the locations, and for 20% of them, the increase was found to be statistically significant (Figure 3a). They were specific mainly to the slightly sultry and sultry conditions classes for eastern and western regions (Figures 3b  and 4). Out of the total data series, 40% indicated downward trends and 17% were characterized by statistical significance (Figure 3a). They were specific mainly to the cold stress and comfortable sensation classes. For most of the cases, the DOP exceeds 150 days/year, with the cool and slightly cool classes exceeding 200 days/year, or even 250 days/year in western Romania; thus, these conditions are mostly accomplished all year round, except for hot summer days, which are characterized by DOP varying from 100 to 150 days/year corresponding to slightly sultry sensations, respectively less than 50 days/year (apart from Constanța) for sultry conditions. However, the DOP values are considerable higher for the extra-Carpathians regions, especially for those located in plain and low hilly or tableland areas. The lowest values characterize the center of the country (Cluj-Napoca and Sibiu weather stations) (Figure 2).

Changes Detected in TeK Index Parameters
Increasing trends were dominant for both parameters considered. Most of them were found to be statistically insignificant, when all the data sets were considered: 53% of the FO series and 72% of the DOP (Figure 3a,c).
The analysis revealed that the FO by comfort/discomfort classes increased for half of the locations, and for 20% of them, the increase was found to be statistically significant (Figure 3a). They were specific mainly to the slightly sultry and sultry conditions classes for eastern and western regions  3b and 4). Out of the total data series, 40% indicated downward trends and 17% were characterized by statistical significance (Figure 3a). They were specific mainly to the cold stress and comfortable sensation classes. In the case of the DOP, only 15% of the data series indicated significant changes (10% increasing and 5% decreasing). In general, the DOP for cool and cold conditions decreased, whereas for the other classes (slightly cool, comfortable, slightly sultry and sultry) it increased, yet not statistically significant. The only region in Romania where a significant change was detected is the southeastern one, located on the seaside of the Black Sea: the length of comfortable and hot stress periods increased over the 56- and 5% decreasing). In general, the DOP for cool and cold conditions decreased, whereas for the other classes (slightly cool, comfortable, slightly sultry and sultry) it increased, yet not statistically significant. The only region in Romania where a significant change was detected is the southeastern one, located on the seaside of the Black Sea: the length of comfortable and hot stress periods increased over the 56year period. A few stationary trends were detected, but they were not spatially coherent (Figure 4).  In the case of the DOP, only 15% of the data series indicated significant changes (10% increasing and 5% decreasing). In general, the DOP for cool and cold conditions decreased, whereas for the other classes (slightly cool, comfortable, slightly sultry and sultry) it increased, yet not statistically significant. The only region in Romania where a significant change was detected is the southeastern one, located on the seaside of the Black Sea: the length of comfortable and hot stress periods increased over the 56-year period. A few stationary trends were detected, but they were not spatially coherent (Figure 4).

FO and DOP Spatial Distribution
For the considered cities, the extreme heat discomfort class of TE (Hot) had a very low frequency at all stations, at a maximum of 2 days/year. Since it seemed irrelevant, we decided not to include it in this study.
The best represented class for the FO was that corresponding to very cold conditions: 90-140 days/year, as average values, depending on the location. The maximum values reached 140-180 days/year. It was followed by the cool and cold sensation classes, with more than 90 days/year, as average values. According to this index, the comfortable or warm stress sensation classes were less frequent, usually less than 15-20 days/year ( Figure 5).
By analyzing the DOP, we found that the most representative class was that characterized by cold sensation one, which was present throughout the entire year, with more than 300 days/year. It was followed by very cold and cool ones covering about 200 days/year. Based on this index, the comfortable and warm conditions were the least frequent and did not exceed maximum values of 60 days/year ( Figure 5).

Changes Detected in the TE Index Parameters
The analyses revealed that most of the detected trends for the TE index increased, reaching up to 50% of the total data series considered for the FO parameter in each class and 53% of the series for the DOP; statistically significant trends represent more than 30% in both cases (Figure 6a,c). The statistically significant decreasing ones were specific to the FO in the case of very cold and cool days, respectively, for the DOP of very cold conditions. at all stations, at a maximum of 2 days/year. Since it seemed irrelevant, we decided not to include it in this study.
The best represented class for the FO was that corresponding to very cold conditions: 90-140 days/year, as average values, depending on the location. The maximum values reached 140-180 days/year. It was followed by the cool and cold sensation classes, with more than 90 days/year, as average values. According to this index, the comfortable or warm stress sensation classes were less frequent, usually less than 15-20 days/year ( Figure 5). By analyzing the DOP, we found that the most representative class was that characterized by cold sensation one, which was present throughout the entire year, with more than 300 days/year. It was followed by very cold and cool ones covering about 200 days/year. Based on this index, the comfortable and warm conditions were the least frequent and did not exceed maximum values of 60 days/year ( Figure 5).

Changes Detected in the TE Index Parameters
The analyses revealed that most of the detected trends for the TE index increased, reaching up to 50% of the total data series considered for the FO parameter in each class and 53% of the series for Atmosphere 2020, 11, x FOR PEER REVIEW 11 of 24 the DOP; statistically significant trends represent more than 30% in both cases (Figure 6a,c). The statistically significant decreasing ones were specific to the FO in the case of very cold and cool days, respectively, for the DOP of very cold conditions. Significant increases were found in cold, fresh, comfortable and warm conditions for both FO and DOP (Figure 6b,d). While the very cold and cool stress classes had the highest frequencies in terms of the number of days over the considered period, this parameter indicated a downward trend for all the locations and for most of them it was detected to be statistically significant. For the fresh, comfortable and warm sensation classes, the great majority of upward trends were found statistically significant (Figure 7). When attention was paid to the DOP, significant negative trends were detected for all the locations considered during very cold days. For eight of them, the decrease was statistically significant and for Significant increases were found in cold, fresh, comfortable and warm conditions for both FO and DOP (Figure 6b,d).
While the very cold and cool stress classes had the highest frequencies in terms of the number of days over the considered period, this parameter indicated a downward trend for all the locations and for most of them it was detected to be statistically significant. For the fresh, comfortable and warm sensation classes, the great majority of upward trends were found statistically significant (Figure 7). When attention was paid to the DOP, significant negative trends were detected for all the locations considered during very cold days. For eight of them, the decrease was statistically significant and for the other classes, significant increase was specific mainly to the eastern and south-eastern cities of the country. For comfortable and warm conditions, significant increase was specific to western cities. However, in the case of fresh, comfortable and warm classes, increasing trend, statistically significant or not, was dominant. Only a few series indicated no change.

FO and DOP Spatial Distribution
This index presents most of the days in the "middle" classes (cool, slightly cool, neutral and hot sensation), the extreme ones (extremely cold and windy, very cold and very hot) lasting less than 10 days/year. The cold conditions had the highest frequency in Eastern Romania (Constanta, Galati and Iasi). The class that covered the majority of days was the neutral one (80-130 days/year). The other three classes (cool, slightly cool and hot) did not differ very much in regards to FO values: the average ones were not higher than 100 days/year. By classes, the minimum FO values ranged from 50 days/year at Timișoara for cool conditions, to 55 days/year for hot days, at Galați and to 70 days/year for the slightly cool class (Figure 8).
The mean DOP varied between 300 and 350 days/year at most weather stations, for the cool, slightly cool sensation and neutral conditions, assuming that they were present over almost the entire year. For the cold and hot classes, it summed around 100-150 days/year. In some cases, the frequency was much higher: cold days at Iași exceeded 300 days/year. The frequency of the extreme classes was less than 50 days/year (Figure 8).

FO and DOP Spatial Distribution
This index presents most of the days in the "middle" classes (cool, slightly cool, neutral and hot sensation), the extreme ones (extremely cold and windy, very cold and very hot) lasting less than 10 days/year. The cold conditions had the highest frequency in Eastern Romania (Constanta, Galati and Iasi). The class that covered the majority of days was the neutral one (80-130 days/year). The other three classes (cool, slightly cool and hot) did not differ very much in regards to FO values: the average ones were not higher than 100 days/year. By classes, the minimum FO values ranged from 50 days/year at Timis , oara for cool conditions, to 55 days/year for hot days, at Galat , i and to 70 days/year for the slightly cool class (Figure 8).

Changes Detected in the H Index Parameters
The trend type for which FO had the highest share was the significantly decreasing one, covering 31% of the data sets, followed by the significantly increasing one (24% of the data sets). The remaining series were equally distributed between not significant upward and downward trends and no trend (15% for each type) (Figure 9a).
Slightly cool, hot and very hot conditions significantly increased in regards to FO in most of the cities, (e.g., 7 out of 10 for slightly cool conditions) ( Figure 10).
For the DOP series, decreasing trends were dominant (58% of the series). Among them, 37% were detected to be statistically significant (Figure 9c) and they characterized especially the cold stress classes, from extremely cold to neutral (Figure 9b). The increasing trends were detected for classes from cool to very hot conditions, for both parameters. The cold sensation class was characterized by significant decrease for all the locations considered (Figure 9b,d). The mean DOP varied between 300 and 350 days/year at most weather stations, for the cool, slightly cool sensation and neutral conditions, assuming that they were present over almost the entire year. For the cold and hot classes, it summed around 100-150 days/year. In some cases, the frequency was much higher: cold days at Ias , i exceeded 300 days/year. The frequency of the extreme classes was less than 50 days/year (Figure 8).

Changes Detected in the H Index Parameters
The trend type for which FO had the highest share was the significantly decreasing one, covering 31% of the data sets, followed by the significantly increasing one (24% of the data sets). The remaining series were equally distributed between not significant upward and downward trends and no trend (15% for each type) (Figure 9a). The DOP indicated negative slopes, equally significant and insignificant for the extremely cold class. The significant decrease was specific to eastern and southeastern Romania. Additionally, no trend characterized four cities, located mainly in the southwestern part of Romania, but this might be attributed to the lack of data (due to extremely low values, the trend could not be calculated) ( Figure 10).  Slightly cool, hot and very hot conditions significantly increased in regards to FO in most of the cities, (e.g., 7 out of 10 for slightly cool conditions) ( Figure 10). The DOP indicated negative slopes, equally significant and insignificant for the extremely cold class. The significant decrease was specific to eastern and southeastern Romania. Additionally, no trend characterized four cities, located mainly in the southwestern part of Romania, but this might be attributed to the lack of data (due to extremely low values, the trend could not be calculated) ( Figure 10).   For the DOP series, decreasing trends were dominant (58% of the series). Among them, 37% were detected to be statistically significant (Figure 9c) and they characterized especially the cold stress classes, from extremely cold to neutral (Figure 9b). The increasing trends were detected for classes from cool to very hot conditions, for both parameters. The cold sensation class was characterized by significant decrease for all the locations considered (Figure 9b,d).
The DOP indicated negative slopes, equally significant and insignificant for the extremely cold class. The significant decrease was specific to eastern and southeastern Romania. Additionally, no trend characterized four cities, located mainly in the southwestern part of Romania, but this might be attributed to the lack of data (due to extremely low values, the trend could not be calculated) ( Figure 10).

FO and DOP Spatial Distribution
For this index, no extreme heat stress days were identified over the entire period considered, in the analyzed cities of Romania. Only four very strong heat stress days were registered at one single location, and thus, due to the extremely low frequency, the two classes (extreme heat stress and very strong heat stress) were not considered for further analysis.
According to the UTCI, in terms of FO, most days were characterized by no thermal stress (140-180 days/year). The slight cold stress days ranged from 60 days/year, at Galat , i, to 80 days/year, at Sibiu; the number of moderate cold stress days varied between 50 days/year, at Timis , oara, and 80 days/year, at Ias , i, and the moderate heat stress characterized, as average values, 30-60 days/year. The remaining classes covered less than 50 days/year for all the locations (Figure 11).
The DOP reached the maximum length for slight cold stress at all the weather stations, except for Bucharest, reaching up to almost 350 days/year, followed by the no thermal stress class. Moderate cold stress was only recorded over 350 days/year at three stations: Botos , ani, Galat , i and Ias , i. It was present elsewhere for less than 200 days/year. The interval when strong and very strong cold stress days can occur did not exceed 150 days at any station. The strong heat stress and extreme cold stress conditions were found to last no more than 50 days/year; however, extreme cold stress occurred at four weather stations over the entire period analyzed (Figure 11).

Changes Detected in the UTCI Index Parameters
For both FO and DOP of the UTCI index, mainly no change was detected in 41% of the FO series ( Figure 12a) and 42% of the DOP series (Figure 12c). This could be explained by the small number of days attributed to extreme thermal (hot or cold) discomfort conditions. Significant changes were found in 38% of data sets considered for FO, respectively, and in 37% of the DOP. The increasing trends were dominant for the heat stress classes (mainly for strong and moderate heat stress), whereas decreasing trends became dominant for the cold stress conditions, especially in the case of strong and very strong cold classes (Figure 12b,d). The DOP also indicated significant increasing trends for the heat stress classes, as well as for slight and moderate cold stress ones. For the no thermal stress class, the statistically not significant increase indicated a higher rate than the significant one. classes covered less than 50 days/year for all the locations (Figure 11).
The DOP reached the maximum length for slight cold stress at all the weather stations, except for Bucharest, reaching up to almost 350 days/year, followed by the no thermal stress class. Moderate cold stress was only recorded over 350 days/year at three stations: Botoșani, Galați and Iași. It was present elsewhere for less than 200 days/year. The interval when strong and very strong cold stress days can occur did not exceed 150 days at any station. The strong heat stress and extreme cold stress conditions were found to last no more than 50 days/year; however, extreme cold stress occurred at four weather stations over the entire period analyzed (Figure 11).

Changes Detected in the UTCI Index Parameters
For both FO and DOP of the UTCI index, mainly no change was detected in 41% of the FO series ( Figure 12a) and 42% of the DOP series (Figure 12c). This could be explained by the small number of days attributed to extreme thermal (hot or cold) discomfort conditions. Significant changes were found in 38% of data sets considered for FO, respectively, and in 37% of the DOP. The increasing trends were dominant for the heat stress classes (mainly for strong and moderate heat stress), whereas decreasing trends became dominant for the cold stress conditions, especially in the case of strong and very strong cold classes (Figure 12b,d). The DOP also indicated significant increasing trends for the heat stress classes, as well as for slight and moderate cold stress ones. For the no thermal stress class, the statistically not significant increase indicated a higher rate than the significant one. As a spatial distribution, the DOP seemed to increase significantly for the strong heat stress in the western and southern regions of the country; for the other regions, no change was detected. Moreover, the no thermal stress class was found to have mainly a positive trend all over the country, except for the northern-most city (Botoșani), for which an insignificant decrease was detected. The upward trends were significant in the extreme western and eastern cities. Statistically significant negative slopes were specific to most cities for the duration of the periods characterized by strong and As a spatial distribution, the DOP seemed to increase significantly for the strong heat stress in the western and southern regions of the country; for the other regions, no change was detected. Moreover, the no thermal stress class was found to have mainly a positive trend all over the country, except for the northern-most city (Botos , ani), for which an insignificant decrease was detected. The upward trends were significant in the extreme western and eastern cities. Statistically significant negative slopes were specific to most cities for the duration of the periods characterized by strong and very strong cold stress. Only two cities experienced no trends (Oradea and Bucharest) for the strong cold stress class, and just a single one (Oradea) had a decreasing downward trend for the very strong cold stress class ( Figure 13). A statistically significant increasing trend for FO in the case of strong and moderate heat stress classes, as well as for the slight cold stress class was detected for the great majority of the locations considered. The significant decreasing trends were specific for the strong, very strong, and extreme cold stress classes ( Figure 13).
A great number of locations with stationary trends detected in the extreme cold stress class, for both FO and DOP (Figure 13), may be explained by the very low number of days characterized by the given conditions corresponding to both parameters.

FO and DOP Spatial Distribution
According to the THI classification, the most days of the year belong to comfortable sensation class, their average number varying between 140 days, at Craiova, and 150 days, at Sibiu. Chill and cold sensation classes characterized, on average, approximately the same number of days, but they varied with the location from 50 to 70 days/year. Very cold days covered between 30 (Constanța) and 60 days/year (Botoșani, Cluj-Napoca and Iași). The most numerous hot days reached their highest mean value of 50 days/year in Bucharest, and the lowest one, 25 days/year, in Sibiu. The mean multiannual number of very hot days did not exceed 20 days/year, on average, whereas their maximum values went up to 50 days/year. Sultry conditions occurred occasionally (less than 5 days/year) (Figure 14). A statistically significant increasing trend for FO in the case of strong and moderate heat stress classes, as well as for the slight cold stress class was detected for the great majority of the locations considered. The significant decreasing trends were specific for the strong, very strong, and extreme cold stress classes ( Figure 13).
A great number of locations with stationary trends detected in the extreme cold stress class, for both FO and DOP (Figure 13), may be explained by the very low number of days characterized by the given conditions corresponding to both parameters.

FO and DOP Spatial Distribution
According to the THI classification, the most days of the year belong to comfortable sensation class, their average number varying between 140 days, at Craiova, and 150 days, at Sibiu. Chill and cold sensation classes characterized, on average, approximately the same number of days, but they varied with the location from 50 to 70 days/year. Very cold days covered between 30 (Constant , a) and 60 days/year (Botos , ani, Cluj-Napoca and Ias , i). The most numerous hot days reached their highest mean value of 50 days/year in Bucharest, and the lowest one, 25 days/year, in Sibiu. The mean multiannual number of very hot days did not exceed 20 days/year, on average, whereas their maximum values went up to 50 days/year. Sultry conditions occurred occasionally (less than 5 days/year) ( Figure 14). Thermal classes with the longest DOP are those characterized by chill and comfortable sensations, which cover more than 250 days/year, as well as the cold one, which mostly accounted for 150 days/year at each location. The frequency of the hot and cold classes can exceed 100 days/year, while the sultry days were quite rare in the region-they only occurred in four cities for less than 10 days/year.

Changes Detected in the THI Index Parameters
Half of the trends identified for the FO data series were increasing, and among them, 36% were found statistically significant (Figure 15a), whereas the other half was shared between decreasing trends and no trend classes: 34% had a downward trend (and among them 21% are statistically significant) and the remaining series indicated no change. From a more detailed perspective, the significant increase was mainly specific to hot and very hot conditions (90% and 80%, respectively, of the locations), whereas a statistically significant decrease was dominant for very cold (90%) and comfortable (60%) thermal conditions. For sultry, chill and cold sensations, the great majority of the series indicated no significant change (Figures 15b and 16).
Further analysis revealed that 62% of the DOP series indicated an increase and 26% of them were detected to be statistically significant (Figure 15c). They were specific mainly to very hot, hot, comfortable and chill conditions (Figure 15d) and to the eastern cities ( Figure 16). The length of the occurrence period showed no change for extreme thermal conditions for the majority of the locations considered: sultry (90%), cold and very cold (40%) (Figures 15d and 16). Thermal classes with the longest DOP are those characterized by chill and comfortable sensations, which cover more than 250 days/year, as well as the cold one, which mostly accounted for 150 days/year at each location. The frequency of the hot and cold classes can exceed 100 days/year, while the sultry days were quite rare in the region-they only occurred in four cities for less than 10 days/year.

Changes Detected in the THI Index Parameters
Half of the trends identified for the FO data series were increasing, and among them, 36% were found statistically significant (Figure 15a), whereas the other half was shared between decreasing trends and no trend classes: 34% had a downward trend (and among them 21% are statistically significant) and the remaining series indicated no change. From a more detailed perspective, the significant increase was mainly specific to hot and very hot conditions (90% and 80%, respectively, of the locations), whereas a statistically significant decrease was dominant for very cold (90%) and comfortable (60%) thermal conditions. For sultry, chill and cold sensations, the great majority of the series indicated no significant change (Figures 15b and 16).     Further analysis revealed that 62% of the DOP series indicated an increase and 26% of them were detected to be statistically significant (Figure 15c). They were specific mainly to very hot, hot, comfortable and chill conditions (Figure 15d) and to the eastern cities ( Figure 16). The length of the occurrence period showed no change for extreme thermal conditions for the majority of the locations considered: sultry (90%), cold and very cold (40%) (Figures 15d and 16).

Discussion and Conclusions
With this study, we intended to present regional differences in changes in bioclimatic conditions in Romania. The detailed analysis for each index trend points out that, in general, both parameters considered in this study, the frequency and length of occurrence period, indicated a similar pattern for the majority of indices and their thermal conditions classes. Since the cities were randomly distributed across the country, we can conclude that, although aimed to develop a study on a local scale, the similar patterns identified lead to the conclusion that we can talk about a regional scale.
From 1961 until 2016, regardless of the calculation method, its purpose or the number of classes, the trend analysis for FO reveals a shift from cold stress conditions to warm and hot ones for each index. However, the most stressful conditions for hot extremes did not indicate significant change. Under these circumstances, one can say that the climate in the big cities of Romania became milder during the cold season and hotter during the warm periods of the year. Further, all indices, except for TE, indicated a general negative trend in the number of comfortable days in terms of thermal sensation. For the DOP, the conclusion was similar, with a longer occurrence interval during the year for comfortable or warm stress classes, except for the H index.
However, based on the FO analysis, comfort (H, UTCI and THI indices) or even cold stress (TE or TeK indices) conditions, were still dominant in the area. Thus, the increasing trends detected will intensify the hot stress, according to some indices, or on the contrary, they will lead to a more comfortable climate, based on other indices.
Due to the previously obtained results indicating significant changes in temperature, and especially in extreme hot temperatures and sunshine hours in Romania over a similar period with that considered for this study, temperature and sunshine seems to be the triggering factors for the identified changes in bioclimatic indices.
Since the great majority of the weather stations are inside the built areas of the selected cities, they are representative for peripheral urban areas with low-rise buildings and green spaces around them. The general assumption is that, for all the indices, under the impact of the urban heat islands, in the high-rise building or dense building areas specific to central areas and to big neighborhoods, the hot stress intensifies, especially during extreme high temperature conditions (such as heat waves), and the cold stress diminishes during the cold season. All considered cities for this study are mid-size and large cities: from more than 100,000 inhabitants to more than 2,000,000 inhabitans (Bucharest). However, most of them shelter between 200,000 and 350,000 people. In their central areas, we can expect that, during the warm season, under the high-rise building impact and rush hours traffic conditions [73,74], the bioclimatic conditions considerably modify compared to those identified in the peripheral areas and become more stressful in terms of hot stress. For the assessment of the real bioclimatic conditions in all types of local climate zones of the cities, installing urban climate monitoring systems is of crucial importance. They will allow identification of the most critical areas in terms of thermal stress (hot spots) as wells as the most comfortable ones. Moreover, the difference between the results of the indices considered for this study leads to the conclusion that a scientific validation by people's perception of these indices on a national scale for the bioclimate conditions of Romania should become a research priority. After validation of thermal stress classes, the index identified as the most relevant should be used for general biometeorological forecasting and considered for implementation in the early warning system for extreme weather events. In the case that different indices are to be found appropriate for different regions, a regional approach should be considered and implemented by the Regional Weather Forecast Centers for a more efficient protection of the population. Under these circumstances, this study could become an extremely useful tool for local and regional authorities in order to adopt the best adaptation measures in terms of thermal stress in the urban areas considered.