Spatial Estimation of Thermal Indices in Urban Areas — Basics of the SkyHelios Model

Thermal perception and stress for humans can be best estimated based on appropriate indices. Sophisticated thermal indices, e.g., the Perceived Temperature (PT), the Universal Thermal Climate Index (UTCI), or the Physiologically Equivalent Temperature (PET) do require the meteorological input parameters air temperature (Ta), vapour pressure (VP), wind speed (v), as well as the different shortand longtime radiation fluxes summarized as the mean radiant temperature (Tmrt). However, in complex urban environments, especially v and Tmrt are highly volatile in space. They can, thus, only be estimated by micro-scale models. One easy way to apply the model for the determination of thermal indices within urban environments is the advanced SkyHelios model. It is designed to estimate sky view factor (SVF), sunshine duration, global radiation, wind speed, wind direction, Tmrt considering reflections, as well as the three thermal indices PT, UTCI, and PET spatially and temporarily resolved with low computation time.


Introduction
In the era of climate change, a growing interest in information about the local impact of global climate change can be seen in politics, public, as well as in science.To assess the local effect, many questions e.g., about heat related mortality [1] and the urban heat island [2] have been worked on.Another important field of interest is the impact of the local morphology on thermal perception and stress of humans (e.g., [3][4][5][6]).
Thermal conditions for humans are not only dependent on air temperature (T a ), but also on moisture (e.g., in terms of vapor pressure (VP)), wind velocity (v) and the short-and longwave radiation fluxes [4][5][6][7][8][9].While T a and VP are rather inherent in space, v and the different radiation fluxes are severely modified by the environment and, thus, highly volatile in space [10].This especially holds for very inhomogeneous environments like urban environments.
As more and more people live in urban areas and cities suffer the most from local effects of global warming, information is needed there the most.Such information is hard to measure and, thus, can only be derived by numeric modeling.
Thermal comfort at street level can be modified by urban planning through changes in building configuration [5,6,11], surface materials [12] and urban green [9,13].To develop and assess adaptation measures, models are required, which are both fast and easy to use, but comprehensive at the same time.
Currently, the ENVI-met model [14,15] is mostly applied for analyzing the thermal conditions in urban environments.While the prognostic model ENVI-met is quite comprehensive, it is not exactly easy to use.The model does require input in a specific format and some input parameters that are hardly available (e.g., the atmospheric conditions in 2500 m above ground level).Furthermore, ENVI-met is quite slow.Depending on model domain size and parameters, a simulation can require one day of computation time for only some hours of model time.As the domain size is limited to 250 on 250 on 40 grids, only rather small areas of interest can be investigated in high resolution of e.g., 1 on 1 m.
While there are more models with, among others, the intention of calculating urban thermal comfort available or under development (e.g., PALM-4U or UMEP / SOLWEIG [16,17]), none of them do meet all the criteria described above.Therefore, a new model needed to be developed.To meet the criteria, the diagnostic SkyHelios model initially developed as a tool for the rapid estimation of the sky view factor and sunshine duration only [18] needed to be extended by a radiation model, a wind model, as well as the routines for the calculation of the thermal indices like the Perceived Temperature (PT) [19], the Universal Thermal Climate Index (UTCI) [20,21], or the Physiologically Equivalent Temperature (PET) [22].The SkyHelios model together with the extensions is addressed as the "advanced SkyHelios model" in the following.
The purpose of this paper can be summarized in terms of three objectives: • Provide a detailed and comprehensive description of the extensions to the diagnostic SkyHelios model,

•
Describe the models capabilities in the course of a small case study,

•
Indicate opportunities and limitations analyzing the results,

•
Identify the most suitable shading type for the reduction of heat stress.

