Integrating CFD-GIS modelling to re ﬁ ne urban heat and thermal comfort assessment

We introduce a


Introduction
Cities around the world are under constant change. Population growth is increasing the demand for residential, commercial and traffic areas and thus leading to progressive surface sealing and urban densification (Heilig, 2012;Chapman et al., 2017). Additionally, Earth's climate is undergoing dramatic changes, globally affecting cities by altering local temperature and precipitation patterns, thereby enhancing the occurrence of dry periods and increasing the frequency of excessive heat events and tropical nights per year (Bastin et al., 2019;Zhao et al., 2021;Fischer et al., 2021;Büntgen et al., 2021). These effects impact the urban climate and increase the urban-rural temperature differences, known as the urban heat island (UHI) effect (Oke, 1982), thus, deteriorate intraurban heat and outdoor human thermal comfort (Zhao et al., 2014;Manoli et al., 2019). An increasing time of exposure to an uncomfortable amount of heat can be particularly dangerous and potentially fatal for the young, unhealthy, and elderly people as well as for those performing intense physical work under prolonged exposure (Rydin et al., 2012;Vicedo-Cabrera et al., 2021).
The climate in the urban canopy layer (i.e., the volume beneath building roof level) is directly affected by land surface characteristics, urban geometry and building materials (Stewart and Oke, 2012), thus differing from the urban boundary layer (i.e., the volume above the canopy layer) which is mainly influenced by thermally modified air from advection and convection (Oke, 1976). Therefore, climatic conditions in the urban canopy layer are most relevant for both the indoor (Sailor, 2014) and outdoor (Pigliautile et al., 2020) thermal comfort. Within the urban canopy layer, buildings represent the primary determinant of the urban surface roughness, which, if increased, leads to a reduction in wind speed (Wever, 2012). Despite this, buildings create displacement of the wind flow, leading to a variety of complex inner-city flow patterns such as channelling through street canyons or cross-canyon vortices, a large variability in wind speed patterns and thus significant modifications in microclimatic conditions (Oke et al., 2017).
A variety of models and tools such as RayMan, SOLWEIG, ENVI-met and Ladybug-Grasshopper (Matzarakis et al., 2010;Lindberg et al., 2008;Bruse and Fleer, 1998;Evola et al., 2020) exist to simulate urban heat and heat stress as well as to analyse the effects of specific parameters on urban climatic conditions. To account for important variables that influence urban climate (e.g., surface roughness, wind velocity, air temperature variations), they resort to the capabilities of Geographic Information Systems (GIS) or Computational Fluid Dynamics (CFD) to conduct simulations. Both contrast each other through their strengths and limitations (e.g., including and representing these variables) as described in Aghamolaei et al. (2021).
Applicability of GIS includes mapping of local climate zones (Quan and Bansal, 2021), calculating land surface temperatures from remotely sensed images (Amanollahi et al., 2016), simulating heat stress comprising the effects of radiative fluxes (Chen et al., 2016) or analysing UHI intensities (Nakata-Osaki et al., 2018). Coupling GIS with other models and tools enhances the capabilities for urban heat studies (Equere et al., 2021;Liu et al., 2017) and offers the potential to map human thermal stress based on, e.g., the standardized Physical Equivalent Temperature (PET), as was shown by Koopmans et al. (2020). However, with urban physics and thermodynamics in its repertoire, CFD is gaining popularity in modelling microclimatic conditions in urban areas (Gromke et al., 2015;Toparlar et al., 2015;Antoniou et al., 2019;Blocken, 2015;Toparlar et al., 2017). CFD is a numerical method to simulate and study the hydrodynamics of fluid flows by solving the governing Navier-Stokes equations (Versteeg and Malalasekera, 2007). Recent CFD-based studies include assessments of wind velocity (Zhang et al., 2021), analysis of heat fluxes (Allegrini et al., 2015), characterisation of microclimatic extremes (Javanroodi et al., 2022) or, at a finer scale, simulations of the drag effect of urban trees (Zeng et al., 2020) or thermal behavior of different ground surfaces (Yang et al., 2013). Currently, limitations such as cumbersome computation time with increasing scale, a significant level of expert knowledge to handle the available software, data acquisition and over-simplifications of underlying physics represent apparent drawbacks that prevent realisation of CFD's full potential (Mirzaei and Haghighat, 2010;Broadbent et al., 2019;Mirzaei, 2021).
Although GIS and CFD-based approaches are common in urban climatology for analysing urban heat, thermal comfort and the UHI effect at different spatial scales (Mirzaei, 2021), a combination of both approaches to conduct urban climate modelling, based on the authors' knowledge, has rarely been pursued. Chang et al. (2018) integrated CFD and GIS to analyse urban ventilation corridors, using the capabilities from CFD to assess wind dynamics and from GIS for spatial analysis. Meng et al. (2017) combined CFD and GIS capabilities to integrate a GIS-based representation of an urban area into CFD to model microclimatic conditions. Both approaches used GIS solely for data input into CFD or to present CFD-generated datasets.
In this study we couple capabilities of CFD and GIS modelling, allowing to account for the complexity of the urban structure and specific surface characteristics on a fine spatial scale using GIS and to include the complexity of flow dynamics within the built environment using CFD. We integrate CFD-generated urban wind and air temperature patterns on a 3D basis into a refined 2D GIS-based micro-and bioclimatic model, introducing a 2.5D approach that combines CFD and GIS capabilities. We use the GIS-based modelling approach introduced by Back et al. (2021) to carry out finescale mapping of Land Surface Temperature (LST), Mean Radiant Temperature (MRT) and Universal Thermal Climate Index (UTCI) in a 2D urban environment. The UTCI is expressed as equivalent air temperature for a given combination of meteorological parameters  and is capable of representing bioclimatic conditions based on an assessment scale of ten categories from extreme cold stress to extreme heat stress (Bröde et al., 2012;Błazejczyk et al., 2013). In essence, we enhance the capabilities of this GIS-based approach by introducing Normalised Difference Vegetation Index (NDVI)-based methods to improve the Bowen ratio (β), emissivity (ε) and substrate heat flux (G) calculations, all of which contribute to deriving LST and UTCI. The approach based on Back et al. (2021) demonstrated its potential to rapidly assess thermal comfort for different urban forms and to identify priority areas for climate change adaptation measures. Although multiple time stamps during the day were analysed based on changing wind speed and air temperature, their variability within the built environment was neglected. Results for diurnal LST and UTCI development can thus only be assigned to a specific urban form, which limits the full potential of the fine-scale dataset to account for the development within the built environment.
To address this drawback, we resort to CFD enabling consideration of the horizontal and vertical variability in wind speed and air temperature patterns in the urban canopy layer. CFD-generated wind and air temperature datasets are integrated with the aforementioned GIS-based modelling approach to calculate LST and UTCI at a fine spatial scale and at different heights above the ground surface. This detailed data represents a significant improvement from substantially approximated velocity and temperature profiles previously used to calculate LST and UTCI. Combining the capabilities of CFD and GIS increases model complexity. Nevertheless, this enables enhanced fine-scale urban micro-and bioclimatic modelling. We explore the capabilities of this approach within a specific case study in the alpine city of Innsbruck, the capital of Tyrol in western Austria, analysing: (1) the effects of horizontal and vertical wind speed and air temperature variability on LST and UTCI, (2) the relationship between varying surface characteristics and wind speed, air temperature, LST and UTCI, and (3) the limitations and strengths of integrating CFD and GIS capabilities for fine-scale urban micro-and bioclimatic modelling.
Building upon Venter et al. (2021), who pointed out difficulties in relying on LST-based UHI assessments, our analysis confirms this statement and demonstrates that LST overestimates intraurban heat, as they do not properly account for the effects of changing wind speed and air temperature patterns in built environments.

