Modeling Diurnal Changes in Land Surface Temperature in Urban Areas under Cloudy Conditions

: Land surface temperature (LST) in urban areas is a dynamic phenomenon a ﬀ ected by various factors such as solar irradiance, cloudiness, wind or urban morphology. The problem complexity requires a comprehensive geographic information system (GIS)-based approach. Our solution is based on solar radiation tools, a high-resolution digital surface model of urban areas, spatially distributed data representing thermal properties of urban surfaces and meteorological conditions. The methodology is implemented in GRASS GIS using shell scripts. In these shell scripts, the r.sun solar radiation model was used to calculate the e ﬀ ective solar irradiance for selected time horizons during the day. The calculation accounts for attenuation of beam solar irradiance by clouds estimated by ﬁeld measurements. The suggested algorithm accounts for heat storage in urban structures depending on their thermal properties and geometric conﬁguration. Computed land surface temperature was validated using ﬁeld measurements of LST in 10 locations within the study area. The study conﬁrmed the applicability of our approach with an acceptable accuracy expressed by the root mean square error of 3.45 K. The proposed approach has the advantage of providing high spatial detail coupled with the ﬂexibility of GIS to evaluate various geometrical and land surface properties for any daytime horizon.


Introduction
Land surface temperature (LST) is one of key quantities defining the urban heat island (UHI) phenomenon during hot, sunny summer days. The sun heats dry, exposed urban surfaces, like roofs and roads, often to temperatures 25 to 50 • C higher than the temperature of the ambient air mass. The LST subsequently increases the air temperature depending on various factors such as thermal properties of urban surfaces and wind conditions. Water and vegetated areas typically have lower surface temperature thus contributing to cooler air temperature [1][2][3].
LST is considered to be a reliable indicator of the UHI, as there is a strong correlation between LST and near-surface air temperature for the heat radiating from the surface to the atmosphere [4][5][6][7]. LST is a key parameter in the physics of land surface processes, combining surface-atmosphere interactions and energy fluxes between the atmosphere and the ground. Properties of urban materials, in particular, solar reflectance, thermal emissivity, and heat capacity influence the LST and subsequently the development of UHI, as they determine how the solar radiation is reflected, emitted, and absorbed [8,9]. The magnitude and spatial extent of surface UHIs may vary substantially throughout the day or year depending on the radiative and thermal properties of urban surfaces and geometry [3,10].
Over the last few decades, LST of large areas has been frequently remotely sensed using thermal satellite sensors such as MODIS, Sentinel 3, ASTER, Landsat 7 ETM+, or Landsat 8 TIRS [11,12]. The spatial resolution of thermal satellite data is usually very low (tens to hundreds of meters) and therefore their use in the detailed intra-urban LST analysis is very limited. The airborne or terrestrial thermal data can be used for monitoring smaller areas with higher spatial resolution [13]. However, these techniques have only limited temporal availability and their use for frequent LST monitoring in urban areas can be very time-consuming and costly.
The increasing availability of high-resolution geospatial data and adequate modeling techniques provide an alternative approach to high-resolution estimation of LST in urban areas. Several studies showed the potential of geographic information system (GIS) tools and 3D city models for the estimation of solar radiation in urban areas [14][15][16][17]. Most available solar radiation tools implemented in GIS can be applied to digital surface models (DSMs) representing the uppermost of the features on land surface including buildings and trees. When using at high spatial resolution, a DSM can approximate even vertical surfaces present in urban areas [18]. The trees often mapped using Light Detection and Ranging (LiDAR) technology also play an important role in accurate estimates of solar radiation and LST in a complex urban environment [19,20].
Recently, Hofierka et al. [21] presented a GIS-based approach to model LST in urban areas with high spatial resolution using the r.sun solar radiation model, a DSM derived from a virtual 3D city model and input data describing varying properties of urban greenery derived from airborne and terrestrial LiDAR data and Sentinel 2 multispectral imagery. The suggested algorithm estimates instant values of LST based on physical principles of the Stefan-Boltzmann law depending on available solar irradiance and other input data describing thermal properties of the material. However, the software implementation is static in time and does not account for thermal storage in urban materials. The main effect of thermal storage is the reduction and shift of the daytime land surface temperature peak, at the expense of heat release in the evening and during the night. The presented clear-sky application is definitely associated with the highest LST values, however, in many cities and regions across the world the UHI phenomenon can be found even under semicloudy and cloudy weather conditions. Therefore the goal of this paper is to present a new GIS-based methodology for diurnal calculation of LST in urban areas under cloudy conditions. Thermal storage effects in urban materials are included in the proposed algorithm via a hysteresis-type function implemented in open-source GRASS GIS [22]. Our calculations are validated by field measurements of LST. Our methodology is applied to the study area located in a midsized city in Central Europe with mild climate conditions.

