Bioclimatic conditions of the Lower Silesia region (South-West Poland) from 1966 to 2017

This work analyses the temporal and spatial characteristics of bioclimatic conditions in the Lower Silesia region. The daily time values (12UTC) of meteorological variables in the period 1966–2017 from seven synoptic stations of the Institute of Meteorology and Water Management (IMGW) (Jelenia Góra, Kłodzko, Legnica, Leszno, Wrocław, Opole, Śnieżka) were used as the basic data to assess the thermal stress index UTCI (Universal Thermal Climate Index). The UTCI can be interpreted by ten different thermal classes, representing the bulk of these bioclimatic conditions. Stochastic autoregressive moving-average modelling (ARMA) was used for the statistical analysis and modelling of the UTCI as well as separately for all meteorological components. This made it possible to test differences in predicting UTCI as a full index or reconstructing it from single meteorological variables. The results show an annual and seasonal variability of UTCI for the Lower Silesia region. Strong significant spatial correlations in UTCI were also found in all stations of the region. “No thermal stress” is the most commonly occurring thermal class in this region (about 38%). Thermal conditions related to cold stress classes occurred more frequently (all cold classes at about 47%) than those of heat stress classes (all heat classes at about 15%). Over the available 52-year period, the occurrence of “extreme heat stress” conditions was not detected. Autoregressive analysis, although successful in predicting UTCI, was nonetheless unsuccessful in reconstructing the wind speed, which showed a persistent temporal correlation possibly due to its vectorial origin. We conclude thereby that reconstructing UTCI using linear autoregressive methods is more suitable when working directly on the UTCI as a whole rather than reconstructing it from single variables. Electronic supplementary material The online version of this article (10.1007/s00484-020-01970-5) contains supplementary material, which is available to authorized users.