Case study description
The case study is a specific area of interest in the city of Innsbruck, the capital of Tyrol in western Austria, situated at 47°16′ N and 11°24′ E in the Inn-valley at an altitude of approximately 574 m above sea level. Our focal area is located north-east of the city centre (black outline in Fig. 1). The degree of surface sealing subdivides the case study roughly into two areas: the northern area features a higher proportion of green areas, trees and single houses, whereas the southern area contains denser built forms, a lower proportion of green areas and less trees. This is illustrated using a land cover classification approach based on Hiscock et al. (2021) applied to the case study area. A Coloured Infrared (CIR)-image, shown in Fig. 1 for the entire city of Innsbruck, represents vegetation health status indicating the proportion of near-infrared light reflected by the vegetation surface. The image was used to calculate the NDVI, which forms the basis for calculating the Bowen ratio, emissivity (depicted in Fig. 1 for the area of interest) and substrate heat fluxes. Urban wind flow around Innsbruck is characterised by the channelling of the Inn valley from west to east, with the prevailing wind direction from west (Karl et al., 2020).

Data acquisition and software set up
We used a Digital Elevation Model (DEM) with an accuracy of 0.5 m, a CIR-image raster with an accuracy of 0.2 m (captured in late August 2016) as well as a building vector layer. DEM and CIR-imagery were provided by the local government "Land Tirol". The building vector layer of Tyrol (ESRI Shapefile format) was freely accessible (data.gv.at). We utilised the commercial software packages ESRI ArcMap v10.8.1 (ESRI, 2019) and Ansys® Fluent -Release 2020 R1, a finite volume method (FVM) solver, to conduct our analyses and evaluate simulation results. In Ansys®, we simulated wind speed and air temperatures, whereas in ArcMap, we calculated LST, MRT and UTCI based on the approach of Back et al. (2021). We measured essential meteorological parameters, such as air temperature (T a ), wind speed (U Wind ) and direction, relative humidity (RH), cloudiness (N) and vapour pressure (V p ), at a meteorological station within the case study area at a height of 6.0 m above ground (Fig. 2). LST measurements of various surface types were conducted on the 15th of July 2020 at around 15:00 using a thermal infrared camera and the associated software to postprocess the generated datasets (INFRATEC, 2015). These images were used for validation between modelled and measured LST. Fig. 2 shows the fields of view (1-4) for LST measurements on-site. Additionally, we approximated temperature around the building envelope by measuring LST of building walls from the four cardinal directions of selected buildings in the case study area using a thermal infrared camera (Fig. 3 b). The latter served as input to the CFD-based wind speed and air temperature simulations. Global (G R ), direct and diffuse radiation were calculated using the ArcMap's Area Solar Radiation tool including slope and aspect information from the DEM as well as radiation parameters according to atmospheric conditions on the selected time of day. Wind speed within the case study area was generally low during the measurement campaign, arriving from the west. As this is the prevailing wind direction in Innsbruck, we chose this wind direction and an initial wind speed of 4 m/s for the simulations using CFD to generate detailed air velocity-temperature contours approximating the conditions at the specified time of day.