Theoretical Concept for Modeling Land Surface Temperature
Generally, a built-up area exhibits a variable thermal pattern depending on various thermal properties of urban surfaces and geometric configuration of urban structures. For an insulated land surface, the equilibrium land surface temperature, T s is obtained from the heat balance equation based on the Stefan-Boltzmann law which describes the total power radiated from a body in terms of its temperature [21,23]: where α is the unitless solar reflectance or albedo of the surface varying from 0 to 1, I is the global solar irradiance incident on the surface in W.m −2 , which is the output from the r.sun model, ε is the emissivity of the surface, σ is the Stefan-Boltzmann constant, 5.6685 × 10 −8 W.m −2 K −4 , T s is the equilibrium surface temperature in K, T sky is the effective radiant sky temperature, h c is the convection heat transfer coefficient in W.m −2 .K −1 and T a is the air temperature in K [24]. This equation neglects thermal storage in urban structures. Thermal storage reduces and delays the daytime land surface temperature peak, at the expense of heat release in the evening and during the night. Grimmond and Oke [25] report that the storage heat flux is difficult to measure or model because of complex three-dimensional structure of urban surfaces and the diversity of material types which constitute the urban structures. According to Grimmond and Oke [25] the storage heat flux ∆Q s is determined as the energy balance residual from a direct observation of net radiation Q * , and the convective sensible Q H and latent heat Q E fluxes: Therefore we can assume that the effective solar irradiance I s instantly converted to LST on a specific urban material using Equation (1) can be expressed as follows: Equation (3) assumes that part of incoming solar irradiance is converted to the storage heat flux ∆Q s thus lowering/increasing the amount solar irradiance instantly converted to LST via Equation (1) depending on the sign of the ∆Q s term. When ∆Q s > 0, the land surface material accumulates the heat, thus effectively lowering the amount of solar irradiance converted to LST, when ∆Q s < 0, the land surface material releases the heat, thus effectively increasing the amount of solar irradiance converted to LST. The ∆Q s term can be simulated by various mathematical functions, however, the most common form is a hysteresis-type function suggested by Oke and Cleugh [26]: where t is time. In our solution, the coefficient a 1 indicates the overall strength of the heat storage dependence on solar irradiance I. The coefficient a 2 describes the strength and direction of the phase relations between ∆Q s and I in association with temporal changes in solar irradiance I. When the product of a 2 and ∂I ∂t is positive, the heat storage process is strong and lowers the effective solar irradiance I s instantly converted to LST. When it is negative, the heat is released from the material and thus increases the effective solar irradiance I. The coefficient a 3 is an intercept term that indicates a constant release of heat when solar irradiance I is not available. The function is intuitive and easy to parameterize, however, the actual values of a 1 , a 2 , a 3 coefficients depend on various factors such as land cover, geometric configuration of urban space and even concrete meteorological situation. Grimmond and Oke [25] analyzed numerous studies with this model for various land cover classes and documented that the values of a 1 and a 2 coefficients are positive in the broad range of 0.01 to 0.85 and the values of a 3 coefficient are negative in the range of −12.3 to −79.9 W.m −2 .

