High-Resolution Modelling of Thermal Exposure during a Hot Spell: A Case Study Using PALM-4U in Prague, Czech Republic

: The modelling of thermal exposure in outdoor urban environments is a highly topical challenge in modern climate research. This paper presents the results derived from a new micrometeorological model that employs an integrated biometeorology module to model Universal Thermal Climate Index (UTCI). This is PALM-4U, which includes an integrated human body-shape parameterization, deployed herein for a pilot domain in Prague, Czech Republic. The results highlight the key role of radiation in the spatiotemporal variability of thermal exposure in moderate-climate urban areas during summer days in terms of the way in which this directly affects thermal comfort through radiant temperature and indirectly through the complexity of turbulence in street canyons. The model simulations suggest that the highest thermal exposure may be expected within street canyons near the irradiated north sides of east–west streets and near streets oriented north–south. Heat exposure in streets increases in proximity to buildings with reflective paints. The lowest heat exposure during the day may be anticipated in tree-shaded courtyards. The cooling effect of trees may range from 4 °C to 9 °C in UTCI, and the cooling effect of grass in comparison with artificial paved surfaces in open public places may be from 2 °C to 5 °C UTCI. In general terms, this study illustrates that the PALM modelling system provides a new perspective on the spatiotemporal dif-ferentiation of thermal exposure at the pedestrian level; it may therefore contribute to more climate-sensitive urban planning.


Introduction
Human thermal environments have attracted the interest of many researchers in recent decades [1], and their conclusions have found applications in a wide range of fields, including public health, medicine, physiology, sport, biometeorology, and engineering [1]. The issue of human thermal comfort in the outdoor urban environment, however, has received relatively little attention until recently [2]. Nevertheless, climate change and burgeoning populations in urban areas have led to conditions under which thermal stress and strain are severely and ever-more-frequently impacting the everyday lives of billions of urban citizens around the globe [3]. The past two decades have therefore seen a rapid acceleration in interest in human thermal comfort in cities, leading to work that has attracted substantial attention on the part of urban climatologists, bio-meteorologists, and urban planners [4]. At the moment, urban planning has come to the point at which it is demanding reliable thermal exposure prediction tools in the interests of health risk assessment and adaptation to climatic change [5][6][7]. The benefits of the application of such research to the well-being and health of urban dwellers have been clearly documented by many studies [8,9]. However, a spatiotemporally complex, scientifically clearly described, reliable, and widely recognized thermal exposure simulation model for urban areas has not yet been made available. Indeed, modelling of thermal exposure in outdoor urban environments remains a topical and very challenging task for urban climate research [10,11]. This paper responds to the challenge by presenting the first results of the PALM modelling system [12], a new and sophisticated micrometeorological model that employs a biometeorology module in the computation of universal thermal climate index (UTCI) [13,14]. It may well contribute to progress in the field. The simulations presented herein were performed for a domain in Prague-Dejvice, Czech Republic. At outset, they were validated against a study by [15].
The particular aims of this investigation were (a) to analyze the spatiotemporal patterns of thermal exposure in a moderate-climate urban area during summer days in the context of key environmental parameters of the human thermal environment-mean radiant temperature (MRT), air temperature (T), humidity (H), and wind speed (v); (b) to provide detailed analyses of the temporal variability of thermal exposure based on the UTCI index in a range of selected public urban places; and finally (c) to discuss further improvements to and application of the model in an analysis of thermal comfort in a complex urban environment and its practical usefulness in urban planning.