Integrating CFD and GIS capabilities
Our approach follows several key steps (Fig. 3 A-G) representing a flow of data between GIS and CFD. To support modelling efforts, on-site measurements to obtain surface temperatures of the building envelope and a simplified 3D building model of the site were required. All steps are explained in the following sections in detail.

Step A: mean surface temperature calculation
We used the GIS-based modelling approach from Back et al. (2021) to calculate mean LST for specific surface types; trees, ground and buildings ( Fig. 3 A), relevant for CFD modelling, as well as to calculate the UTCI. The modelling approach combines a fine-scale surface classification, comprised of eight different surface classes, thermal characteristics (global radiation, direct radiation and diffuse radiation), surface characteristics (emissivity and Bowen ratio) and meteorological input data. Based on this data and well-established physical relationships in the model setup, the approach used an adaptation from Matzarakis et al. (2010) and Bröde et al. (2010) to first evaluate LST, followed by the MRT and finally the UTCI. LST calculation depends on atmospheric radiation, long-wave radiation flux density emitted by the surface and net all-wave radiation flux density, which were calculated (Eqs. (1)-(5)), based on Matzarakis et al. (2010) and previously described in Back et al. (2021): where T a is air temperature [°C] and B represents an approximation to the substrate heat flux [W/m 2 ]. The latter depends on net all-wave radiation flux density (Q) using the relationship: Q is a function of atmospheric radiation (A) and long-wave radiation flux density emitted by the surface (E) and is calculated using Eq. (3): where G R represents global radiation. Q incorporates E, which, in turn, incorporates LST (Eq. (4)). E is thus calculated using: where σ is the Stefan-Boltzmann constant (5.67•10 −8 [W/m 2 K 4 ]). To calculate Q and solve Eq. (4), A is calculated using Eq. (5): where V p is vapour pressure [hPa] and N is cloudiness [Okta]. LST and E are the two unknown variables, which are solved for using the system of Eqs. (1) and (4). The output is a map depicting LST, with the system of Eqs.
(1) to (5) solved for each cell across this map (Fig. 3 A). After solving this system and determining the mean values of different surface types, we use the following temperatures as input values to model wind velocity and air temperature distribution in the CFD simulations ( Fig. 3 E): ground

Fig. 2.
Dataset from the meteorological station within the case study area (T a -air temperature, V p -vapour pressure, U wind -wind speed, N -cloudiness and RH -relative humidity) and fields of view (1-4) for LST measurements using a thermal infrared camera.

2.3.2.
Step B: approximation of temperature data of the building envelope The GIS-based modelling approach is only capable of calculating surface temperatures visible from a bird's eye view, thus, excluding temperatures of the building walls. To address this drawback, we conducted on-site measurements of surface temperatures of specific buildings in the area, using a thermal infrared camera (Fig. 3 B). Mean surface temperatures of the four cardinal directions of the building envelope were used as input values in the CFD simulations to model wind velocity and air temperature distribution across the case study area (Fig. 3 E). Following temperatures were assigned: north (23.0°C), east (29.0°C), south (42.0°C) and west (27.0°C).