Implementation in GRASS GIS
The implementation of the suggested algorithm in GRASS GIS is based on two calculation steps. The first step includes a calculation of solar irradiance under real-sky conditions for all time horizons representing diurnal changes starting from the sunrise and ending at sunset times using the r.sun command in GRASS GIS [27]. In case of cloudy weather conditions, the beam component of solar irradiance is attenuated by clouds. In the r.sun command, this attenuation is expressed via the coeff_bh parameter in a form of raster file [22]. The range of values is within the <0,1> interval, with 0 indicating a complete attenuation of beam irradiance while 1 indicates no attenuation. This requires an estimation of cloudiness for a given moment as a spatially distributed parameter. This can be estimated using high-resolution satellite data or approximately via a direct observation. Then, using shellscripting in GRASS GIS [28], we can compute solar irradiances for each time horizon using a cycle between sunrise and sunset times. In the second step, we calculate the heat storage effects via the calculation of effective solar irradiance (Equation (3)) using the r.mapcalc command, a map algebra tool in GRASS GIS [22].
The key parts of the shell script code are presented in the Supplementary Materials of this paper.
The heat storage effects are represented by a hysteresis-type function (Equation (4)) with a 1 , a 2 and a 3 coefficients that need to be determined in advance. These coefficients can be set as a single value valid for the whole area or as spatially distributed parameters represented by rasters. The selection of time step in the diurnal calculation depends on the user. The script can be also modified to include a more accurate numerical estimation of the solar irradiance derivative ∂I ∂t . The outputs of the second step are diurnal time-series of raster files representing heat storages as well as effective solar irradiances. Finally, the computed effective solar irradiances are then used as input solar irradiance files in the lst.stefan-boltzman.sh script [21] to compute the LST values for each time horizon.

The Study Area and Data Collection
This study focuses on the central part of the Košice City in Eastern Slovakia ( Figure 1) as an example of urban area typical for moderate climate of Central Europe that experiences hot summer days and associated UHI phenomenon [29].
ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 4 of 14 tool in GRASS GIS [22]. The key parts of the shell script code are presented in the Supplementary Material of this paper. The heat storage effects are represented by a hysteresis-type function (Equation (4)) with a1, a2 and a3 coefficients that need to be determined in advance. These coefficients can be set as a single value valid for the whole area or as spatially distributed parameters represented by rasters. The selection of time step in the diurnal calculation depends on the user. The script can be also modified to include a more accurate numerical estimation of the solar irradiance derivative t I   .
The outputs of the second step are diurnal time-series of raster files representing heat storages as well as effective solar irradiances. Finally, the computed effective solar irradiances are then used as input solar irradiance files in the lst.stefan-boltzman.sh script [21] to compute the LST values for each time horizon.

The Study Area and Data Collection
This study focuses on the central part of the Košice City in Eastern Slovakia ( Figure 1) as an example of urban area typical for moderate climate of Central Europe that experiences hot summer days and associated UHI phenomenon [29]. The study area (4 km 2 ) was mapped in 2016 and 2017 during several surveying campaigns [30]. The collected data include photogrammetric, satellite, airborne and terrestrial laser scanning data. In this study, we used a DSM with 0.5 m cell size derived from a 3D city model with a level of detail (LoD) 2 consisting of buildings, urban greenery and terrain ( Figure 2) as well as spaceborne multispectral imagery of the Sentinel 2A satellite to derive albedo values for the LST model [21]. The study area (4 km 2 ) was mapped in 2016 and 2017 during several surveying campaigns [30]. The collected data include photogrammetric, satellite, airborne and terrestrial laser scanning data. In this study, we used a DSM with 0.5 m cell size derived from a 3D city model with a level of detail (LoD) 2 consisting of buildings, urban greenery and terrain ( Figure 2) as well as spaceborne multispectral imagery of the Sentinel 2A satellite to derive albedo values for the LST model [21]. The land surface temperature was measured in 10 locations within the study area with various environmental settings (e.g., roofs, walkways, roads, parks) to reflect the diversity of the study area ( Figure 1). The Comet Pt1000TG7/E temperature probe with a data logger was used to collect the LST data with a measurement accuracy of +−0.15 °C . We have selected 25 June 2020 to measure LST from 6:00 to 19:00 zonal Central European Time. In two locations we measured LST continuously with a 30-s time step, in the other locations with poor permanent accessibility, such as parking lots, walkways, we measured LST every hour continuously for 10 min with the first measurement at 6:00. The data from the instrument's data loggers were processed in the Microsoft Excel software. During the measurements, cloudiness was estimated in four categories (clear, semiclear, semicloudy and cloudy) by visual observation of the presence of beam irradiance on the measuring site.