Thermal Exposure in Urban Areas
Heat exchange between the human body and the environment is driven by radiation, convection, evaporation, and conduction [1]. The quantity of heat exchanged depends upon four key environment-related parameters associated with thermal exposure: air temperature (T), humidity (H), wind speed (v; more precisely, air velocity), and radiation (r, frequently summarized as the integrating variable of mean radiant temperature (MRT)). Two human-related parameters are also critical: activity and clothing [16]. The interactions between these six factors are intrinsic to human thermal exposure [17], which may be expressed by several thermal/biometeorological indices. Among others, the UTCI [13,14] and PET (physiologically equivalent temperature) [18], in particular, are the most frequently used indices in studies addressing thermal comfort in urban areas [2]. A systematic review of the existing indices can be found in [2].
Within a city, the thermal, moisture, aerodynamic, and radiation conditions are markedly different from those in the countryside [19]. Enhanced heat accumulation in a city results from the prevalence of horizontally and vertically oriented impervious surfaces giving rise to intensive absorption of short-wave and long-wave radiation and its manifold reflections [20]. Furthermore, the enclosed spaces between buildings lead to the absorption of emitted long-wave radiation (in other words, to the reduction in heat loss through radiation). A further source of extra heat is related to human activities in urban environments: anthropogenic heat emissions (e.g., from transportation, industry, and airconditioning), together with the lack of evapotranspiration in the absence of soil and vegetation. All these factors result in a complex spatiotemporal pattern of T, H, v, and MRT in urban areas [11,[21][22][23][24][25][26] and increased human thermal exposure in outdoor urban environments [10,27].
In the past, studies of thermal comfort in urban areas have employed in situ measurements of the environmental parameters T, H, v, and r (usually in an approximate calculation of MRT) to calculate one of the thermal indices [28]. However, as indicated above, results based on in situ measurements lack sensitivity to the enormous spatial variations in human thermal exposure in urban areas [29]. It is therefore difficult to obtain data for the spatial representation of thermal exposure in a whole urban area [10,18]. This problem may be addressed by means of modelling techniques.

Modelling Thermal Exposure in the Urban Environment
During the last two decades, substantial progress in modelling urban microclimate processes has been associated with computational fluid dynamics (CFD) models, e.g., MI-TRAS [30], ENVI-met [31], MUKLIMO_3 [32,33], UrbClim [34], and most recently the PALM modelling system (e.g., [12,35,36]). However, the modelling of urban canopy-layer processes at a fine-scale resolution remains a challenging task, where uncertainty is an inevitable part of the modelling process [10,12].
Difficulties and uncertainties are particularly associated with the heterogeneity of surface characteristics (e.g., roughness, albedo, emissivity, heat capacity, plant canopy, and others) as well as a variety of human activities (e.g., anthropogenic heat and air pollution). Together, these result in high variability among elements of the urban atmosphere (e.g., T, H, MRT, and wind components). The length scales of surface heterogeneities in urban areas lie in the 1-100 m range. Therefore, models with input grids far below the kilometer scale are required to resolve the characteristic atmospheric features triggered by a specific urban environment. In view of this fact, atmospheric modelling of urban areas usually considers a limited number of relevant processes; otherwise, it is performed at coarser resolutions that cannot properly describe the conditions within street canyons. One of the greatest challenges lies in which approach to take for turbulence simulation. This may be formulated as some repeatedly asked question: is it sufficient to use a parameterized model of turbulence (e.g., Reynolds-averaged Navier-Stokes (RANS)) or is it more efficient to utilize a large-eddy simulation model (LES) with explicitly resolved turbulence [37]? Modern, increasing computational capabilities would appear to render the application of more detailed LES-based models more attractive [12].
LES models solve the three-dimensional prognostic equations for momentum, temperature, humidity, and other scalar quantities (such as chemical species). PALM-4U is the first open-source meteorological model based on LES principles, and it includes most urban canopy-related processes. This study utilizes several modules embedded in PALM, in particular, a land surface (LSM), plant canopy (PCM), and radiation model that employs external radiation forcing from the Weather Research and Forecasting (WRF) model in our setting. Other PALM-4U components were also applied: the Cartesian topography, building surface model (BSM, formerly urban surface model-USM; see [15]), the radiative transfer model (RTM) [36], the human biometeorology module (BIO) [38], and online chemistry modules. The setup also utilized online and offline domain-nesting PALM-4U features. A new approach to MRT calculation in RTM (see Section 4.2.1) facilitates a BIO module with more precise calculation of derived biometeorological indices such UTCI (as well as PET), which have been found to be the most suitable for application in obstacle-resolving meteorology models [16]. This provides a new perspective on the spatiotemporal differentiation of the physical level of thermal comfort at the pedestrian level by means of validated model results [15].

Domain Description
The study area is located in the northwest part of the center of Prague, the capital city of the Czech Republic ( Figure 1). The model was configured in two nested domains. The parent domain measures 4000 × 4000 m, height 3000 m, at a resolution of 10 m with vertical stretching from 250 m to 20 m. The purpose of such a domain is to allow proper simulation of convection effects and description of the orography of the area. The nested child domain, used for the next analysis, measures 1440 × 1440 m, height 250 m, at a resolution of 2 m. Its purpose is to model the area of interest at a resolution that permits representation of the processes within street canyons (see [12,15]). The parent domain includes highly complex terrain in its northwest and southwest parts. The child domain may be described as a densely built-up area with specific urban conditions created by a traffic circle (Vítězné náměstí) in combination with two boulevards: Evropská/Čs. Armády running south-toeast and Jugoslávských Partyzánů/Svatovítská running north-to-south ( Figure 1). Moreover, this section of the town has been identified as a "hot-spot" and vulnerable area to high temperatures during heat waves in previous studies (e.g., [39]). The eastern and, partially, the southern parts of the child domain consist of the historical residential areas typical of Prague-Dejvice, with a combination of old and new buildings and a variety of other urban components (such as gardens, parks, and parking lots). Most of the southwestern and northeastern parts of the domain are sparsely built-up. Specific to the area are green intra-blocks with gardens and trees, usually with pervious surfaces. In contrast, the historical centers of Prague usually feature impervious intra-blocks (  The red points are located in the east-west (E1-4) and north-south (N1-4)-oriented streets with or without a few trees; the green points (G1-4) are located in densely planted streets; the blue points indicate densely planted gardens in courtyards (T1-2) and "piazzettas" (R1-4); the pink points are placed on impervious surfaces, often public markets (M1-2), closed courtyards (C1-2), or beer gardens (O1-2); and the yellow points are located on open grass surfaces (P1-2) or on children's playgrounds (P3 is on sand, and P4 is on artificial grass). The vector geodata for the left maps were provided by the ESRI (Environmental Systems Research Institute, Europe NUTS 0 Boundaries).

Meteorological Conditions
The summer episode from 19 to 21 July 2018 was used for this analysis [15]. Synoptically, the weather was influenced by a high-pressure ridge centered over northern Germany on 19 July 2018, moving to the Baltic Sea/northern Poland on 21 July 2018. The daily maximum temperature exceeded 30 °C only on 21 July (30.1 °C); the previous days' maxima were 27.6 °C and 28.1 °C, respectively (the daily regime of measured air temperature is in Figure A10 in Appendix D). The beginning of the period was partially cloudy, largely mid-altitude altostratus forming in the morning and early afternoon of 19 July. The period between the afternoon of 19 July and the late afternoon of 21 July was mostly clear, with some high-altitude cirrus. The late afternoon and evening of 21 July was cloudy, mostly with low-level cumulus. The highest values of relative humidity occurred at night (93%), dropping to 27% during the day. The wind above roof level was northwesterly, light, mostly below 3.2 m·s −1 , and often as low as 1.6 m·s −1 . The peak measured wind speeds of 3.1-4.5 m·s −1 were observed in the afternoon of 19 July 2018. The time of sunset was 19:02 UTC on 19 July 2015, and the time of sunrise and solar noon on 20 July were 03:15 and 11:08 UTC, respectively (for more details, see [15]).

GIS Data
BSM and LSM require the input of a number of material input parameters to resolve the energy balance equations. These detail the surface in terms of albedo, emissivity, roughness length, thermal conductivity, and capacity of the skin layer and the volume in terms of thermal capacity and volumetric thermal conductivity. Similarly, RTM requires detailed descriptions of the albedo and emissivity of surfaces. When operating at the high resolution involved in the test case herein (2 m), ground surfaces and the nature of walls and roofing materials become very heterogeneous. Any bulk parameter setting would therefore be inadequate. Detailed settings of these parameters, wherever possible, were the only option. To obtain these data, a supplementary, on-site data collection campaign was conducted to create a database of the fine geospatial data. This included information concerning the wall, ground, and roof materials; colors; roughness; and window fraction, which were used to estimate surface and material properties. Each surface was described in terms of material category, albedo, and emissivity. Categories were assigned parameters estimated on the basis of surface and storage material composition and thickness. The parameters of LSM surface layers were assigned according to the properties of the corresponding PALM LSM category. BSM materials utilized a special set of categories and their parameters; the first two layers were inferred from the properties of the material near the surface, which may differ from the rest of the volume. Each tree was described by its position and vertically stratified leaf-area density (LAD). The trunk and crown were measured separately; estimated perimeters and heights were available for each, while crowns were categorized by shape. LAD decreases towards the center of the crown as a result of the higher branches-to-leaves ratio that occurs there. Building heights were available from the Prague 3D model, maintained by the Prague Institute of Planning and Development. All descriptions of surfaces, materials, trees, and their properties were collected in GIS formats and then preprocessed into the PALM-4U input files (according to the PALM input data standard) corresponding to the particular domain setup. For details of input data collection and processing, see [15].

PALM-4U Model Setup
The dynamic core of the PALM-4U model was configured with the Wicker and Skamarock 5th-order advection scheme [40] and the multigrid pressure solver [41,42]. The radiative fluxes were adapted from the WRF simulation, and their interactions with the urban canopy layer, including the resolved plant canopy, were modelled by the embedded RTM model [36]. The surface energy balance for the individual surfaces (vegetation, pavement, buildings, and water) was calculated by the LSM and BSM components; see [12]. The dynamic and energy processes arising out of resolved trees and shrubs were modelled by PCM.
The initial and boundary conditions for the parent domain of the PALM-4U simulations were obtained from the WRF model. WRF (version 4.0.3) was run on three nested domains with horizontal resolutions of 9 km, 3 km, and 1 km and at 49 vertical levels. The inner domain contained 84 × 84 grid points. The configuration was standard, but parameterizations were chosen in light of decreasing possible discrepancies that might arise from the boundary conditions. NOAH LSM (National Centers for Environmental Prediction/Oregon State University/Air Force/Hydrologic Research Lab model) and RRTMG (is Rapid Radiative Transfer Model) radiation were used in all simulations. The Yonsei University scheme was chosen for parameterization of the planetary boundary layer (PBL). Apart from all of the above, no other urban parameterization was used in the WRF model (all urban processes were simulated in PALM). MODIS (Moderate-resolution imaging spectroradiometer) land-use categories were used without modification. WRF output data were collected from overlapping runs with 12 h durations, initialized from the Global Forecast System (GFS) operational analyses and predictions. The first six hours of each run served as a spin-up, and the following hours (7)(8)(9)(10)(11)(12) were used for the generation of boundary conditions for offline nesting. As the WRF-supplied air flow lacks a turbulent part, any additional turbulence was added at the border of the parent domain by PALM's embedded Synthetic Turbulence Generator (STG) module in every PALM model time-step.

Mean Radiant Temperature Calculation
MRT is of particular importance in determining thermal exposure. With the new angular discretization scheme in place, RTM version 3.0 also added a method for calculating the mean radiant temperature (MRT) [36]. MRT at an arbitrary point in space is defined as the temperature of an imaginary object for which that object would be in radiative equilibrium with its surroundings, meaning that the absorbed irradiance would be equal to the emitted radiance.
Considering both long-wave (LW) and short-wave (SW) radiant fluxes for a hypothetical object with emissivity and SW albedo , MRT can be derived from the equilibrium equation together with the Stefan-Boltzmann law: where , is the mean SW irradiance of the hypothetical object and , is its mean LW irradiance.
This novel calculation method facilitates the simulation of either a spherical blackglobe thermometer or an arbitrary generic object of specified albedo, emissivity, and shape. The clothed human body constitutes a special case of such a shape in simulation; it is tall and narrow, and it therefore affected by radiation from the sides proportionally more than by radiation from above. The MRT values calculated in RTM in the form of average SW and LW irradiances for the hypothetical human body were then used to compute other elements of human thermal comfort in a newly added biometeorology module (BIO; see [12,38]).
The angular discretization scheme in RTM calculates MRT by means of weighted averaging of SW and LW irradiance from all discretized directions. The weights applied to particular directions were used for parametrization of the shape of the simulated object. RTM currently offers three different built-in body shape parametrizations, selectable by the variable mrt_geom. For mrt_geom = 0, the weights are identical in all directions, simulating a spherical object (a black globe thermometer if = 1 and = 0). The second scheme (mrt_geom = 1) was recommended by the authors of the BIO module, based on their earlier implementation of MRT in the BIO module [38]. In this scheme, the weight formula for zenith angle theta ( ) is given as follows: with the default values of parameters a = 0.88 and b = 0.12 prescribed by the BIO module authors. However, this formula lacks a scientific description and evaluation. Moreover, its geometrical representation is not adequate for simulation of a human being. It is asymmetrical vertically (the highest weight is slightly above the horizon), and directions below ≈10 degrees from nadir have zero weight. These issues necessitated the implementation of another shape parameterization scheme. This scheme (mrt_geom = 2) simulates a spheroid S specified by the vertically oriented symmetry half-axis a and half-axis b. Because the symmetry axis is vertical, the weight for a specific direction depends on only the zenith angle and is proportional to the area of the orthographic projection of S from that direction. The orthographic projection shape is the ellipse with half-axes and b, where b is independent of and the area of is directly proportional to . To calculate , a vertical cross section is taken through the centre of S. In this cross section, S is replaced by ellipse E with half-axes a and b and is the line segment formed by orthographic projection of E from zenith angle . The equation of a generic tangent line to ellipse at point , is given as Points , lie on ellipse ; hence, / + / = 1. The tangent of the angle = − between a horizontal line and tangent is equal to its slope: Coordinates , of the tangent point, given angle , would be as follows: This gives two solutions for two tangents, and , from which the positive for and negative for are selected. Then, is given as The origin-crossing line parallel to tangent is given as The distance between two generic parallel lines is described parametrically as + + = 0 equals | | , so the distance between and its origin-crossing parallel equals and the distance between and is double that: The parameters a and b can be specified using configuration parameter mrt_geom_params; the default values are a = 0.88 and b = 0.12. The results in this study were obtained by means of the ellipsoid scheme (mrt_geom = 2) with default values of a and b.

Universal Thermal Climate Index calculation
Calculation of UTCI is based on the Fiala-UTCI multi-node dynamic of human heat transfer and regulation model [43] and includes a clothing model reflecting behavioral adaptation of clothing insulation by the general urban population in response to actual outdoor temperatures [44]. UTCI is then defined as the isothermal air temperature of the reference condition that would elicit the same dynamic response of the physiological model for any combination of T, H, v, and r [13]. UTCI have been already previously successfully employed in a bioclimate simulation in Prague during a heat wave period at the regional scale [45]. In this study, UTCI with ellipsoidal MRT parameterization set in RTM (see Section 4.2.1, MRT) was calculated using the BIO module [38] implemented in PALM-4U. The analyzed outputs represent 10-min averages of the calculated UTCI. These values were calculated at every time step of the model, which varies in a range from 0.1 to 1.0 s depending on prevailing meteorological conditions.

Model Validation
Although the primary aim of this study lies elsewhere, it is directly relevant to summarize the results of a preliminary validation check made against an existing PALM modelling study of the same area [15]. This was performed by an observation team in situ, checking the model against the actual conditions measured "on the ground." First, the WRF results with atmospheric soundings for Praha-Libuš (WMO ID 11520), the meteorological station closest to our domain of interest, were validated. The profiles modelled generally correspond well with station measurements, with some notable exceptions in the surface layer. This is very important to the study herein because the boundary conditions for the PALM simulations were taken from the lower levels. The model tends to show less of a diurnal range, underestimating stability by night and instability during the day. The differences were typically between 1 and 2 K, with only 12 UTC on 21 July exhibiting a difference of ≈4 K (see Figure 8 [15]). The global radiation from the WRF simulations for the summer campaign show very good agreement with the observations. Only in the noon period on 21 July were the maxima in WRF around 25 W.m −2 higher, quite possibly due to an underestimation of the cloud cover in the WRF simulation (see Figure 8 [15]).
Wall and surface temperatures were measured by FLIR SC660 for the 66 ground (horizontal) and 73 wall (vertical) evaluation points representing various surface types [15]. It is not possible to describe such a number of results in detail, so only the most important points will be presented. Generally, at vertical surfaces, the modelled surface temperature agrees fairly well with the observed one in the case of traditional buildings constructed in brick or building blocks; the differences here were typically between 0 K and 3 K. Modern buildings with complex wall structures, e.g., mosaics of metal, glass, and/or plastic, typically exhibit lower agreement. Specific cases involved buildings with glassy or glass-like material surfaces or buildings with high percentages of window fractions. A glassy surface usually reflects the sky or the opposite building(s). Consequently, the derived surface temperature primarily represents the surface temperature of the reflected object (wall, pavement, tree crown, or sky) and not the observed object itself. This type of building is not possible to validate by this mode of measurement.
The differences between observations and models for horizontal surfaces were usually lower, typically between 0 K and 2 K. The validation results demonstrated that modelling of real grassy surfaces in urban environments constitutes a definite challenge. The energy balance of grassy areas may depend strongly on soil-water content, assumed plant canopy, or their direct connection with deeper soil horizon layers, largely unknown to this study. Also revealed was the close dependency of the modelled surface temperature for surfaces in street canyons upon dual rows of trees and on proper representation of leafarea density (LAD), which was especially challenging where crowns form an umbrellalike covering. The model underestimated street surface temperature by up to 5 K at one study point, although such surfaces were modelled well in other tree-shaded locations. Large tree crowns tend to become arranged in clusters with free space between them. At the same time, the modelled cells are homogenous and the clusters are randomly distributed. More detail concerning surface temperature validation may be found in Section 4.3 in [15].
A comparison of the modelled and observed wall heat fluxes at ground-floor level (see Figure 30 [15]) showed similarities in the daily cycles of a similar amplitude. The model slightly overestimated the observed values (by approximately 5-10 W.m −2 ), while the corresponding modelled surface temperature agreed fairly well with the observations. The modelled wall heat flux at first-floor level (see Figure 31 [15]) showed a pronounced daily cycle, while the observed wall heat flux showed only a weak daily cycle of significantly lower amplitude. The modelled surface temperature, however, showed a lower amplitude (higher night-time and lower day-time temperatures) compared with the observations. This behavior suggests an underestimation of the walls' insulation properties for the modelled wall configuration.
Air temperature was validated at three different locations (see Figure 32 [15]). Daytime temperatures fit well, with differences lower than 1 K in general. Only one station, located in Sinkuleho kolej, was temporarily affected by tree canopy in its immediate surroundings, and its differences reached 3 K in maxima. Nighttime temperatures were affected by marked cooling on 19-22 July, which was not captured well by LES, probably due to a misrepresentation of stable conditions [15]. Differences during the coldest part of the night may reach 2.5 K. The final validated variable was wind speed. Wind speeds of under 1 m·s −1 were observed in urban canyons during the measuring campaign (see Figure 33 [15]). The model results generally agree well with the observations.

Spatiotemporal Pattern of Thermal Exposure in Urban Area
For a better understanding of the spatiotemporal pattern of thermal exposure, it is worth considering the map of land-cover categories for Prague-Dejvice classified to PIDS classes and with marked tree positions (see Figure 1).
At the beginning of the simulation, in the late-night hours before sunrise (00:00-03:00 UTC; +02:00 CEST), the typical spatial UTCI/thermal exposure pattern of the "late-night phase" in an urban environment may be observed (Figure 2, left). No thermal stress is predicted in the area of interest (see Figure A9 in Appendix D). Nevertheless, the highest values of UTCI at this stage are usually close to the walls of buildings, arising out of the sensible heat flux and LW radiation emitted from them. Thermal radiation from the walls is at its most intense in narrow alleys and in smaller courtyards with dense tree canopies. Dense trees exert a strong influence on the radiative cooling of surfaces. Conversely, low UTCI values are usually associated with open spaces with few or no trees during the latenight phase (mainly open grass areas, squares, and parking lots). Such low UTCI values are mainly attributable to a lower MRT associated with more rapid radiative cooling due to a large sky-view factor (mainly in cases of clear skies). Moreover, air flow and lower air temperatures play a minor role in reducing MRT, as opposed to UTCI.
During sunrise, the spatial pattern of thermal exposure changes swiftly. Gradually, moderate heat stress is expected to be experienced in more than 50% of the area (see Figure  A9 in Appendix D). At this time, known as the "transient morning phase", the spatial pattern is determined primarily by direct sunlight (i.e., the component of MRT linked to short-wave radiation). In open spaces, this change creates a west-to-east gradient with the highest UTCI values in the immediate vicinity of east-facing walls, with a sharp transition between irradiated and shaded areas (Figure 2 right). UTCI values generally rise rapidly in wider street canyons running east-west, while they remain low in streets with an approximately north-south orientation due to persistent shading. Lower UTCI values in the shadows of the trees also become apparent. UTCI values in courtyards rise slowly due to persistent shading, and these sites exhibit the lowest UTCI values at this time of the day.  Figure A8 in Appendix C. The land-cover categories in Prague-Dejvice, classified by the PALM input data standard (PIDS) classes with tree positions, may be found in Figure 1.
The onset of the third phase in the city, known as the "morning phase" of the spatial thermal exposure pattern, may be observed between 06-07 UTC (+02:00 CEST). At this stage, heat stress is already expected to be experienced in most of the area and the percentage of area under strong heat stress gradually increases. UTCI equalizes progressively in north-to-south and west-to-east streets (Figure 3, left). In north-to-east streets, longwave radiation, the result of higher surface temperature, increasingly contributes to MRT values, while MRT rises rapidly due to direct short-wave radiation. The highest UTCI values are bound to the surroundings of the south-to-east and the sunlit south faces of buildings. Low UTCI values in the shade of vegetation and buildings are already clearly evident. Closed courtyards with a lower proportion of higher vegetation begin to heat up rapidly, leading to a steep rise in UTCI values.
The highest daily averages of UTCI values are reached during the "afternoon" phase of the spatial thermal exposure pattern, which begins between 10 and 11 UTC (+02:00 CEST). It is noticeable that the UTCI peaks occur approximately 2-3 h earlier than peaks in air temperature. At this stage, up to 70% of the area experiences strong heat stress (see Figure A9 in Appendix D). Such high UTCI values reflect the daily maximum of MRT values as a consequence of high surface temperature and its associated emission of longwave radiation, combined with high values of short-wave radiation remaining from direct sunlight. The maximum UTCI values shift from south-facing to southwest-facing walls on buildings within directly irradiated areas, where in some cases, very strong heat stress is predicted (Figure 3, right). However, a strongly developing turbulent flow results in increased spatiotemporal variability of the UTCI pattern. The values within individual courtyards vary widely depending on the quantity and location of their vegetation. Areas of relatively lower UTCI values in the shade of vegetation are clearly evident during the afternoon phase. The lowest UTCI values may be observed on land with larger, compact groups of abundant trees with high, thick crowns (gardens and parks).  Figure A8 in Appendix C. The land-cover categories in Prague-Dejvice, classified by the PIDS classes with tree positions, may be found in Figure 1.
The "evening transition" phase of the spatial thermal exposure pattern in the city typically begins between 14 and 15 UTC (+02:00 CEST). Progressive shading is the dominant factor in spatial UTCI variability. Consequently, the share of area under strong heat stress conditions gradually decreases (see Figure A9 in Appendix D). In most areas, the UTCI pattern exhibits an east-west gradient, with a sharp break at the shadow border (Figure 4 left). The highest UTCI values are bound to neighborhoods in which the walls of the buildings face west. The courtyards shade progressively and occur in localities with relatively lower values of UTCI.
The "early night phase" of the spatial thermal exposure pattern begins to develop between 18 and 19 UTC (+02:00 CEST) immediately after sunset (Figure 4, right). The highest UTCI values are strongly linked to narrow streets and especially to courtyards with a high proportion of impervious surfaces, typically asphalt or concrete; however, no thermal stress is expected to be experienced in most of the area (see Figure A9 in Appendix D). This is mainly due to high MRT values arising out of a combination of sensible heat flux from the buildings, long-wave radiation, and limited ventilation. On the other hand, MRT values in open localities with lawn/grass decrease rapidly; UTCI values in these localities steadily approach those expected for parks and gardens at any given time. The spatial UTCI distribution pattern of the early night phase works towards the pattern characteristic of the late-night phase (see above).  Figure A8 in Appendix C. The land-cover categories in Prague-Dejvice, classified by the PIDS classes with tree positions, may be found in Figure 1.

Streets (Canyons)
Thermal comfort on streets is largely dictated in the daily course of UTCI by street orientation (east-to-west/west-to-east and north-to-south/south-to-north streets) and by the location of the point analyzed (south/north and east/west).
A north-to-south street is characterized by a single daily maximum of UTCI values. In the morning to noon phase, UTCI maxima for the west sides of streets (near east-facing walls; points N1 and N3) arise before noon and the maxima for the east sides of streets (near west-facing walls; points N2 and N4) shift to the afternoon. Finally, both sides of the streets show similar UTCI values at noon ( Figure 5). However, the results herein suggest that UTCI maxima on the east sides of streets (near west-facing walls) reach values of up to 5 °C higher than those on the west sides of streets (near east-facing walls); thus, the intensity of solar heating is the same as that of the already-heated surface (warmed by short-wave reflection and long-wave emission during morning hours). For a street running east-to-west, the difference between the north (irradiated) and south (largely shaded) sides of the street are characteristic. On the south sides of streets (near north-facing walls; points E2 and E4), two daily UTCI maxima occur, bounded by short spells of direct insolation-lower in the early morning and higher in the late afternoon ( Figure 5). On the northern side of the irradiated parts of east-to-west streets (near south-facing walls), the UTCI values are substantially higher and there is a single maximum of UTCI following the period of maximum insolation (points E1 and E3; see Figure 5).
In general, the highest UTCI values in urban canyons occur during the afternoon on the insolated north side of east-to-west streets and on the east side of streets running northeasterly. These parts exhibited over 38 °C UTCI for several hours, a figure that generates very high heat stress. However, model simulations indicate that, when direct irradiation ceases, UTCI values on both sides of these streets show relatively small and con-stant differences in MRT and UTCI during the night. At the same time, predicted fluctuations in the daily course of UTCI values during the daytime are generally higher for eastto-west streets than for north-to-south streets. This phenomenon is predominantly attributable to turbulent air flow in an urban canyon: higher air velocity on the irradiated side leads to the ascent of hot air, then cooler and warm air above the roofs mix, and colder air descends on the shaded side (see Figure 6). At the same time, Figure 6 is proof that insulation of a south-oriented wall can modify the vortexes in street (see the following discussion in Section 6.1).

Figure 6.
Ten minute averages of the vertical wind component proving the modification of thermodynamics in an urban canyon (w, positive = upwards) between points E3 and E4 in an east-to-west canyon at 10:00 (left) and 13:30 (right) UTC on 20 July 2018: the view runs east to west (the buildings are grey, while the one on the right has a south-facing wall). Axis units are expressed in meters.
Moving beyond the effects of the side of the street, it is important to note that the results herein further indicate that the localities of points E1 and E3 returned relatively higher values of UTCI than the localities of points E2 and E4. In other words, the thermal and radiative properties of the northern, irradiated part of the street and adjacent buildings have a strong influence on thermal comfort not only in the irradiated part but also in the shaded part of the street. This effect is more pronounced for east-to-west streets than for north-to-south streets. Nevertheless, considering the results for points G1-4 (see UTCI index on Figure A1 in Appendix A and MRT on Figure A7 in Appendix B), a substantial street-level cooling effect may be attributed to the trees, which substantially reduces daily UTCI maxima and increases nighttime UTCI minima only slightly.

Courtyards
Four representative points were selected: two in courtyards with dense trees (T1 and T2) and the other two in largely paved courtyards with little vegetation (C1 and C2). The results for the daytime show that the UTCI values in the courtyards with trees were 5-10 °C lower than those in the paved courtyards. The trees mitigate spells of moderate heat stress substantially, thus avoiding the occurrence of strong heat stress (Figure 7). Maximum UTCI values occurred in the afternoons. However, this tendency may alter with respect to the irradiance of the courtyard. Small and closed courtyards are usually shaded by neighborhood buildings, and insolation is usually lower as a result. Interestingly, during the first part of the night, the UTCI values remained slightly higher in paved courtyards, a response to higher surface temperatures arising out of heat accumulation in the grounds and walls of buildings. Higher levels of long-wave radiation emitted by the tree crowns in comparison with the radiation from the sky contribute to higher MRT and consequently to higher UTCI values in the course of night in a similar fashion to the situation in street canyons (see Section 5.2.1). During the second part of the night, the UTCI differences between paved courtyards and green courtyards with trees vanish, while phases of higher UTCI values in courtyards with trees may appear. Nevertheless, in all cases, UTCI values during the night remain substantially below the level of thermal stress (Figure 7).

Figure 7.
Moving averages (n = 5) of UTCI in selected locations: C points are located in closed courtyards, and T points are located in courtyards with densely-planted trees. The MRT values for the same times appear in Figure A3 in Appendix B.
The daily maxima of the UTCI values in open public places are bounded by relatively long periods of strong or even very strong heat stress (Figure 8). Only in piazzettas partly shaded by trees did the daily UTCI maxima occur before noon and were heat-stress conditions substantially reduced. During the night, no heat-stress conditions were anticipated in any of the public spaces investigated. The results further suggest quite small differences in UTCI (≈2-5 °C) between paved and grassy open spaces (playgrounds). These were more pronounced at the beginning of the simulation period since soil-moisture content decreased progressively in response to continuous hot, dry weather during the simulated period.

Thermal Exposure in a Moderate-Climate Urban Environment
This study demonstrates the spatiotemporal variability of human thermal exposure and its key environmental parameters by means of the newly developed PALM-4U model and its biometeorology module. It utilizes the well-established UTCI index in a pilot study covering a moderate-climate urban area in Prague-Dejvice during a hot spell.
Six phases of spatial pattern of UTCI in urban environments during hot sunny daysnight, sunrise transitional, morning, afternoon, sunset transitional, and early night-are described. Those phases are especially driven by variability in radiative heat, here represented by MRT as a key factor in thermal exposure [45][46][47]. The importance of innovative approaches to MRT, as employed in this study, is thus highlighted. The new shape parameterization scheme described in Section 4.2.1 constitutes a clear, scientific description and is more appropriate to representation of the human body than previous approaches. It is capable of providing new perspectives on MRT accuracy in micro-scale simulations as well as of providing a valuable source of data for validation in future studies.
The results also emphasize that the spatiotemporal variability of the wind has important effects on thermal exposure in urban areas [48,49]. The effects on thermal exposure of the complexity of turbulence in street canyons [35,50] as well as in cross-ventilated building areas [51][52][53] have been particularly difficult to quantify in the absence of information concerning wind velocity components in street canyons. In urban canyons without trees, only the primary vortex with its leeward side near the overheated surfaces and the windward side near shaded surfaces has been modelled with any regularity, following the description that appears in [54]. However, during typical sunny weather conditions, when the wind speed above a roof falls below 4 m·s −1 , a nonuniform heat flux warms the air near the walls. This may lead to changes in the strength and shape of the primary vortex (see Figure 17 in [35]). The intensity of the vortex changes in response to street orientation, surface materials, and physical properties, including window fraction.
On the other hand, this study suggests that, in moderate-climate urban areas, the heterogeneity of humidity generally has only a minor effect on the spatiotemporal variability of thermal exposure and thermal comfort (for details, see e.g., [55][56][57]). Similarly, the role of the much-addressed spatiotemporal variability of temperature fields in urban areas [58][59][60] is implied to have a comparatively lower effect on thermal exposure in urban areas than MRT; however, further research and analyses are needed.
Focusing more specifically on particular urban structures, the highest thermal exposure is to be expected near the irradiated north sides of east-to-west streets [18,[61][62][63]. However, the results herein suggest that the same level of heat exposure, or even higher, may also be observed around north-to-south streets during particular periods of the day. On the other hand, thermal exposure in courtyards with largely paved surfaces during the daytime is lower than that in open public spaces and compared to exposure in the irradiated parts of streets in the same period. These results need to be further verified and analyzed with respect to the complexity of the specific microclimate of courtyards [64]. It is important to note that heat exposure in streets increases where there is a high percentage of window fraction in the buildings and/or where highly reflective paints predominate on both sides of the street (irradiated and shaded) [64]. It is appropriate to hypothesize that this effect may be attributed to the strong reflection of short-wave radiation. This would then be followed by the absorption of reflected radiation by pavements, which are usually of lower albedo and higher roughness (this needs to be verified by a specially targeted validation campaign).
That trees have a substantial role to play in the reduction of thermal exposure in street canyons and open public places is evident both from the results herein and previous studies. Model simulations estimate the cooling effect of trees within a range from 4 °C to 9 °C UTCI during the day, and they anticipate only minor UTCI increases during the night.
These numbers tally closely with empirical studies conducted in various cities in the central European region [49,64,65]. As a general rule, it could be said that trees decrease thermal stress during the day but may increase it during the night [64,66,67]. Turning to the combination of street orientation and reduction of heat exposure attributable to trees, the model herein lies in agreement with the observations made by [68], anticipating a more marked cooling effect from trees on north-to-south streets than on east-to-west streets. However, the size, shape, height, crown spread, LAD, and the location of trees with respect to the investigated point also play important roles [69][70][71][72].
The presence of grass rather than artificial paved surfaces in open public places leads to less pronounced cooling effects in a range between 2 °C and 5 °C UTCI. These numbers probably remain, however, somewhat overestimated when a comparison is made with empirical studies performed in this region [49,73,74]. The cooling effect of the grass in the urban environments depends heavily on soil moisture and other conditions (e.g., the structure of the soil layers and the LAD). A decrease in the soil moisture, in particular, leads to an increase in surface temperature. This effect has been modelled [15] and observed [75]. Furthermore, on summer days, grass on high-quality soil and subsoil layers typically exhibits higher soil moisture than surfaces with only a shallow, separate soil layer (e.g., non-irrigated green tram-lines) [15].

Model and Simulation Limits, Validation, and Further Development
This study applied the PALM modelling system, SVN revision 4508 [76]. The model was validated [15] and provides similar or slightly higher values than the observations made in comparable central European environments (e.g., [49,77]). Nevertheless, both the model itself and its configuration still have some limitations. An exhaustive description of the model limitations appears in [15]; only those that may affect biometeorology are mentioned here.
First, the results herein demonstrate that urban macro-structure (geometry and orientation) in combination with the radiative and thermal characteristics of surfaces, both horizontal and vertical, have a substantial effect on human thermal exposure in built-up areas [15,64]. However, acquiring data concerning the radiative and thermal characteristics of surfaces in a real urban environment at the required levels of complexity and accuracy within reasonable financial limits constitutes a difficult challenge. We put substantial effort into the collection of detailed data, including the performance of an extensive field campaign, but a certain level of generalization and model uncertainty is still present in the input data. At the same time, the sensitivity tests carried out on the results of the input data allow a level of detail appropriate to modern modelling techniques [64]. It should also be noted that these modelling simulations were primarily validated with respect to air temperature, surface temperature, wind velocity, and concentration of air pollution (NOX and PM10) [15]. Further, the model was configured without PALM's building energy model (BEM), and the inner temperature of buildings was considered constant for the purposes of the simulation. In summer, a constant indoor temperature may certainly be perceived as a simplification, particularly for buildings without air conditioning in which the indoor temperature partially follows the diurnal cycle of the outdoor temperature. It needs to be taken into account that, for the simulated domain, almost no information concerning the quantity and placement of A/C systems was available. On the other hand, the majority of the buildings in the study area are old apartment houses, presumably without central A/C systems. There were no visible individual A/C units on the walls. Thus, it may be deduced that few A/C systems were to be expected in the domain.
The current method of discretization of terrain and buildings in PALM also introduces a limitation. PALM utilizes the Cartesian model grid, actually meaning that the entire volume of each grid cell contains either atmosphere (free air or air with a plant canopy included) or an obstacle (terrain or building). As a consequence, every model surface that is made up of a boundary between obstacle and air has its normal parallel to one of the grid axes. If the modelled domain contains uneven terrain or walls, the discretization may create artificial steps that affect UTCI, the most transparent bioclimatic index. Nevertheless, a comparison with the model results employing PET and other bioclimatic indices apart from UTCI could prove useful.
Regardless of the fact that further direct comparisons would be beneficial, the results of this study demonstrate the substantial potential of the model in the study of spatiotemporal patterns of thermal exposure in urban areas as well as in analyses of the influence of blue, green, and grey features on thermal exposure. The PALM-4U model, with its biometeorology module, covers thermal exposure (i.e., physical level of thermal comfort) in urban areas at an unprecedented level of complexity and resolution. Its application may facilitate substantial progress in research into thermal comfort and resulting heat stress in urban areas. Vital urban planning applications follow from this. However, the quality of the model output depends heavily on the quality of input data. Further long-term interdisciplinary cooperation in developing models for urban climate research is required as well as collaboration in the preparation and organization of the detailed spatial input data essential to most potential applications of the models.
Author Contributions: Coordinate the study, conduct the literature review and prepare the underlying concept of the paper, J.G., M.L.; participate in the model development and test, perform the simulation, and interpret the results, J.R. and P.K.; simulation processing, visualization, and interpretation of the results, J.G. and M.L. design the paper and discuss results, E.S.K. All coauthors participated in the text preparation and its revisions. All authors have read and agreed to the published version of the manuscript.
Funding: Financial support was provided by the Operational Program Prague-Growth Pole of the Czech Republic project "Urbanization of weather forecast, air-quality prediction, and climate scenarios for Prague" (CZ.07.1.02/0.0/0.0/16_040/0000383), which is co-financed by the EU. Coauthor J.G. was supported by the Czech Academy of Sciences under the program for research and mobility support of starting researchers (MSM100302001).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
Input data for a model are available on our archivation server (see http://hdl.handle.net/11104/0315416). PALM SVN code version is cited in litarature and archived in DOI repository.

Acknowledgments:
The terrain-mapping campaign of building properties was co-financed by the Strategy AV21 project "Energy interactions of buildings and the outdoor urban environment", financed by the Czech Academy of Sciences. We would like to thank Jiří Cajthaml and the students of the Faculty of Civil Engineering of the Czech Technical University, Prague, for their help with the terrain-mapping campaign. The PALM simulations, and pre-and postprocessing were performed on the HPC infrastructure of the Institute of Computer Science of the Czech Academy of Sciences (ICS), supported by the long-term strategic development financing of the ICS (RVO:67985807). Parts of the simulations were performed on the supercomputer of the IT4I Czech supercomputing center, supported by The Ministry of Education, Youth, and Sports from the Large Infrastructures for Research, Experimental Development, and Innovations project "IT4Innovations National Supercomputing Center-LM2015070". The WRF and CAMx simulations were performed on the HPC infrastructure of the Department of Atmospheric Physics of the Faculty of Mathematics and Physics of Charles University, Prague, supported by the Operational Program Prague-Growth Pole of the Czech Republic project "Urbanization of weather forecast, air-quality prediction, and climate scenarios for Prague" (CZ.07.1.02/0.0/0.0/16_040/0000383), which is co-financed by the EU. Tony Long (Carsphairn, Scotland) helped work up the English.

Conflicts of Interest:
The authors declare no conflict of interest.