Step C: urban form approximation
To simulate the case study, we prepared a simplified 3D model of the urban infrastructure (Fig. 3 C). We have retained most of the infrastructure features including accurate dimensions and placement of the majority of the buildings (based on surveyed city maps), approximate number and height of trees using a digital elevation model of the region and size and extent of open areas. However, we do assume rooftops to be flat as opposed to the sloped rooftops (e.g., shed, gable or hipped roof) as generally seen in Innsbruck's architecture. This reduces the complexity of the model as well as simplifies and ensures that a consistent and manageable mesh can be generated for CFD simulations. Another simplification made is the removal of hedges and smaller vegetation as their minimal occurrence. Trees were also generalised as simple cuboids. The model of the area is padded with empty volume to allow for the flow to develop before interacting with the urban structure and to minimise the development of the non-physical reverse flow resulting in erroneous results.

Step D: mesh generation
Generating an adequate mesh entails dividing the complete volume into smaller, more simplified units. For each of these units, the Navier-Stokes equations are iteratively solved by the Ansys® FLUENT solver. By controlling the size of edge and length we can control the granularity (number of elements) and thereby the quality of the mesh, as shown in Table 1.
A vital aspect of meshing is to select an optimum degree of regiment such that the errors fall within an acceptable range while maintaining  low computational requirements. To this end, we conducted a grid convergence study based on the guidelines provided by Celik et al. (2008). According to the used meshes, the velocity at approximately the centre of the area of interest was chosen as a parameter for this study. Details of each of the tested meshes are summarised in Table 2. The average cell size (h), is computed based on the number of cells N used to discretise the space according to the formula in Eq. (6), where V is the volume of each cell. The grid refinement ratio (r 2,1 ) between the fine and medium grid sizes is calculated as the ratio between h1 and h2 representing the average cell size of the two meshes respectively: The order of convergence (p) of about 9.789 is computed using Eq. (8): where ε 32 = ϕ 3 − ϕ 2 and ε 21 = ϕ 2 − ϕ 1 with ϕ n being the solution of the n th mesh. The approximate relative error can thus be calculated as: The error was found to be 6.32 %. With this information, we can then calculate the Grid Convergence Index (GCI) using Eq. (10), This solidifies the numerical accuracy for the fine mesh as within acceptable limits for a CFD simulation.

Step E: CFD simulations
With Ansys® FLUENT, we generate a steady-state simulation of the case study area and created a detailed air velocity-temperature contour for our area. The contours generated are dependent on the urban infrastructure along with tree canopy within the vicinity, and the temperature of various surfaces (i.e., ground, trees, rooftops, and walls in different directions), previously described in Sections 2.3.1 and 2.3.2. Such a fluid-structure simulation allows us to incorporate the vertical and horizontal variability in wind speed and air temperature patterns in the calculations of the LST and UTCI, described next.

Step F: translation of CFD data into GIS
We transferred CFD data into GIS using a six-step approach in ArcMap 10.8.1, comprising the tools: (1) point to raster, (2) geo-reference, (3) raster to point, (4) clip, (5) natural neighbour and (6) extract by mask. Wind velocity and air temperature datasets were converted from the CFD software into text format (.csv file) with assigned x and y coordinates. In ArcMap 8.1.3 we added the x and y datasets with attached wind speed and air temperature values and converted our point dataset into a spatial raster. This enabled georeferencing accordingly to the coordinates of the case study area. Subsequently, the data were reconverted into a point dataset to enable interpolation using the natural neighbour technique. This operation created a raster dataset, from which buildings were extracted using a defined mask.

Step G: refined micro-and bioclimatic modelling
To enhance speed, accuracy and capabilities of the GIS-based model, we included NDVI-based methods to improve the Bowen ratio (β), emissivity (ε) and substrate heat flux (G) calculations, all of which were used to calculate LST and UTCI. Using the NDVI to calculate these parameters enhances capabilities to more profoundly account for variability in surface characteristics. Additionally, the surface classification used to assign emissivity and Bowen ratios in the original model is no longer necessary, leading to a reduction of the data sets required to run the model. Finally, we introduced a sky view factor (SVF)-based approach to refine the calculations for downwards orientated long-wave radiation in urban areas within the GIS-based modelling approach. In addition to the atmospheric radiation, we included long-wave radiation emitted by the urban structure and, thus, approximated total downwards orientated long-wave radiation in urban areas. We used the CFD-generated wind and air temperature datasets translated into the enhanced GIS-based modelling approach to calculate LST and UTCI at a fine spatial scale and at different heights above the ground surface. To calculate LST, we used wind velocity and air temperature datasets at the height of 0.2 m. To calculate MRT and UTCI we used the datasets at the heights of 0.2 m and 1.75 m. How we calculate these metrics in the GISbased approach is explained in subsequent paragraphs.  Matzarakis et al. (2010) to calculate downwards orientated long-wave atmospheric radiation (A). This represents a good approximation for A, when studying large areas on a regional scale or flat rural areas on a regional to microscale. However, when the focus is set on specific urban areas, urban structure needs to be considered, as A changes with the amount of visible sky for every considerable location on the surface indicated by the Sky View Factor (SVF). The SVF describes the proportion of visible sky above one particular observation point. Accurate determination of the SVF is crucial to account for the effects of the built structure (e.g., building height and geometry) on net radiation in urban areas and is therefore essential in describing urban climatology at the microscale (Dirksen et al., 2019). Additionally, long-wave radiation is emitted by the urban structure and contributes to total downwards orientated long-wave radiation in urban areas (R L ). Other components such as the reflectance of long-wave radiation at the façades are neglected here. Therefore, R L is approximated using Eq. (11),