Input Data
Solar irradiance is the most important input parameter in the LST model (1) presented above. Spatial distribution of solar irradiance in urban areas depends on many factors, such as the Earth's geometry, urban morphology, and atmospheric conditions. These factors can be described by a set of equations creating a complex solar radiation model such as r.sun [27]. The r.sun model is one of the most widely used GIS-based solar radiation models that is implemented in the open-source environment of GRASS GIS [28,31]. The r.sun module in GRASS GIS calculates raster-based maps of The land surface temperature was measured in 10 locations within the study area with various environmental settings (e.g., roofs, walkways, roads, parks) to reflect the diversity of the study area ( Figure 1). The Comet Pt1000TG7/E temperature probe with a data logger was used to collect the LST data with a measurement accuracy of +−0.15 • C. We have selected 25 June 2020 to measure LST from 6:00 to 19:00 zonal Central European Time. In two locations we measured LST continuously with a 30-s time step, in the other locations with poor permanent accessibility, such as parking lots, walkways, we measured LST every hour continuously for 10 min with the first measurement at 6:00. The data from the instrument's data loggers were processed in the Microsoft Excel software. During the measurements, cloudiness was estimated in four categories (clear, semiclear, semicloudy and cloudy) by visual observation of the presence of beam irradiance on the measuring site.

Input Data
Solar irradiance is the most important input parameter in the LST model (1) presented above. Spatial distribution of solar irradiance in urban areas depends on many factors, such as the Earth's geometry, urban morphology, and atmospheric conditions. These factors can be described by a set of equations creating a complex solar radiation model such as r.sun [27]. The r.sun model is one of the most widely used GIS-based solar radiation models that is implemented in the open-source environment of GRASS GIS [28,31]. The r.sun module in GRASS GIS calculates raster-based maps of solar irradiance (W.m −2 ) in all three global solar irradiance components (beam, diffuse and reflected) for clear-sky and also real-sky (cloudy) conditions. The computation accounts for sky obstructions (shadowing) by local land surface features represented by a DSM [22]. Using a high-resolution DSM, the r.sun model can provide a sufficiently accurate estimation of solar irradiance on urban surfaces [21].
The LST model (1) implemented in GRASS GIS as the lst.stefan-boltzman.sh shell script [21] requires an initial estimation of LST, which is based on measured air temperature and weather conditions (sunny, overcast day, etc.). The GRASS GIS implementation assumes that the values of ambient air temperature T a and radiant sky temperature T sky are also known beforehand, e.g., from meteorological stations within the city. The official meteorological measurements of air temperature in the Košice City were used for the ambient air temperature parameter of the model (ranging from 16.6 to 25.5 • C on 25 June 2020). The radiant sky temperature as a single value valid for the whole area can be estimated using one of the available approximation methods published, e.g., in [32], depending on cloudiness.
In this study, we use simple cloudy-sky direct temperature models where T sky = T a − 6 under cloudy conditions [33]. Other temperature models can be used to improve the accuracy of the assessment depending on the available data.
The other input parameters of Equation (1) are spatially distributed raster data layers: solar reflectance (albedo) α, thermal emissivity of the surface ε, and convection heat transfer coefficient h c . Albedo and emissivity values can be estimated from multispectral satellite data [34] or assigned to individual land cover objects based on published data. We used the estimates of albedo and thermal emissivity derived from Sentinel 2 multispectral satellite data presented in [21]. The convection heat transfer coefficient h c is usually the most difficult parameter to estimate. It depends strongly on wind speed and direction, geometry of the building and surroundings objects, roof height above ground level, building material texture (roughness) and surface to air temperature difference [9]. More accurate estimates of h c require complex modeling techniques using 3D city models, data on building materials and 3D wind simulations [21]. Our estimates of h c presented in Table 1 are based on expert judgment supported by our previous research in this area [21]. The modified values reflect a different meteorological situation on 25 June 2020 and modified land cover classes (Figure 3).