Materials and Methods
The SkyHelios model is a micro scale model for the calculation of sky view factor (SVF) and shading in complex environments using the graphics engine MOGRE [18,23].The model does run on 64 bit Windows machines.
SkyHelios supports many different spatial input file formats that also may be combined (mostly formats supported by the geospatial data abstraction library/openGIS simple features reference implementation GDAL/OGR [24]).However, some functionality can only be used with vector input files (e.g., polygon shapefiles (3D, or containing a height field) for buildings and point-feature shapefiles for trees.All spatial input should be projected in a metric projection, which needs to be specified by its spatial reference ID (SRID) number.
For modeling thermal conditions for humans in complex environments, apart from SVF, a quite sophisticated radiation model, as well as a wind model is required.Both parts are added to the SkyHelios model forming the "advanced SkyHelios model" and are introduced in more detail below.T a and VP are currently considered to be static in space by the advanced SkyHelios model.However, as the SkyHelios model is completely diagnostic, they are variable in time depending on the model input.

Radiation Modeling
The radiation calcuations in SkyHelios are performed on vector-basis and therefore can be run for any point within the model domain.Most of the parametrizations used equal the ones used in the RayMan model [25,26].The most relevant ones are described in the following section.
In contrast to most other models, the SkyHelios model does not consider the urban surroundings in terms of different surfaces, but in terms of pixels in the fisheye image (compare to Figure 1) weighted by the sine of the distance to the images center.Other parameters like the incident diffuse shortwave radiation at a specific location are estimated for the upper hemisphere at once.

Sky View Factor
All radiation calculations in SkyHelios are based on the local sky view factor (SVF).SVF is the fraction of the visible sky, seen from a certain point ( [27], p. 353).It is dimensionless and ranges from 0 to 1, where 0 means that the sky is totally covered by terrain or obstacles, while 1 stands for a free sky.In SkyHelios, first, a fisheye image is rendered for the current location and elevation by the graphics engine (e.g., Figure 1).SVF is determined by distinguishing transparent and colored pixels in the image.Transparent pixels are counted as free sky, while all others are considered as covered by obstacles.As in the real environment, the Fish-eye is a half-sphere, not all of the pixels should have the same influence on the SVF.Therefore, a dimensionless weighting factor ω proj (Equation (1)) is used to consider the projection that adjusts the impact of a pixel by the sine of the zenith angle ϕ ( • ) This results in a spheric SVF.If the planar SVF is desired, another correction by ω planar (dimensionless) needs to be performed (compare to Equation (2)).It increases the impact of objects close to the ground by the cosine of the azimuth angle (counted from the ground to the top): The accuracy of the sky view factor calculations has been assessed by [28].

Shortwave Radiation
Global radiation (G) consists of the direct solar irradiation (I) and the diffuse shortwave irradiation (D) ( [26,27], p. 14) (all energy flux densities in W m 2 ).Both of them are dependent on several parameters and therefore need to be modeled individually.
For a perfect clear sky condition (with no clouds and no horizon limitation), G can be calculated directly using Equation (3) proposed by [10,29,30].In SkyHelios, Equation ( 3) is also used to derive an initial global radiation (G 0 ): . ( Equation ( 3) requires the initial direct solar radiation I 0 , an energy flux density in W m 2 , the solar zenith angle ϕ ( • ), the actual air pressure pr in hPa, as well as the one at sea level for a standard atmosphere (pr 0 = 1013 hPa) and the Linke turbidity factor T L .

Direct Shortwave Irradiation
According to [29], I can be estimated as a function of I 0 ( W m 2 ), ϕ ( • ), T L (dimensionless), the relative optical air mass r opt (dimensionless), the vertical optical thickness of a standard atmosphere δ opt (dimensionless as well), pr in hPa, and cloud cover cc in octas (0 = clear sky to 8 = overcast sky, Equation ( 4)): This of course only holds for unshaded conditions.Under shaded conditions, I is usually assumed to equal 0 W m 2 .The relative optical air mass (r opt ) can be estimated by (Equation ( 5)) [31].
In Equation (5), β s describes the solar elevation angle in • .Using β s and r opt , δ opt can be estimated following an approach by (Equation ( 6)) [32]: Direct shortwave reflections by the surrounding obstacles are considered evaluating the surfaces lighting factor in terms of the blue channel of the fisheye image (compare to Figure 1).The reflection is estimated for any pixel in the fisheye as the blue color value (considering the orientation of the sun to the target surface) times the direct shortwave incident radiation and is considered to be isotropic.The shortwave reflections are therefore not calculated directly as proposed by e.g., [33][34][35] for it is found to be too time-consuming at run-time, but abbreviated from the graphic engines scenes lighting algorithm results.The directional lighting of the scene thereby considers the orientation of any surface to the light source (the sun) as the cosine of the angle between light source direction (a vector into the sun) and the surfaces normal vector.

Diffuse Shortwave Irradiation
According to [36], D can also be calculated as the sum of isotropic (D iso ) and anisotropic diffuse radiation (D aniso , both in W m 2 ).The isotropic part can be calculated by Equation ( 7): This equation requires the direct solar irradiation assuming a clear sky with no clouds I clear ( W m 2 ).The anisotropic component D aniso ( W m 2 , Equation ( 8)) can be approximated by a similar equation if the sun is visible: For the case of the sun covered by horizon or obstacles, D 0,aniso becomes 0 W m 2 .
For non clear sky conditions, a linear correction according to [36] can be applied (Equation ( 9)).It considers the cloud cover (cc) in octas: For a completely covered sky (cc = 8/8), Ref. [36] proposes a simplified equation: In that case, global radiation can be approximated by scaling the initial global radiation (G 0 ) by 0.28 and the local sky view factor (Ψ S ) [36].
Shortwave diffuse reflections are considered in SkyHelios by estimating the diffuse incident radiation based on Equations ( 7)- (10) and scattering it isotropical scaled by the surfaces albedo.

Longwave Radiation
All surfaces are emitting longwave radiation according to the Stefan-Boltzmann law.It is calculating the longwave radiation flux density P lw ( W m 2 ) emitted by a perfectly black radiator surface (s) at a given surface temperature T s (K).The Stefan-Boltzmann law is modified by including an emission coefficient lw (dimensionless, Equation (11)) to be able to apply for non-perfectly black surfaces ( = 1.0), e.g., for humans with an lw of approximately 0.97 ( [10,37], p. 151): Elements within the equation are the longwave radiation flux density P lw in W m 2 , the Stefan-Boltzmann constant of 5.67 × 10 −8 W m 2 •K 4 (σ), the radiating surface area A s in m 2 , and its surface temperature T s in K.
Equation (11) is the basis for the estimation of both the longwave emissions by the sky (P lwA ), as well as the surroundings.The downwards longwave radiation is calculated based on T a and a virtual surface area determined from VP and modified by cloud cover (cc) if available (Equation ( 12)).
Finding the appropriate value for lw can get very hard when it comes to longwave irradiation emitted by the atmosphere [38].While there are empirical values available for quite some time [39][40][41], the determination of lw for a cloudy sky remains quite complex [38,[42][43][44]: The lw of the current surrounding urban structures (buildings, plants) is derived from the object's color in the fisheye image (compare to Figure 1).The surface temperature is determined iteratively evaluating short-(direct and diffuse) and longwave radiational gain as well as longwave emissions.The local T a thereby serves as an initial guess.Soil heat flux is approximated to be 19% of the surfaces energy balance Q if Q ≥ 0.0 W m 2 .For cases with Q ≤ 0.0 W m 2 , a negative ground heat flux of 0.32 • Q is considered (please refer to [26] for details).

Mean Radiant Temperature
The mean radiant temperature T mrt ( • C) is one of the most important input parameters for all sophisticated thermal indices applied in human-biometeorology.It is an equivalent surface temperature, which summarizes the effect of all the different short-and longwave radiation fluxes [29,37,45,46].
T mrt is defined as the surface temperature of a perfectly black, isothermal surrounding environment, which leads to the same energy balance as the current environment [10,37,47,48].
By including the shortwave radiation flux density P sw,s ( W m 2 ) calculated by the diffuse solar irradiation and the diffuse reflected global radiation D s ( W m 2 ) multiplied by the shortwave absorption coefficient (1.0-Albedo, α abs,s ) into Equation (11), the total radiation flux density to and from a surface, P s ( W m 2 ), e.g., the human body can be calculated by Equation (13): Dividing the environment of a person p into a number n of isothermal surfaces i and considering a projection factor Pr to correct the relative surface size of p and s as well as the clothing clo of p, T mrt in W m 2 can be calculated following the principle of equal radiation fluxes caused by the actual and the reference environment (Equation ( 14)): Solving Equation ( 14) for T mrt results in Equation (15), which is perfectly applicable by numerical micro scale models: Most numerical models in urban biometeorology, however, are using further simplifications [17,25,26].

Wind Modeling
Wind data is estimated in SkyHelios based on a diagnostic wind model based on the principles by [50].It considers up to four stream modifications for each individual obstacle.In agreement to [50], a upwind stagnation zone, a downwind recirculation, a downwind velocity deficite zone, as well as a street canyon vortex can be considered.However, updated parametrizations are implemented to allow for improved precision.

Upwind Stagnation Zone
For most obstacles, an upwind stagnation zone according to [51] is considered.The maximum windwards extension is determined by Equation ( 16): In Equation ( 16), L f describes the maximum windward extension of the stagnation zone, h is the obstacles height, and wi stands for the obstacles width (all units in meters).
Ref. [51] also introduces a modified power law profile to determine the average wind speed component ū (m/s) at any vertical level z for within the front eddy zone containing a factor that reduces the wind component perpendicular to the obstacle (compare to Equation ( 17)): In addition, a vortex zone was included into the front eddy zone [51].Its streamwise length L f v (m) is determined by Equation ( 18): The ellipsoidal vortex zone is then filled by a parametrization gathered from fitting experimental wind-tunnel data [51].The trigonometric Equations ( 19) and ( 20) are used to calculate the average vortex components in the horizontal ( ū, Equation ( 19)) and the vertical ( w, Equation (20), both in the m/s) direction: ū(z) w(z) l v and h v are defined as the current length and height of the vortex in meters.

Downwind Recirculation Zone
Ref. [52] introduced an improved parametrizations for the velocity deficit zone in the lee of an obstacle.The shelter model by [52], which is also adopted in the QUICK-URB model calculates a Gaussian velocity deficit pattern using Equations ( 21)-( 23): u d represents the velocity deficit in the lee of an obstacle in m/s, x, y, and z are the stream coordinates in x-, y-, and z-directions, W the width, and H the height of an obstacle in m, U(H) the mean wind speed at the top of the obstacle in m/s based on the upstream power law profile, and C D is the drag coefficient.Γ is defined as 0.6 • c 2 a .The similarity coordinates η, and ζ, as well as the vertical coefficients c a , and a g are calculated according to the following equations: κ in Equation ( 26) is the von Kármán constant of 5.67 The main advantage of the integration of the shelter model by [52] is the streamwise velocity deficits that can be calculated more accurately than by the original parametrization [53].

Street Canyon Parametrization
As a response to the overestimation of the width of a street canyon vortex by the original parametrization (compare to [54]), a modified street canyon model was implemented in SkyHelios.One significant difference to the original parametrization according to [50] is the transition zone at the ends of the vortex formed by vertical wedges.
The vertical wedges are modified after Equation (28), where d sc represents the distance of a point to the upwind obstacle (in meters) and u rt the wind speed component in rooftop height (m/s): In addition, the criteria to detect a street canyon has been modified according to [54].It is now based on the length of the upstream obstacles recirculation zone (L R ).Whether a street canyon is set or not is distinguished according to Equation (29).l represents the obstacles streamwise length: For the central part of a street canyon, Ref. [54] proposes a streamwise speed modification according to Equation ( 30): Equation ( 30) calculates the street canyon modification from the distance to next wall d nw , the street canyon width d SCw (both in meters) and the function F SC , stated by Equation ( 31): The only new parameter in Equation ( 31) is the crosswind distance to the street canyons center d SCccw in m.
For further details on the parametrizations, please see [55].
The wind field calculated by the given functions most likely contains some divergence.Assuming incompressible air, this divergence has to be minimized in order to get a valid wind field.Mathematically, this is performed by minimizing the functional for the scalar H in Equation (32): In Equation (32), α h and α v are horizontal stability factors in s/m, u, v and w are the stream components in m/s, u 0 , v 0 and w 0 are representing the initial stream components in m/s while dx, dy and dz are the grid spacings in m.

Area of Interest and Data
The advanced SkyHelios model is tested for two study areas in Freiburg, southwest Germany (approximately 47 • 59 N, 7 • 51 E).The place of the old synagogue is located westwards from the inner city of Freiburg between the main university buildings I and II (KG I and KG II), the university library and the theater (compare to Figure 2, bottom left).The Institutes Quarter is a city quarter north of the city center mainly consisting of institute buildings of the Albert-Ludwigs-University Freiburg.It covers an area of approximately 700 m on 500 m (compare to Figure 2, top right).
At the place of the old synagogue, some tree information is also available.The point-feature shapefile available consists of 21 trees.To test the radiation calculation capabilities of the SkyHelios model, the attribute fields "opticalDen", "albedo", and "emissivity" (for the optical density, surface albedo and surface emission coefficient) have been assigned arbitrary values in a wide range.For example, the optical density ranges from 0.06 to 1.0, where 0 is transparent and 1 is opaque.
The urban climate station Freiburg, a meteorological background station [56] at the top of a highrise building within the "Institutes Quarter" provides data that can be used as a roof-top reference.The stations records are covering a 13-year period from 1 September 1999 to 30

Results and Discussion
In spite of there being a lot of results generated during this study, only some of them are presented here.However, the selected results are sufficient for the identification of the most appropriate shading type in terms of reduction of heat stress as well as for demonstrating the capabilities of the advanced SkyHelios model.For a spatial result, the example of PET in a height of 1.1 m on 1 July 2008 at 1:30 p.m. LST (Figure 3) was selected.The index is strongly influenced by both radiation and wind speed [6,57] and thus shows the impact of both.
Looking at Figure 3, the strongest reduction of PET by approximately 8.0 • C (compared to PET 38.3 • C at undisturbed locations) can be found in locations shaded by solid obstacles (mostly buildings).This is mainly caused by a reduction in T mrt of almost 22.5 • C due to the obstacles blocking direct solar radiation.Global radiation is reduced by around 650 W m 2 .This is in quite good agreement with studies for summer days in Freiburg, e.g., [58].Obstacles with transparency are causing a slighter reduction.For example, the trees with an optical density of 0.8 are reducing PET by approx.6.3 • C, those with an optical density of 0.6 are causing a reduction of 5.4 • C. Global radiation is reduced from 857 W m 2 in unshaded locations to 393 W m 2 and 334 W m 2 in the same locations, respectively.The effect on PET thereby strongly depends on wind speed.During a measurement campaign in August 2007, a measuring point below a tree was found to be 4.6 • C (PET) cooler on average than one not located below a tree [59].However, as the shortwave transitivity of the tree is unknown, wind speed is slightly higher and radiation fluxes within complex urban environments are generally quite inhomogeneous in space, results are hard to compare.This is also shown by a PET reduction of 3.0 to almost 10.0 • C was found comparing a measurement point under a tree with four other ones without trees for the 2nd of August 2001 in Freiburg [25].However, that study applied the RayMan model, which does not account for transparency.In a significantly warmer setting, Ref. [60] found a reduction in PET of up to 15.6 • C due to (solid) trees for a summer case in Lisbon.The dependence of the reduction on general thermal conditions is indicated by the maximum reduction found for winter conditions was 2.7 • C in the same study [60].Wind speed in the model domain was quite low within the whole model domain at 1:30 p.m. on the 1st of July 2008 with only 2.4 m s in a height of 10 m above ground level.This causes even lower wind speed at pedestrian level (1.1 m above ground) where wind speed was reduced to 0.9 m s in (spatial, arithmetic) average.However, this comprises wind speed in locations without obstacles input (in the northwest and in the southeast of the combined domain where wind speed is approximately 1.4 m s constantly.Wind speed therefore is mostly found to be lower close to the obstacles, but can also be increased due to corner flow and channeling to up to 3.8 m s in some locations.It thereby needs to be noted that wind speed (in the height of 1.1 m) is only slightly decreased by most of the trees.Only the small trees (e.g., in the center of the Place of the Old Synagogue) cause a stronger reduction in wind speed of around 0.1 m s at that level.However, in some locations, an increase by the same amount can be found due to the air current avoiding the tree crown.
Increased wind speed at building corners does reduce PET by up to 2.0 • C in the results.Stagnation caused by buildings, on the other hand, does lead to an increase in PET of up to 15.4 • C in this example.Both effects are clearly present in the results.However, the reduction of PET by increased wind speed can only be found in quite small areas, while low wind speed increases PET in rather large areas (almost all the red areas in Figure 3 are caused by low wind speed).
In spite of the size and the high resolution of the model domain, the model run time of almost one hour on a below average machine can be considered to be quite fast (e.g., compared to the prognostic model ENVI-met [15]).However, to the authors' knowledge, there is currently no other model available that can calculate similar results for a domain size and resolution considered in this study.The model performance in terms of computation time can therefore hardly be compared.
The low computational effort is partly by utilizing the graphics hardware, but also due to the diagnostic model design that allows for avoiding solving transport equations.This comes at the price of unknown previous conditions resulting in e.g., ground or wall heat storage can only be parameterized assuming rather constant conditions.This will e.g., lead to an underestimation in surface temperature in the time after sunset [26].
Some parts of the advanced SkyHelios model have been validated in the past (e.g., the SVF model [28]).Other parts do equal those of the very well validated RayMan model (e.g., the radiation model without the reflections part [25,26]).Furthermore, the results generated by the advanced SkyHelios model are generally plausible and in basic agreement to findings by other studies (e.g., [6,[61][62][63]).However, further validation of the whole model and comparison to on-site measurements is required to assess the model's accuracy.

Conclusions
Results show that heat stress for humans in urban areas can best be prevented by providing shading and, at the same time, not reducing wind speed too much.While buildings are solid obstacles in terms of radiation and therefore can provide more comfortable conditions in the shade, they are solid obstacles to wind at the same time.This mostly causes a reduction in wind speed leading to more uncomfortable conditions.The most comfortable conditions can be found below the trees.They are providing shade and are hardly causing wind speed reductions at the pedestrian level.This of course only holds as long as they are large enough (the small trees in the center of the place of the old synagogue are hardly casting any shadow, but are causing wind sheltering as the crowns are only slightly above the pedestrian level) and as long their crowns are dense enough.The most suitable shading type is therefore found to be provided by individual large trees with dense crowns.
The advanced SkyHelios model is capable of estimating wind speed and -direction, the mean radiant temperature, as well as the thermal indices PT, UTCI, and PET for large areas of interest in high resolution as shown by the results.This allows for analyzing larger areas of interest in more detail by using average office computers, which makes it a very valuable tool for all users working on spatial and temporal dimensions in the field of human thermal biometeorology.

Figure 1 .
Figure 1.Fisheye image showing the upper hemisphere with trees and buildings as generated by the SkyHelios model in production mode.The colors and opacity correspond to different shortwave albedo, longwave emissivity, direct radiation factor of the surfaces (including direct shortwave reflections).The checkerboard background was added to visualize the objects opacity.
April 2013 in hourly resolution.The parameters used in this study are air temperature (T a ) in • C, vapour pressure (VP) in hPa, global radiation (G) in W m 2 , wind speed (v) in m/s and wind direction (WD) in • .To keep it short, only results for the 1 July 2008 at 1:30 p.m. local standard time are presented here.The 1st of July 2008 was a clear summer day with T a of 28.2 • C, vapour pressure of 17.8 hPa, an incident wind speed of 2.4 m/s in 10 m height from 244 • and global radiation of 887 W m 2 at 1:30 p.m.

Figure 2 .
Figure 2. Screenshot of SkyHelios main window showing a combined model domain consisting of two areas of interest in Freiburg, Southwest Germany: the "Institutes Quarter" and the "Place of the Old Synagogue".

igure 3 .
Physiologically Equivalent Temperature (PET) on 1 July 2008 at 1:30 p.m. in a height of 1.1 m above ground level.The calculations consider both the "Institutes Quarter" (upper right) and the "Place of the Old Synagogue" (lower left) together in one large area of interest.