Downwards orientated long-wave radiation
The proportion of long-wave radiation emitted by the urban structure, E•(1-SVF), is divided equally representing a simplification of a vertical façade. Furthermore, for simplicity, we assume that surface temperature of the façades are equal to that of the ground. Therefore, E in Eqs. (11) and (12) is defined as long-wave radiation flux density emitted by the surface and is equal for both ground surfaces and façades.
In contrast to Eq. (4), here in Eq. (12) we replaced A with R L in Eq. (12) to calculate E. Long-wave radiative fluxes from the surfaces and the atmosphere are weighted based on the SVF, which is calculated using the freely available software System for Automated Geoscientific Analysis -SAGA (Conrad et al., 2015).

Fine-scale Bowen ratio calculation
Rigo (2006) introduced a second-degree polynomial regression function to calculate the Bowen ratio from NDVI band values based on in-situ data conducted during the Basel Urban Boundary Layer Experiment (BUB-BLE) in 2002. However, this function does not cover the entire NDVI range, only accounting for NDVI band values between −0.2 and 0.4. To account for a wider range of NDVI band values, including nonvegetated surfaces and therefore consider the intraurban variability, we established an exponential regression between the Bowen ratio (β) and the NDVI based on the existing second-degree polynomial regression function after Rigo (2006) and on literature-based boundaries, denoted by Eq. (13): To obtain Bowen ratios for NDVI band values between −0.5 and −0.2 and between 0.4 and 1.0, we fitted a curve to the existing second-degree polynomial regression function based on scientific literature-based boundaries for Bowen ratio values in urban areas (Grimmond and Oke, 2002;Christen and Vogt, 2004;Oke et al., 2017) (Fig. 4). NDVI band values closer to −1 indicate water body surfaces, which are not addressed in this study. Bowen ratios range from 0.1 to 8 (Fig. 4), which agree with those occurring in urban areas (including cropland) (Grimmond and Oke, 2002, Christen and Vogt, 2004, Oke et al., 2017. Within NDVI band values of −0.2 and 0.4, the second-degree polynomial regression function from Rigo (2006) (blue dots in Fig. 4) coincides well with our exponential function (R 2 = 0.998).

Substrate heat fluxes
Substrate heat fluxes (G) (Eq. (14)) was calculated based on a relationship with NDVI for urban areas (Kustas and Daughtry, 1990;Parlow, 2003;Rigo and Parlow, 2007). This approach assumes that substrate heat fluxes decrease with an increase in biomass, as indicated by the NDVI (Rigo and Parlow, 2007).
where Q represents net all-wave radiation and is calculated using Eq. (15).
which differs slightly to Eq. (3) in that we substitute A with R L in Eq. (15) to calculate Q, now including a SVF-based approach to calculate A and, additionally, the proportion of long-wave radiation emitted by the urban structure.

Surface emissivity
The surface emissivity (ε) is essential in urban surface temperature calculations (Mitraka et al., 2012;Kikon et al., 2016). To improve emissivity values for different surface characteristics we used a combination of two approaches: (1) we associate emissivity values from scientific literature (Coutts and Harris, 2012;Kotthaus et al., 2014;Song and Park, 2014;Mandanici et al., 2016) to specific surface classes based on a modified classification approach from Hiscock et al. (2021) for NDVI band values between −0.50 and 0.135 and (2) we used an approach based on Van De Griend and Owe (1993), also described by Zhang et al. (2006) and Liu and Zhang (2011) Conversely to Zhang et al. (2006) and Liu and Zhang (2011) who proposed a range between 0.157 and 0.65, we found vegetation signals between NDVI band values 0.135 and 0.157 in this specific case study. Based on the Bowen ratio we found latent heat to exceed sensible heat (Bowen ratio below 1) with NDVI band values above 0.135. We used this threshold to distinguish between vegetated and nonvegetated surfaces. Eq. (16) therefore calculates the emissivity for NDVI band values between 0.135 and 0.65. Table 3 represents emissivity values according to the NDVI band values −0.5 to 0.65.

Adapted LST and UTCI calculations
After adapting the system of equations including enhanced input values for the Bowen ratio, emissivity and substrate heat fluxes, we integrated the CFD-generated wind and air temperature datasets into the enhanced GISbased modelling approach and calculated LST and UTCI at a fine spatial scale using the datasets at the heights of 0.2 m and 1.75 m. We adapted LST and UTCI calculations by providing enhanced input values as described in the previous sections. LST was calculated using Eq. (17): We replaced the approximation of substrate heat fluxes (B) (used in Eq. (1)) with the NDVI-based method to calculate substrate heat fluxes (G). Enhanced emissivity and downwards orientated long-wave radiation values come to bear in calculating E (Eq. (12)), as well as Q (Eq. (15)). Again, both values E and LST are the two unknowns evaluated iteratively by using the system of Eqs. (12) and (17). MRT was calculated based on Back et al. (2021) using Eq. (18): where I ⁎ is the radiation intensity of the sun on a surface perpendicular to the incident radiation direction, a k is the absorption coefficient of the Fig. 4. The introduced exponential regression function (black dashed line) based on values using a second-degree polynomial regression function (Rigo, 2006) for NDVI band values between −0.2 and 0.4 (blue dots) and the corresponding equation to calculate the Bowen ratio from NDVI values. The vertical axis depicts the Bowen ratios, and the horizontal axis depicts the NDVI range used in this study. irradiated body surface area of short-wave radiation (standard value for the human body of 0.7), ɛ p is the emission coefficient of the human body (standard value of 0.97) and f P is the surface projection factor. Finally, we calculated UTCI using Eq. (19): ( 1 9 ) where P Vapour is the water vapour pressure [kPa].

Effects of wind and air temperature variability on LST and UTCI
The horizontal and vertical variability in wind speed and air temperature has substantial effects on LST and UTCI. Wind channelling is visible in the open spaces and street canyons with west-east orientation, leading to an increase in wind speed and consequently to a decrease in air temperatures (Fig. 5). For street canyons with north-south orientation, both variables show an inverse relationship. Especially courtyards and leeward areas where the wind is blocked and deflected by buildings, wind speed decreases and leads to an increase in air temperatures (Fig. 5). Such areas show higher LST and UTCI values and can be considered as thermal hot spots (Fig. 7). Similar to buildings, trees lead to a deceleration in wind speed, however, decrease air temperatures in leeward areas. This effect intensifies with increasing tree cover, leading to a decrease in UTCI values (effect is less pronounced for LST values), and indicating thermal cool spots with enhanced human thermal comfort (Fig. 7). With an increase in height, wind speed increases, while air temperatures decrease (Fig. 6). This is illustrated in Fig. 6 for a location in an open street canyon. This specific location in the case study area is characterised by an NDVI of −0.08, a Bowen-ratio of 1.98, a SVF of 0.68 and a surface temperature of 47.01°C. The decrease in wind speed near ground due to an increase in the surface roughness is highlighted. However, visible in our analyses, built environments have more complex wind and air temperature patterns especially in the lowest 2 m above the ground surface of the urban canopy layer (Fig. 5). Whilst values at a height of 1.75 m are particularly relevant for the human thermal comfort, values at lower heights, e.g., 0.2 m in this case study, can be more relevant for animals. Differing in wind and air temperature patterns, both heights show slightly different areas of maximum UTCI values (Fig. 7), with UTCI values at lower height more influenced by the surface temperatures due to lower wind speed.

The role of varying surface characteristics
The diversity of surface characteristics is expressed using NDVI band values ranging from −0.5 (sealed surfaces) to 0.65 (vegetated surfaces). LST is highest with lowest wind speed and lower NDVI values (Fig. 8). Decreasing NDVI band values lead to increasing differences between LST and air temperatures. Additionally, with decreasing wind speed LST increases at a lower rate for NDVI band values 0.3 and 0.5 (vegetated surfaces) compared to NDVI band values −0.5 to 0.1 (sealed surfaces). This effect is visible in Fig. 8 in the contrast between trends of air temperature and LST respectively. Comparing all four variables shows a clear dependence of air temperature to wind speed, the strong interaction of both variables with LST and their variability with different surface characteristics.
We observe similar interactions and dependencies between wind speed, air temperature, NDVI and UTCI as we did for LST, however less pronounced (Fig. 9). The variations in NDVI band values strongly correlate with LST and UTCI. High UTCI values indicate hazardous locations of extreme heat, posing potential health risks. Such conditions can develop in areas with an increased degree of surface sealing (lower NDVI band values), lower wind speeds and higher air temperatures. In contrast, lower UTCI values can be found in shaded and densely vegetated areas (higher NDVI band values), indicating thermal cool spots.

LST-based studies overestimate thermal discomfort
Our results highlight the effects of varying wind speed and air temperatures at different heights in built environments and, thus, coincide with the findings from Venter et al. (2021). They found LST-based UHI to overestimate heat stress and the contribution of urbanization to the local temperature, when compared to crowdsourced air temperatures. The complexity of built environments leads to both increasing and decreasing wind speeds vertically as well as horizontally and on small time scales. Consequently, this variability leads to changes in air temperatures and subsequently in LST and UTCI values. Our results show that the human thermal comfort (UTCI) can vary substantially within the first 2 m above the ground surface in the urban canopy layer. UTCI values at elevated heights (1.75 m) are more influenced by changing wind speed and consequently air temperatures. With wind speed generally lower close to the ground surface, UTCI values at lower heights (0.2 m) are more influenced by rising air temperatures which in turn are stronger interacting with LST. As we have demonstrated the strong correlation between surface characteristics, wind speed, air temperature, UTCI and LST (Figs. 8 and 9), we conclude that LST represents the thermal conditions near ground level rather than those elevated from the ground (1.75 m), which are, however, more representative of the human thermal comfort. We confirm that studies solely based on LST overestimate heat stress and, thus, neither contribute to a holistic view of Fig. 6. Wind speed, air temperature and Universal Thermal Climate Index (UTCI) with changing height above ground at one location in an open street canyon in the case study area. The vertical axis depicts the height above the ground surface, and the horizontal axes depict wind speed, air temperature and the UTCI respectively. urban heat development nor towards comprehensive heat assessments and hot spot analysis.

Limitations and strengths of the integrated approach
Due to the disadvantage in computational time in CFD modelling, we calculated micro-and bioclimatic conditions only for one specific time step of the day considering the associated meteorological conditions. Future work should include modelling over a diurnal cycle to include the effects of changing surface and building wall temperatures on the air velocity-temperature contour and the effects of changing initial wind speeds and air temperatures on their distribution across the case study area. We compared modelled and measured LST and found a good agreement between the measured datasets and the modelled datasets using the integrated CFD-GIS modelling approach (see Section 3.4.2). To validate other variables, such as wind speed or air temperature distribution across Fig. 8. Correlation between wind speed, air temperature, Land Surface Temperature (LST) and varying surface characteristics represented by specific NDVI band values (0.5 to −0.5). The vertical axis depicts the temperatures (LST and air temperature), and the horizontal axis depicts the wind speed. Fig. 9. Correlation between wind speed, air temperature, Universal Thermal Climate Index (UTCI) and varying surface characteristics represented by specific NDVI band values (0.5 to −0.5). The vertical axis depicts the temperatures (UTCI and air temperature), and the horizontal axis depicts the wind speed. different heights, would, however, require a very extensive measurement campaign in order to compare measured variables with modelled variables at such a high resolution predefined by the CFD simulations. Furthermore, given that most variables are used in the model setup to calculate LST, we believe that this comparison is a good proxy for the other variables. Nevertheless, with this study, we demonstrated the added benefit of combining CFD and GIS and explored the implications of the variability in wind speed and air temperature on urban heat and thermal comfort assessments.

Comparing CFD-integrated and non-integrated modelling
CFD-integrated (Fig. 10 b) and non-integrated (Fig. 10 a) modelling shows stark differences in magnitude and spatial variability of LST. Nonintegrated CFD modelling implies that wind speed and air temperature are set to constant 4.0 m/s and 23.7°C respectively for every single raster cell defined in our study area. CFD-integrated modelling, however, accounts for the additional horizontal and vertical variability in both wind speed and air temperature, leading to a wider range of and, perhaps, more realistic LST patterns (Fig. 10 b). With an initial moderate air temperature (23.7°C) and wind speed (4.0 m/s) the distribution of both variables throughout the case study leads to substantial variations, resulting in air temperatures and LST as high as 38.5°C and 60.0°C respectively. Nonintegrated CFD modelling prevents LST temperatures from exceeding 37°C, underestimating measured surface temperatures. This is clearly visible when comparing modelled and measured datasets (Fig. 12), with CFDintegrated LST modelling showing good agreement with measured LST from thermographic surveys.
Furthermore, non-integrated CFD modelling (Fig. 10 d) consequently neglects the effects of horizontal and vertical wind speed and air temperature variations on the human thermal comfort (UTCI). With wind speed and air temperatures unchanged for every raster cell (Fig. 10 c), maximum UTCI only reaches a value of 29.6°C. The homogeneous distribution of UTCI values across the different surface characteristics shows that too much emphasis has been placed on surface temperatures alone, neglecting how varying air temperatures, wind speed and surface characteristics influences thermal comfort. Using the CFD-integrated approach, these effects become more pronounced (Fig. 10 b, d), thereby highlighting the importance of horizontal and vertical wind speed and air temperature variability for studies on human thermal comfort.

Comparing modelled and measured LST
Comparing modelled with measured LST shows overall good agreement in terms of temperature range, extreme values and distribution (Fig. 11). Unfortunately, as calculations of modelled LST are based on a dataset from 2016 and measurements were conducted in 2020, vegetation health status is not directly comparable. Notwithstanding, no visible changes can be observed in the placement of vegetated areas and trees, buildings and mobility infrastructure. We see vegetated areas appearing cooler than sealed surfaces in the modelled dataset, however not with the same magnitude observed during the measurement campaign (Fig. 11). Comparing measured LST of dry grass with modelled LST of vegetated surfaces shows similar temperatures (Fig. 11). Similar to measured LST on very bright surfaces (road markings), we see lower LST values in the modelled dataset, however, again a mismatch in magnitudes (Fig. 11 Detailed Section from field of view 1). This is explained by the fact that the lowest surface emissivity in our modelling approach has a value of 0.85. To represent the effect of white paint on the roads, lower emissivity values will need to be specified for these surface types. Modelled LSTs for sealed surfaces, especially roads, are in good agreement with measured data showing temperatures around 50°C.

Improving downward orientated long-wave radiation
We calculated total downward orientated long-wave radiation in urban areas as explained in Section 2.3.7. Comparing spatial distribution and magnitudes of this variable with downwards orientated long-wave atmospheric radiation defined by Eq. (5) shows significant differences (Fig. 12). Neglecting the SVF (Fig. 12 a) leads to higher emphasis on the air temperature within the initial equation calculating atmospheric radiation (see Eq. (5)). This results in a spatially distributed atmospheric radiation that is similar to air temperature patterns across the study area of interest (Fig. 12 a). Including the SVF more strongly emphasises the urban structure and influence of building and vegetation heights on total downwards orientated long-wave radiation (Fig. 12 b). Air temperature still affects atmospheric radiation (with Eq. (5) incorporated in Eq. (12)) but these patterns are less pronounced (Fig. 12 b). When comparing both approaches to calculate LST and UTCI, we observe differences in both spatial distribution and magnitude. LST mean values decrease from 36.8°C to 35.1°C, while UTCI mean values decrease from 31.0°C to 30.5°C across the study area of interest when using the SVF-based approach. Both LST and UTCI values decrease with decreasing SVF, leading to larger differences in LST and UTCI values within tree alleys and street canyons between both approaches. These differences are less pronounced in open spaces. Additionally, using the SVF-based approach refines the determination of net all-wave radiation, which, in turn, is fundamental for further calculations.

Conclusions
We have introduced a 2.5D approach coupling the capabilities of CFD and GIS modelling for fine-scale urban environmental studies. The approach enables comprehensive analysis of the effects of wind speed and air temperature variability on the surface temperatures (LST) and the human thermal comfort (UTCI -Universal Thermal Climate Index) in the urban canopy layer. Using the approach in a specific case study area in Innsbruck, Austria, we demonstrated the relationships between wind speed, air temperature, surface characteristics (based on the NDVI), LST and UTCI and highlighted the importance of the variability in wind and air temperature patterns and the consideration of different heights to locate thermal hot spots in urban environments.
Results show that UTCI values can vary substantially within the first 2 m above the ground surface in the urban canopy layer. With wind speed generally lower near the ground surface, the strength of the correlation between air temperature, UTCI and LST decreases with increasing height. As such, LST correlates with thermal conditions near the ground surface rather than with those at elevated heights (1.75 m), which are, however, more representative of the human thermal comfort. Our evidence supports the hypothesis that studies solely based on LST overestimate human thermal discomfort and, thus, neither contribute to a holistic view of urban heat development nor towards comprehensive heat assessments and hot Fig. 11. Measured (1-4) and CFD-integrated modelled LST for the specific location within the case study area of interest. spot analysis. Considering these findings, we urge caution when relying on LST-based studies alone for urban heat studies, as this could lead to misinterpretations of hot and cool spots in cities and thus provide misleading insights into the anticipated effects of adaptation measures.
In this study we have successfully demonstrated the capabilities of combining CFD and GIS and introduced an integrated approach capable to improve intraurban studies towards comprehensive assessments of urban heat at a very fine scale. The results of this study and further applications of the integrated approach can provide a better understanding of the interactions between wind velocity, air temperature distribution, varying surface characteristics and micro-and bioclimatic conditions in the urban canopy layer and can contribute to sustainable urban planning fostering resilient cities.

Data availability
Data will be made available on request.

Declaration of competing interest
The authors declare that they have no competing interests.