Results
Simulating diurnal changes in LST requires a dataset representing geospatial and temporal variations of various input factors. We have prepared a dataset consisting of a DSM with 0.5 m spatial resolution representing urban structures such as buildings and trees, land cover representing thermal properties of urban surfaces, albedo and meteorological data such as cloudiness, ambient temperatures and wind speed describing the meteorological conditions in the city on 25 June 2020.
In the first step of the proposed methodology, we calculated real-sky solar irradiance using the r.sun command in GRASS GIS for the selected time horizons with one-hour time step starting at 5:00 and ending 19:00 Central European Time. The calculation included the shading effects of DSM features. Cloudiness was estimated in 10 locations within the study area ( Figure 1) by visual inspection during field measurements of LST in four categories: clear, semiclear, semicloudy and cloudy ( Table 2). During the day, mostly altocumulus and cirrus-type of clouds were present with frequent changes in available direct sunlight at measurement spots. Less transparent cumulus clouds developed in the afternoon. At 17:00 all locations experienced a rain event. The cloudiness affects the beam component of solar radiation. To describe this effect during the period of LST measurements, we assigned the following values of the coeff_bh parameter: clear=1.0, semi-clear=0.75, semi-cloudy=0.5, cloudy=0.25. As the r.sun module requires this parameter as a raster file, we interpolated the recorded values in 10 locations using the v.surf.rst interpolation module in GRASS GIS to 15 raster files of the coeff_bh parameter used by the r.sun module for the respective time horizon. In the calculation of real-sky irradiance, we used a DSM, albedo derived from the Sentinel 2A satellite data [21] and Linke's coefficient of turbidity of 4.3 suggested by the manual page of the r.sun module for urban areas [22,35].
In the second step, we calculated heat storage terms and effective solar irradiance for each time horizon starting from 6:00 to 19:00. This second step requires raster files of solar irradiance computed in the first step as well as coefficients of the hysteresis-type function (Equation (4)). Based on the range of values reported by Grimmond et al. [36] and Grimmond and Oke [25] we used the coefficients a 1 = 0.1, a 2 = 0.1, a 3 = −30. This implies rather weak hysteresis effects which, in our opinion, better reflect quickly changing conditions such as cloudiness and available beam radiation. The resulting 14 raster files of the effective solar irradiance were then used in the lst.stefan-boltzman.sh script published by Hofierka et al. [21] for calculation of LST for each time horizon. This LST model requires several input data such as albedo, thermal emissivity, convection heat transfer coefficient and ambient air temperature (see Section 2.4).
The methodology presented in Section 2 produces a single raster file for each time horizon representing the spatial distribution of LST in the city. Figure 4 shows a spatial distribution of LST over the study area in three selected time horizons (8:00, 14:00 and 18:00). The spatial distribution of LST at 14:00 ( Figure 4) shows urban areas with high LST during clear-sky conditions and strong LST gradients between urban greenery and condensed built-up areas with roofs in the city center. The LSTs in the morning and late afternoon are much lower and with lower gradients between different land cover classes, however, the LST at 18:00 is affected by the rainfall event at 17:00 (Table 2). These LST maps clearly show that urban greenery areas have much lower and relatively stable LSTs during the day regardless of the amount of solar irradiance.
17 h CD-rain CD-rain CD-rain CD-rain CD-rain CD-rain CD-rain CD-rain CD-rain CD-rain  The graph showing the simulated diurnal changes in LST in 10 selected locations is presented in Figure 5. The colors of the lines are grouped according to land cover classes. The highest values of LST are on roofs followed by asphalt roads, parking lots and walkways. The difference between LST on roofs and air temperature is in the range of 20-25 Kelvins during the midday hours with the The graph showing the simulated diurnal changes in LST in 10 selected locations is presented in Figure 5. The colors of the lines are grouped according to land cover classes. The highest values of LST are on roofs followed by asphalt roads, parking lots and walkways. The difference between LST on roofs and air temperature is in the range of 20-25 Kelvins during the midday hours with the highest amount of solar radiation. This is the case even during a cloudy day with mild air temperatures. Urban greenery areas have LSTs only slightly higher than the air temperature. Due to the inclusion of heat storage modelled by a hysteresis-type function (Equation (4)), LSTs are slightly higher in the afternoon. Diurnal changes in LST are strongly affected by cloudiness (Table 2). In our model, the cloudiness affects the amount of beam irradiance which often presents more than 50% of the global irradiance during clear-sky conditions. Point 8 situated on asphalt walkway (Figure 1) is affected by shadows cast by surrounding buildings leading to lower LSTs in the afternoon ( Figure 5). Point 4 is located in a small park under large trees and therefore it experienced LST values very close to air temperatures almost during the whole day. The simulated LSTs are compared to measured LSTs in Figure 6. In comparison to the simulated LSTs ( Figure 5), the increase of LSTs in the morning is slower indicating that LSTs are more affected by heat storage than assumed in our model. At 17:00, the measured LSTs are affected by a brief rain event that strongly lowered measured LSTs. This event was only weakly represented in the simulation. The reduction of beam radiation only partially explains the strong decrease of LSTs and these are more affected by other physical processes and changes of thermal properties of urban surfaces due to the presence of water. The overall root mean square error (RMSE) of our simulation is 3.45 K and Pearson's correlation coefficient is r = 0.91 (Figure 7). The simulated LSTs are compared to measured LSTs in Figure 6. In comparison to the simulated LSTs ( Figure 5), the increase of LSTs in the morning is slower indicating that LSTs are more affected by heat storage than assumed in our model. At 17:00, the measured LSTs are affected by a brief rain event that strongly lowered measured LSTs. This event was only weakly represented in the simulation. The reduction of beam radiation only partially explains the strong decrease of LSTs and these are more affected by other physical processes and changes of thermal properties of urban surfaces due to the presence of water. The overall root mean square error (RMSE) of our simulation is 3.45 K and Pearson's correlation coefficient is r = 0.91 (Figure 7). event that strongly lowered measured LSTs. This event was only weakly represented in the simulation. The reduction of beam radiation only partially explains the strong decrease of LSTs and these are more affected by other physical processes and changes of thermal properties of urban surfaces due to the presence of water. The overall root mean square error (RMSE) of our simulation is 3.45 K and Pearson's correlation coefficient is r = 0.91 (Figure 7).