Introduction
Bioclimatology finds applications in many fields such as climate change (Wu et al. 2019), health research (Bröde et al. 2018), epidemiology (Di Napoli et al. 2018), military (Galan and Guedes 2019), and urban planning or even to determine the attractiveness of tourist places such as coastal towns or health resorts (Błażejczyk and Kunert 2011b;Ge et al. 2017). However, at a time when meteorological weather forecasts can be modelled anywhere on the planet, index for representing comprehensive bioclimatic conditions (de Freitas and Grigorieva 2017). The World Meteorological Organization (WMO) officially promoted the use of UTCI as the most suitable tool for determining bioclimatic conditions at the international symposium in April of 2009(WMO 2009). In the last decade, Polish scientists have used UTCI in different parts of Poland and they confirmed that this index is well suited for describing the local bioclimatic conditions (Błażejczyk et al. 2010(Błażejczyk et al. , 2012(Błażejczyk et al. , 2013Chabior 2011;Kuchcik et al. 2013;Okoniewska and Wie ¸cław 2013;Nidzgorska-Lencewicz 2015;Bryś and Ojrzyńska 2016;Rozbicki 2016, 2018). UTCI has also been used worldwide. In Europe, for example, Novak (2013) investigated its use for regions of the Czech Republic. From the perspective of tourism, researchers have analysed changes of UTCI in Greece, Luxemburg, and Hungary Nastos 2011, 2013;Nemeth 2011). Outside Europe, Bröde et al. (2012b) have made predictions using UTCI in South Brazil, Coutts et al. (2016) in Australia, Maciejczyk et al. (2017) in the Arctic, and Ndetto and Matzarakis (2015) in Tanzania.
In this paper the UTCI was used to determine the biometeorological conditions of the Lower Silesia Voivodeship (south-west Poland). We present a broad statistical analysis of the bioclimatic conditions in terms of the UTCI, and use a stochastic modelling approach to test its reconstruction either as a whole or starting from single meteorological variables that define it. We show how statistics of biometeorological conditions are related to the available synoptic stations in the Lower Silesia region. The significantly strong statistical spatial relationships seem to support the modelling of bioclimatic conditions for those locations that do not have a history of meteorological measurements. However, the stochastic analysis appears to be suitable for modelling the UTCI directly rather than starting from its meteorological components.

Location and available data
The meteorological data for seven synoptic stations ( Fig. 1) used in the article have been provided by the Institute of Meteorology and Water Management (IMGW) by way of the climate R package (Czernecki et al. 2020). Air temperature, water vapour pressure, wind speed, and cloud cover for the years 1966-2017 were used. For the calculation 1200 UTC observations (in total 132,951 observations;18,993 per station) were used (Błażejczyk and Kunert 2011b). All analysed synoptic stations are part of the WMO database and they fulfil WMO standards (Jarraud 2008) including measurements of the wind speed at a height of 10 m above ground level. Leszno and Legnica stations had gaps of measurements in cloud cover for the years 1966-1977 (about 12% of data for these stations). For the UTCI, calculations gaps for the cloudiness were approximated from the sun's duration (Bryś et al. 2019).

Universal Thermal Climate index
The UTCI is a function of air temperature, T a ( • C); wind speed, v (m · s −1 ); water vapour pressure, e (hPa); and mean radiant temperature, T mrt ( • C) (Błażejczyk et al. 2013): UT CI = f (T a, v, e, T mrt ).
( 1 ) The mean radiant temperature was calculated with the Bioklima software (Błażejczyk 1996), which uses the following mathematical relationships: where R is the absorbed solar radiation by the human body (W · m −2 ), I rc is the coefficient reducing convective and radiative heat transfer through clothing, L g is the ground radiation (W · m −2 ), L a is the atmosphere's back radiation (W · m −2 ), s h is the emissivity coefficient for humans (0.95), and σ is the Stefan-Boltzmann constant (5.667· 10 −8 W· m −2 · K −4 ). The absorbed solar radiation (R) was calculated using the SolAlt model based on cloudiness (N [%]) and position of the sun (hSl [ • ]); detailed formulas have been reported by Błażejczyk (2005). The range of limiting conditions in order for UTCI to be applicable are as follows: wind speed between 0.5 and 20 m · s −1 (Novak 2013), air temperature from −50 to 50 • C, mean radiant temperature and air temperature difference T mrt − T a from −30 to 70 • C, and relative humidity that needs to be higher than 5% (Bröde et al. 2012a). In the analysed period, only the wind speed exceeded the prescribed limiting conditions. About 2% (i.e. 2768 obs.) of the data were greater than 20 m · s −1 and about 4% (i.e. 5420 obs.) of the data were lower than 0.5 m · s −1 . When wind data exceeded the upper limiting conditions, the UTCI reached its utmost limits. Therefore, following Novak (2013), it was decided to cut off wind data at the limiting conditions and include them in the analysis. The calculated UTCI values correspond to the conditions of a man walking at a speed of 4 km/h (which is equivalent to metabolic changes of 2.3 MET). Activities like mountain hiking will represent a higher metabolic rate; thus, people would have a slightly higher perception than the UTCI prediction.
The scale of thermal stress classes was taken from the physiological model proposed by Bröde et al. (2012a), which is shown in Table 1.

ARMA model
Several works in the literature have reported on the significant relationships between UTCI and various climatic oscillations (Owczarek 2019). For example, there is a significant link between atmospheric circulations (e.g. North Atlantic Oscillation, NAO) and climatic conditions in Poland (Marsz et al. 2019). Since climatic oscillations have successfully been modelled by means of stochastic models, we propose to use AutoRegressive Moving Average models  (Maidment and et al 1993) to decompose bioclimatic conditions (according to UTCI) into their fundamental linear components (i.e. temporal correlation and noise). AutoRegressive Moving Average (ARMA) models (3) have widely been used in hydrology (Salas et al. 1985;Haltiner and Salas 1988;Maidment and et al 1993). ARMA(p,q) stands for Auto-Regressive, AR(p) and Moving Average, MA(q). Consider a time series of N equally spaced observed data (e.g. temperature) y t , with sample average μ = E(y), statistically stationary (i.e. showing no temporal trends), and temporally correlated. An ARMA(p,q) model may be written to represent such series as: with p autoregressive parameters φ(1), ..., φ(p) and q moving average parameters, θ(1), ...., θ(p). The noise t in Eq. 3 is an uncorrelated Gaussian process with zero mean and unit variance (Maidment and et al 1993). Alternatively, the time series can first be detrended, and then standardised by subtracting the mean and then divided by the standard deviation. The standardised data with the removed trend can then be checked for the most suitable ARMA model by computing the sample autocorrelation function (ACF)   (Brockwell et al. 1991). The sample ACF is defined as : where t and s are indexing the time series and k = t − s is the lag at which the correlation is calculated. The PACF is defined by the equation : where: E y 1 ,y n (Y ) is best mean square predictor of Y. If, for example, ACF (4) is nonzero only at lag zero and PACF (5) tails off at lag 2, then the most suitable model is ARMA(0,2), which is coincident with an MA(2). On the other hand, if PACF is nonzero only at lag 0 and ACF tails off at lag, 1 for example, then the best suitable model is ARMA(1,0), which is AR(1). Nonzero ACF and PACF up to a certain lag will otherwise define an ARMA (p,q) model. The purpose of having fitted a suitable ARMA model is to use it in order to decompose the standardised signal into its correlated and uncorrelated components (Brockwell and Davis 2016). Similarly, the ARMA model can then be used for generating a statistically equivalent time series as well for forecasting. The goodness of the model can be examined for example by using the Akaike Information Criteria (AIC) (Akaike 1974), or better still by the modified AIC called AICc which also accounts for overfitting and overparameterisation (Brockwell et al. 1991). AICc is defined as: where: σ 2 is the maximum likelihood estimator of the noise variance, N is the number of observations, and p, q are For each meteorological component (air temperature, vapour pressure, cloudiness, and wind speed) and the computed UTCI from the daily data, we divided the entire record into homogeneous periods showing an almost linear trend, which could then be easily removed. The next step in the decomposition process was to remove the seasonality by subtracting the monthly mean as computed by averaging out the daily data within the same month. Data were then standardised and the most ARMA(p,q) were estimated using the R software (Team and et al 2013). The estimated autoregressive model was then used for generating new data.

Annual and seasonal variability of UTCI in Lower Silesia for the years 1966-2017
Six out of the seven presented stations show comparable biometeorological conditions. The UTCI annual average values are the highest in Opole in 70% of the times for the analysed period   (Fig. 4).
UTCI values oscillate for the seven stations almost synchronously, which reflects the high Pearson correlation coefficients computed between all stations ( Table 2). The high correlation coefficients confirm the comparability of UTCI oscillations in the studied region. Table 3 shows that UTCI values can become extreme (i.e. according to the classification given in Table 1), reaching high values during winter and low values during summer. However, apart fromŚnieżka station, the occurrence of such extremes is generically very low (see Fig. 3). Additionally, the monthly picture offered by Fig. 2 for all stations confirms that there are only few values representing cold stress in summer and heat stress in winter. The case ofŚnieżka station needs some further discussion which will be had in the next section. Here, it suffices to say that the frequency and yearly distribution of the frequent extreme UTCI conditions is probably ascribable to its geographical location, elevation, and aspect ratio (Figs. 3 and 4).

ARMA results
The results of the ARMA model are presented only for the air temperature in Wrocław for the sake of clarity and synthesis. However, all model results are presented in Table 4. All time series were divided into three almost equally long periods : 1966-1983, 1984-2001, and 2002-2017 to which linear trends were fitted. As far as temperature at Wrocław station is concerned, goodness of fit returned R 2 =0.197, R 2 =0.552, and R 2 =0.732 for the respective aforesaid time periods. The computed linear trends for each period and at each station were thus used to detrend the series (Fig. 5) (Maidment and et al 1993).
The next step in the decomposition process was standardisation, which we obtained by subtracting the daily mean, here assumed to be equal to the monthly one, and then dividing it by the standard deviation of the corresponding month. In doing so, we assumed that daily changes are negligible compared with monthly ones, and the resulting process also allowed removing the effect of seasonality. After standardisation, we reasonably assumed the statistical stationarity and homogeneity within the three time periods and then computed both ACF and PACF for the whole length of the series. Both ACF and PACF are shown in (Fig. 6).
ACF tails off gradually over several lags, where as PACF becomes 0 almost immediately after lag 1. This suggests that an autoregressive model order 1, namely AR(1), could be sufficient to remove all correlations in the signal thereby Missing models of cloudiness are due to gaps in the missing data which did not allow for a stochastic model obtaining the uncorrelated residuals and good correspondence between the modelled and sample data (Fig. 7). For a specific station, the AR(1) scored an AICc value of 0.1871, which suggests a good performing model (Sakamoto et al. 1986). The present model can also be made to generate new data, which must be de-standardised and summed according to the original monthly mean and linear trend in order to obtain the statistically equivalent time series (Fig. 8) to the real data. It is important to appreciate how the predicted temperatures (in red) are all contained within the grey belt, which has a width that is equal to ± 1 standard deviation (σ 1 ). The σ 1 according to so-called three-sigma rule of thumb statistically represent at least 70% of the measured data (Wheeler et al. 1992). Thus, the red signal is clearly just one possible realisation of the stochastic process (Fig. 8).
Whilst the proposed linear stochastic approach was shown to be successful in decomposing all time series into deterministic and noise components, this was not the case for the wind speed. For this variable (Table 4), the decomposition of the time series did not lead to satisfactory results for any ARMA model order that we tested. This may have been caused by the vectorial origin of the variable itself, where changes in wind direction with similar  . 8 Reconstruction (continuous line) and forecast (dotted line) of the UTCI compared with the measured series with σ 1 range (grey area) using Wrocław as example magnitude might have produced ambiguity in the data. This generates some spurious temporal correlation, which cannot simply be removed by linear autoregressive models. The UTCI was then computed from the data (blue line in Fig. 8) with the purpose of comparing it with its stochastic model, or with the one obtained from single stochastic models for each meteorological component. Figure 8 shows that fitting a stochastic model on the UTCI after applying the decomposition steps already described produces a stochastic model that successfully removes all correlations and separates the noisy component from the correlated one (green line). This model was built using all UTCI data up to 2017 and then using 2018 data as a prediction.
We also tried to compute the UTCI starting from the stochastic models of each single meteorological component including the wind speed, of which not all correlations could successfully be removed as already described. This resulted in a fluctuating signal (red line in Fig. 8) very similar to the green one in terms of visual appearance in the time domain. This would also result in similar probability density Fig. 9 Autocorrelation functions for 2 approaches of modeled UTCI compared with the UTCI calculated from the data using Wrocław as example distributions, even though the autocorrelation function might be substantially different (Fig. 9).

Discussion
The average annual UTCI values in the Lower Silesia region range from 2.9 • C in Leszno to 5.6 • C in Opole (Table 3). Slightly higher values can be found in the literature. Rozbicka and Rozbicki (2018) showed that the annual average ranged from 3.5 to 9.3 • C for the period 1998-2015. According to Ma ¸kosza (2013), the average UTCI values in the Lubuskie province change between 6.1 • C in Zielona Góra and Gorzów Wielkopolski and 8.1 • C in Słubice. Błażejczyk et al. (2014) calculated an annual UTCI value of 4.1 • C in Warsaw and 5.7 • C in Kołobrzeg, for the period 1991-2000. However, those studies were carried out over shorter periods and in more recent years, where the influence of ongoing climatic changes may have affected the result, thus explaining the slightly higher values. This notwithstanding, all studies correspond to the same heat stress class of "slight cold stress." The situation forŚnieżka is significantly different, it being a mountain station, where cold stress at low temperatures is considerably increased by the high wind speeds (Bröde et al. 2012a). TheŚnieżka Mountain has the most commonly used hiking trails in Karkonosze and knowledge of the local bioclimate can be useful for supporting future tourism in the area. Due to frequent high wind speeds, the average UTCI value is as low as − 24.7 • C. It is worth stressing that cutting off the wind speed to 20 m· s −1 (Novak 2013), which we imposed in order for the UTCI formula to be applicable, has not affected the interpretation of UTCI in terms of thermal comfort. For those regions, the limiting wind speed of 20 m· s −1 already represents high cold stress classes that require additional protective measures, especially from the wind (Błażejczyk 2011a). This is an important piece of information especially for tourism applications, particularly if such a combination of meteorological events leading to extreme cold happens during the spring and summer seasons. Most of the time, tourists might be unprepared for exposure to such extremes (e.g. as occurred inŚnieżka in May where for the analysed period, 30% of the thermal conditions corresponded to strong cold classes (Fig. 2). These findings may also have applications in analysing the bioclimate in arctic windy conditions (Araźny et al. 2019). Each station shows a weak upward trend for UTCI, which is in line with the air temperature and solar radiation trends that occur in Wrocław for the analysed period (Bryś 2010(Bryś , 2013. All synoptic stations also show UTCI seasonal patterns that reflect the Polish annual climate (Marsz et al. 2016(Marsz et al. , 2019. Two heat waves were recorded for the "very strong heat stress" class in 1994 and in 2015. For example, for Opole in 1994, the heat wave persisted for 6 days, i.e. from July 28 to August 2. In Wrocław, Leszno, and Legnica, the heat wave was also classified as "very strong heat stress" and lasted 4 days, i.e. from July 29 to August 1. In Jelenia Góra and Kłodzko, values were 1-2 • C UTCI below the "very strong heat stress" class threshold. In 2015, a strong heat wave lasted from of August 7-12 in Opole, while it was cooler by 5-6 • C of UTCI in other cities. In both heat waves, Jelenia Góra and Kłodzko had the lowest values of UTCI. The 2015 heat wave was also presented for Warsaw by Rozbicka and Rozbicki (2018) where the period August 3-16 was analysed. The number of "strong heat stress" class events for all stations in the analysed period was 1452, which corresponds to a relative frequency of about 1%. Therefore, events leading to heat stress are rather rare and we can assume that in the studied area there is low risk of thermal comfort disruption caused by thermal classes of the warm type (Bröde et al. 2012a). As far as cold thermal classes are concerned, there were 5613 observations of "extreme cold stress" inŚnieżka (i.e. 29% of the time), and 1 observation in Leszno and Legnica in 1996 (0.005%). "Very strong cold stress" occurred 7288 times (5.5% of all station observations) in 5406 days over 52 years. In terms of individual stations, 457 times (2.5%) occurred in Jelenia Góra, 722 times (4%) in Kłodzko, and 664 times (3.5%) in Legnica and Leszno, 661 times (3.5%) in Opole, 339 times (1.5%) in Wrocław, and 3878 times (20%) iń Snieżka. In summary, "extreme cold stress" represents 4.3% of all observations (mainly inŚnieżka), "very strong cold stress" 4.8%, "strong cold stress" 14.8%, and "moderate cold stress" 22.6%. The "slight cold stress" class occurred with a frequency of 17.4%; "no thermal stress" 29.7%, which was the most common thermal class for the analysed period and area, "moderate heat stress" 4.9%, and "strong heat stress" and "very strong heat stress" 1.5% (Fig. 3). These results are similar to those registered in other parts of Poland (Kuchcik et al. 2013;Błażejczyk et al. 2014;Rozbicka and Rozbicki 2018).
We performed a stochastic modelling of the meteorological variables and of the UTCI, using the results for air temperature in Wrocław station as an example. Removing the trends (Fig. 5) and seasonality was necessary in order to use the ARMA models (Fig. 7). The ARMA models were estimated according to the sample auto-and partial-correlation functions (Fig. 6). All models were successful in removing all correlations from noise for all components except for wind speed (Table 4), as shown for example by the comparison with theoretical quantiles for the temperature signal (Fig. 7). As shown in Fig. 8, the prediction underestimates the actual air temperature which is the result of a very warm year in 2018 compared with the past 52 years. The modelled UTCI signal correctly oscillates within the ± σ 1 variability band as the calculated one (blue line) up to the year 2017 (Fig. 8). However, the model also appears to return lower values for 2018 according to the results described for the temperature above. The higher statistical variability for the year 2018 could not be captured by the stochastic model fitted with only the variability of previous years given the ergodicity assumption preserved by the model. It is also worth recalling here that the linear properties of a fluctuating time series are completely identified not only by its probability density function, but also by its temporal correlation. Indeed, this is precisely the property that escapes from a comparison of the UTCI computed from the individually modelled components and the UTCI directly computed from data or from its stochastic model (Fig. 9). Although the magnitude of the variability might be the same for the two modelling approaches, the autocorrelation function of the UTCI reconstructed from the stochastic modelling of each single variable (Fig. 9, right-hand panel) shows a completely different decaying tail, suggesting "too long" a correlation compared with the UTCI obtained from the data (Fig. 9, central panel) and that from its direct stochastic modelling (Fig. 9, left-hand panel). This phenomenon might well be the result of the way the residual correlation in wind speed propagates and is amplified by the UTCI formula, which uses monomial power law relationships of up to the 6th degree for wind speed and its multiplication with other variables.
The presented approach based on stochastic models that have found broad applications in hydrology (Maidment and et al 1993) and the economy (Milani et al. 2017) is therefore appealing for use in bioclimatology as well. However, autoregressive models have successfully recognised and decomposed linear properties (i.e. temporal autocorrelation) of the time series. Somehow counterintuitively, our results suggest that in order to predict UTCI or generate statistically equivalent data, a direct autoregressive model of the UTCI would be more successful than starting from individual models for the single components that make up the UTCI. As all analysed stations are also spatially correlated (Table 2), stochastic models allow for generating statistically equivalent dataset that can be used to reconstruct the bioclimatic conditions in the region.

Conclusions
This study analysed the bioclimatic conditions of Lower Silesia (SW Poland) for the period 1966-2017. A longterm analysis of the bioclimate of Lower Silesia showed many similarities to other parts of Poland (Kuchcik et al. 2013;Ma ¸kosza 2013;Błażejczyk et al. 2014;Rozbicka and Rozbicki 2018). The most favourable biometeorological conditions occur between April and October. The most common thermal class is "no thermal stress," which is also consistent with cases found in the literature. In the case of windy conditions as inŚnieżka located at the top of the mountain, calculation of the UTCI requires that wind speed be cut off at its uppermost limiting value. Although this may well be questionable as it introduces a strong nonlinearity and discontinuity in the resulting time series, its final effect on UTCI may still be correct if the index falls into those classes representing extreme stress conditions. This was the case for the stations analysed in this study, even though we would still recommend that this forced approach should be verified on a case-by-case basis. Forecasts for the air temperature and UTCI for Wrocław were presented using a comparison of the measured data with those of the corresponding stochastic models. However, the presented ARMA modelling approach has shown itself not to be applicable for wind speed, specifically in successfully removing all temporal correlations from the residuals. Accordingly, the UTCI calculated from the single modelled components failed to reproduce some of the linear statistical properties of the real signal, that is, its temporal autocorrelation. On the contrary, the autocorrelation model of the UTCI signal seemed successful and can thus be used for reconstruction and forecasting purposes. This work paves the way for future studies that aim, for instance, not just at clarifying the propagation of the temporal correlation across the formula for the UTCI, but also at reconstructing the dynamical system governing regional bioclimatic conditions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.