Discussion
Diurnal changes in LST are affected by various factors with solar irradiance being a dominant factor. The beam component of solar irradiance is strongly affected by cloudiness and morphological features such as buildings and trees casting shadows. The urban morphology can easily be represented by a 3D city model or a digital surface model. However, cloudiness is much harder to estimate. Even standard meteorological observations of cloudiness at meteorological stations are relatively rough estimates unless more accurate pyrheliometer data on beam solar irradiance are available. Thus the cloudiness estimation plays a major role in the accuracy of the LST model. In the presented case study, we used a simple in situ estimation of cloudiness during the period of LST measurements that helped parameterize the LST model with acceptable results reasonably well approximating the diurnal changes in the measured LST values.
The comparison of LST measurements in 10 locations ( Figure 6) with modeling results presented in Figure 4 and Figure 5 showed that the heat storage model needs a better parameterization

Discussion
Diurnal changes in LST are affected by various factors with solar irradiance being a dominant factor. The beam component of solar irradiance is strongly affected by cloudiness and morphological features such as buildings and trees casting shadows. The urban morphology can easily be represented by a 3D city model or a digital surface model. However, cloudiness is much harder to estimate. Even standard meteorological observations of cloudiness at meteorological stations are relatively rough estimates unless more accurate pyrheliometer data on beam solar irradiance are available. Thus the cloudiness estimation plays a major role in the accuracy of the LST model. In the presented case study, we used a simple in situ estimation of cloudiness during the period of LST measurements that helped parameterize the LST model with acceptable results reasonably well approximating the diurnal changes in the measured LST values.
The comparison of LST measurements in 10 locations ( Figure 6) with modeling results presented in Figures 4 and 5 showed that the heat storage model needs a better parameterization reflecting various thermal properties of urban surfaces. It is apparent especially in the morning hours, during rapid changes in cloudiness (15 h), or a rainfall event. While several papers [25,36] suggested a range of values for coefficients a 1 , a 2 , a 3 associated with specific land cover classes, there is a lack of reliable reference values for individual urban surface materials that could be directly used in this LST model. Ideally, the thermal storage response of each land cover class to solar irradiance over different time periods should be measured and used in a spatially distributed form based on the land cover map of the area. Other factors, such as convective heat transfer coefficient h c or albedo, can also exhibit a diurnal variation according to actual meteorological and physical conditions of the surface materials. Clearly, there is a need for more research to understand the temporal variation of most input parameters.
In the presented methodology, we have neglected the anthropogenic factor. The anthropogenic heat release can increase LST significantly in many urban locations. In general, this factor is important especially in industrial areas and during winter. For example, Feng et al. [37] report about 0.66 • C increase of LST for selected Chinese regions. They also report that this factor is relatively small during summer when solar radiation is strong and relatively large during winter with weak solar radiation. Christen and Vogt [38] estimated the contribution of anthropogenic heat flux for Basel in Switzerland at the annual average of +20 W.m −2 . In our case study, we can expect a relatively small contribution from anthropogenic heat to LST variations since we are considering the summer period with moderate temperatures and our study area contains central parts of the city with mostly residential and commercial (shops, offices, etc.) areas. However, during winter days when the heating infrastructure in the buildings is turned on, this can lead to an underestimation of LST. It should be noted that the anthropogenic component can be incorporated into the proposed methodology by adding anthropogenic heat flux release from buildings and roads to the left-hand side of Equation (3). This would require more detailed information on anthropogenic heat release from traffic on streets and buildings.
In comparison to LST derived from thermal satellite sensors, this methodology has a clear advantage of high spatial and temporal resolution depending mainly on availability of geospatial data such as DSM and land cover. Solar irradiance, a key component in our approach, can be estimated for any daytime moment using the r.sun solar radiation model. The disadvantage of this approach is the dependence on various input parameters describing thermal properties of urban surfaces such as the convective heat transfer coefficient h c , emissivity ε or albedo α. Rough estimates of these input parameters for different land cover classes affects the accuracy of the model. Also, this model is applicable only to daytime calculations of LST when solar radiation is available. Nevertheless, the presented results showed that this methodology is sufficiently flexible and robust to simulate daytime diurnal variation of LST in urban areas at unprecedented high spatial and temporal resolution that cannot be achieved by standard remote sensing techniques.
The implementation in an open-source GIS provides a necessary flexibility for different uses and integration with other geospatial tools and models. Using this model, different land cover scenarios or meteorological conditions can be evaluated. Such analytical tools could greatly help with a formulation of UHI prevention and mitigation measures for many cities.

Conclusions
LSTs in urban areas exhibit strong varying spatial and temporal patterns that need to be studied with sufficient spatial and temporal resolution. Geospatial methods and tools can improve our ability to predict this phenomenon and mitigate the unwanted UHI effects at various scales.
In this study, we developed a GIS-based algorithm for daytime diurnal calculation of LST in urban areas under cloudy meteorological conditions based on the Stefan-Boltzmann law and a heat storage model. The suggested methodology is implemented using shell scripts in open-source GRASS GIS. It uses geospatial solar radiation and map algebra tools, a high-resolution digital surface model representing urban structures, spatially distributed data representing various thermal properties of urban surfaces and meteorological conditions. The key component in our solution is a calculation of real-sky solar irradiance using the r.sun solar radiation model in GRASS GIS for any selected time during the day. The calculation accounts for attenuation of beam solar irradiance by clouds estimated by visual inspection or field measurements. The heat storage model is based on the hysteresis-type function with parameters depending on thermal properties of urban land cover classes.
The results of our LST modeling in the city of Košice showed great spatial and temporal variations that cannot be captured by current thermal satellite sensors. The modeling results were compared to LST measurements in 10 locations across the study area. The acceptable accuracy of our results was confirmed with Pearson's correlation coefficient of 0.91, RMSE = 3.45 K and a diurnal pattern following the real dynamics of the phenomenon. The accuracy of the LST model can be further increased by better estimation of key input parameters such as cloudiness or thermal properties of urban surface used by a heat storage model as well as taking into account a contribution from the anthropogenic heat release component of the LST model.