Modeling the intra-urban nocturnal summertime air temperature fields at a daily basis in a city with complex topography

Detailed knowledge about the intra-urban air temperature variability within a city is crucial for the implementation of adaptation strategies to counteract the negative effects of urban heat stress. Various methods to model urban-rural temperature differences exist, but they often only cover certain periods (heatwave, hot day) or meteorological conditions (sunny and calm) due to computational limitations or limited data availability. Here, we present a land use regression approach to model nocturnal air temperature fields for every single night of the summers 2018 to 2020 in a city with complex terrain (Bern, Switzerland). Furthermore, we investigate the applicability of different model structures and straight-forward computable GIS variables to model cold air drainage, which exerts an important influence on the local-scale climate of cities with complex terrain. The geostatistical models are calibrated with in-situ data of a dense low cost air temperature measurement network and high resolution spatiotemporal (land use and meteorology) data, which are all publicly available. The resulting land use regression models are capable to model and map intra-urban air temperature differences with a good model performance (R 2 : 0.65–0.71; RMSE: 0.69–0.76 K). Evaluations with data from additional measurement stations and periods (summer 2021) show that the models are able to estimate different meteorological and spatial conditions, but that the representation of small-scale topographic features remains difficult. However, the comparatively low computational and financial effort needed to calculate nocturnal air temperature fields at daily basis enable new applications for cities with restricted resources for various areas of interest, such as urban planning (e.g. effect of heat mitigation policies) or heat risk management (e.g. analyze small-scale


Introduction
Mean summer air temperatures as well as the likelihood for air temperature extremes such as very hot days (T max � 35˚C) and tropical nights (T min � 20˚C) are increasing across Europe due to the ongoing anthropogenic warming [1,2]. Cities are particularly prone to this shift in summer air temperatures due to the so-called urban heat island (UHI) effect which leads to excess temperatures in urban environments [3], especially during night [4]. The UHI effect is mainly caused by high heat capacities of the used construction material, large fractions of impervious surfaces, low albedos due to the rather dark surfaces and a lack of vegetation and ventilation [3]. Since various negative impacts for cities regarding public health (e.g. increased fatigue, cardiovascular problems), environment (e.g. increased vulnerability of urban trees) or economy (e.g. cooling costs, reduced productivity of people) may arise from urban heat excesses, adaptation strategies to counteract UHIs should be implemented [5].
To plan and implement adequate adaptation measures, detailed knowledge about the air temperature distribution within a city is crucial. Since observational networks in urban environments imply relatively high costs [6] and national weather services typically do not maintain measurement stations within the city [7] different methods to assess and map urban temperatures have been applied in the recent past. The most widespread approaches include satellite-based land surface temperature analyses (which is not comparable with air temperature since the UHI effect is the largest during the day; e.g. [8,9]), numerical urban climate model simulations (e.g. [10,11]), or analyses based on self-designed urban air temperature measurement networks (e.g. [12][13][14]). Often, only limited timespans with homogeneous synoptic conditions favoring large UHI intensities (e.g. hot days or heatwaves) are investigated due computational limitations (e.g. process based numerical urban climate models) or limited data availability (e.g. satellite images with clouds). However, for several applications such as the assessment of heat risk for the urban populations (e.g. [15]) or the planning of UHI mitigation strategies (e.g. [16]), timespans also covering a variety of synoptic conditions are needed. Additionally, extended modeling periods would also enable the possibility to calculate and investigate threshold statistics such as the number of tropical nights in a specific neighborhood.
The combination of in-situ data from a low cost urban air temperature measurement network with publicly available spatial data offers the opportunity of conducting daily models and maps of entire summer periods with geostatistic land use regression (LUR) models with relatively low computational and financial effort [17]. Since the formation and the magnitude of the UHI effect depends to a large extent on meteorological preconditions [18], it is crucial to incorporate meteorological variables, if a larger timespan is analyzed [4,17]. Many studies have been conducted using meteorological variables as predictors of UHI intensities, of which wind speed and cloud cover (solar radiation) were often identified as the meteorological variables having the greatest impact on the UHI [4,[18][19][20][21]. In addition, other variables such as relative humidity [20,21] or daily temperature range [4] were identified as predicting variables. Precipitation is rarely investigated, since many studies focus on weather conditions that exacerbate urban heat, which implies that days with precipitation are often excluded from the analyses [4,14], although the reducing effect of precipitation on the UHI intensity is evident [22].
In addition to meteorological and land use patterns, cold air drainage (CAD) might be another important factor shaping the air temperature distribution within a city [23,24]. Such cooling airflow systems at local scales are mainly caused by drainage due to hilly terrain [25] and depend on the interaction of such topographical features and urban morphology. In existing LUR studies, such specific cold air variables are often not incorporated (e.g. [4,26,27]) or the focus is set on ventilation paths (e.g. with the proxy front area index, [28,29]) which do not account for CAD caused by terrain features, since the spatial connectivity between cold-air source areas (production) and urban areas (impact) is not considered [24]. However, approaches to model topographic proxies for CAD exist and can be grouped in two categories: First, flow accumulation proxies that take into account that cold air is denser than warm air and that thus the flow and accumulation of cool air parcels can be modeled similar to water [17,30]. Second, relative elevation proxies that compare the position (elevation) of a grid cell with the surrounding grid cells and expect cool air to accumulate in concave landforms [31][32][33][34]. However, some of those studies were conducted in the field of ecology focusing on environmental conditions outside of urban environments [30][31][32][33]. The applicability of those proxies in urban environments has thus received limited attention and comparisons between the different approaches are lacking.
The aim of this study is to model nocturnal UHI intensities and map summertime nighttime air temperatures of Bern, Switzerland for every night during the summers 2018 to 2020, using a LUR approach that combines publicly available land use and meteorological data. Based on data of a dense urban measurement network [35] and previous LUR models [17], we expand the models from heatwave means to daily fields and evaluate different model structures and CAD variables to estimate the nocturnal air temperature variability in a city with complex topography.

Study site
With 134'000 inhabitants, Bern is a medium-sized city in Switzerland [36]. The city is located in the mid-Western part of the Swiss Plateau at a mean elevation of 550 m.a.s.l. and is characterized by its complex topography: Several hills, including two which reach more than 850 m. a.s.l. in the northeast (Bantiger, 947 m.a.s.l.) and in the south (Gurten, 858 m.a.s.l.), various depressions, and a river (Aare) shape the natural topography of Bern and its surroundings [23] (Fig 1). Important regional effects of the complex topography are air temperature inversions with different thickness layers that evolve in summer during every second night [37] and the canalization of the wind due to the blocking of the nearby mountain ranges Jura and Alps which leads to a domination of southwestern and northeastern wind systems [23] (Fig 1A). During calm nights, also weak katabatic winds from southeast occur [23]. Bern experiences warm humid summers with a mean summer air temperature of 17.4˚C and a mean summer precipitation amount of 333 mm throughout the climatological reference period 1980-2010 [38].

Data
Urban air temperature data. From 2018 to 2021, a dense urban air temperature measurement network is operated in the city of Bern throughout summertime. It consists of 70 to 90 stations, which are installed from mid-May to mid-September and record the air temperature at an interval of 10 min. The stations are mounted at a height of approximately 3 m at freestanding poles and cover the microclimatic and topographic heterogeneity of the city. Even though the stations consists of low cost devices (about 65 CHF per station) and are only passively ventilated, they are able to capture nighttime temperatures (22 pm to 6 am) with only small biases when compared to professional actively ventilated sensors (mean bias: -0.12 to 0.23 K, [35]). During daytime, a mean positive bias of 0.61 to 0.93 K is observed, which implies that the uncertainties due to the measurement devices might be larger than the observed daytime urban heat island [17,35]. Due to these constraints, we focus on the nocturnal temperature data in this study. Hence, for the calibration of the models, we use mean nighttime temperatures (22 pm to 6 am) of 61 stations that were placed at the same location from 2018 until 2020 (276 days; Fig 1B). Data from 17 additional stations in 2018 as well as data of the remaining 55 calibration station in summer 2021 is used for the evaluation (92 days each).
Meteorological reference data. The official measurement station in Bern maintained by the Federal Office of Meteorology and Climatology (Meteoswiss) is located about 5 km north

PLOS CLIMATE
of the city center in Zollikofen at an elevation of 553 m.a.s.l. (Fig 1B). It is a WMO-certified station, which should per definition not be influenced by urban heat [39]. Out of the 12 different meteorological variables recorded, we use air temperature at two meters as the basis temperature for the models as well as relative humidity, precipitation, global solar radiation and wind speed and wind direction as explanatory variables (Table 1), which are averaged over a 24 h (6 am to 6 am) and a nighttime 8 h interval (22 pm to 6 am). The data is publicly available and can be downloaded on the website of Meteoswiss [40].
Land use variables. In order to model the influence of the land use on the urban nighttime air temperatures, 14 land use and 4 additional CAD variables were calculated (Table 1), which are shortly described in the following. For more detailed information about the computation of the land use variables, see [17].
Six land cover parameters were derived from the cantonal cadastral survey [41] Cold air drainage variables. In total four different CAD variables of two different groups were derived. The first group was calculated with the commands "flow direction" and "flow accumulation" in ArcGis (version 10.7.1) [46] using a high-precision digital elevation model (DEM) of Bern [42]. With these features, the flow from a cell of a raster to its neighboring cells is calculated (according to the height of the cell) and then accumulated over the entire area. This approach was originally developed for hydrological analyses [47]. Applied to CAD, it follows the idea that cold air behaves similar than water: cool dense air follows the topography [33] and accumulates in the steepest descents [30]. This method was already tested for Bern in an earlier study. Although it worked well in the model evaluation, the amplitude of cooling appeared to be overestimated at some non-monitored locations when mapping it [17]. In this study, different DEM resolutions (25 to 250 m) were tested and logarithmic variables were calculated. One of the variables is calculated using a DEM (Flow Accumulation; FA) and one with buildings and vegetation placed on the DEM (Flow Accumulation with buildings; FAB) ( Table 1). This method implies that the buildings block CAD, but has the disadvantage that CAD may also form at the rooftops of buildings.
The second group of CAD variables was calculated based on comparisons of the elevation of a grid cell with its neighboring cells. The first approach compares the elevation of a site with the minimum elevation within a specified radius in order to distinguish locations where air would drain away (high value) or cold air could pool (low value). This so called "Relative Height" REH is the third CAD variable used here and has been tested in earlier studies [31,32]. If not the minimum, but the mean elevation of the surrounding grid cells is subtracted from a grid cell, the "Topographic Position Index" (TPI) is calculated [34]. Negative TPI values indicate concave positions (valley) where cold air drainage is favored, positive TPI values indicate convex positions (hilltop) where no cold air drainage is expected. TPI is the forth CAD variable tested in this study.
Buffer zone radii. In LUR, temperatures at a specific location are expressed as functions of land use characteristics in circular areas surrounding the location (buffers). Depending on the land use variable, the radius of that area is chosen differently [48]. Here, we tested buffer radii values from 25 to 1000 m. The selection of the best fitting buffer zone radii per land use variable was conducted with linear regression modeling (see [17] for detailed description, Table 1).
Model structure. The study aims to develop spatiotemporal geostatistical models for estimating the nighttime (22 pm to 6 am) air temperature fields in Bern with the antecedent mentioned data. We use the mean air temperature difference of a station of the urban network to the reference station in Zollikofen as the response variable of the model. Although the definition of stations to be "rural" and "urban" is problematic and hence the use of term "UHI intensity" not fully appropriate [49], we use that term for these air temperature differences due to simplicity reasons. To calculate absolute nighttime air temperature fields, we add the observed mean nighttime air temperature of Zollikofen to the modelled UHI intensities.
As first step of the modeling process, we analyze the independent influence of meteorology and land use on the UHI intensity to decide upon which variables are used for the further modeling. Whereas the investigation of the meteorological variables is important for the estimation of potential UHI intensity in general (e.g. strong UHI intensity expected after a sunny day and a calm, dry night), the land use variables are used to estimate the urban heterogeneity and intra-urban distribution of temperature. As second step, three different model structures to combine meteorological and land use variables are tested: First, we build an additive model (ADD ; Table 2), where meteorological and land use variables are combined with a multiple linear regression (MLR) model similarly to the approach of [28] applied to Honkong. Second, a two-step multiplicative structure is tested (MULT, Table 2) according to the diagnostic UHI equation model for cities in northwestern Europe [4]. Here, a meteorological MLR model is first fitted, independently from the land use variables. Based on that meteorological model, a best fit of constant is calculated for every station, based on the mean absolute error in order to get the equation "constant × meteorological model". As a last step, another MLR is conducted to represent the constant with land use variables that will have been selected in the previous chapter [4]. As a third model structure, we test a novel approach of an interactive model structure (INT, Table 2). The idea of that approach is that possible interactions between meteorology and land use variables are considered [50]. Non-significant interaction pairs (p > 0.001) are omitted from the model. Although the computation of the model is rather simple, the structure gets more complicated (Table 2).

Mapping
To illustrate the outputs of the different LUR models, we map the modeled air temperatures with a 50 m x 50 m resolution for the greater area of Bern using R [51]. We show the predicted air temperature distribution of two meteorologically different nights with conditions favoring and hindering strong UHI intensities.

Evaluation
To evaluate the models, additional station data is used. We use data from 17 additional stations that were only installed during summer 2018 (Fig 1B) as well as data of the remaining calibration stations that were installed during summer 2021 (55 out of 61). Explained variance (R 2 ), root mean square error (RMSE) and mean bias (MB) are calculated in order to assess the performance of the models. Additionally, we compare the ability of the different models to predict T Sij represents the observed air temperature at station j of the urban temperature network for the night i. T Z represents the observed nighttime air temperature at the reference station in Zollikofen. M denotes the temporal and L the spatial variables of whom n and m exist. a t represents the slopes of the temporal, b s the slopes of the spatial and c q the slopes of the combined predictors. https://doi.org/10.1371/journal.pclm.0000089.t002

PLOS CLIMATE
nighttime air temperatures over complex topographical features with additional data from three locations at the Gurten hill south of the city at elevations of 848, 694 and 631 m.a.s.l. mounted in summer 2021, as well as with data from three stations which were mounted during summer 2019 in a small-scale topographic depression within the built up area ("Egelsee") in the east of the city (Fig 1B).

Independent selection of meteorological and land use variables
Meteorology. For this analysis, the mean nighttime UHI intensity per day of all 61 stations is calculated (n = 276) and compared with the corresponding value of the meteorological variables (Fig 2). A positive correlation of global solar radiation (r = 0.58; Fig 2A) and UHI intensity, as well as negative relations between precipitation and UHI intensity (r = -0.47; Fig  2B) and relative humidity and UHI intensity (r = -0.49; Fig 2D) can be observed. The influence of wind is less evident, a slight increase in UHI intensity with stronger wind is found (r = 0.06; Fig 2C).
Further analyses on the wind showed the important role of the wind direction. Enhanced winds from north (average between 270 and 90˚) is associated with an increase of the UHI intensity (r = 0.41), while winds from south result in a decrease (r = -0.63). This feature is the most prominent if nighttime winds are analyzed (Fig 2E). Subsequent analysis on precipitation showed that the occurrence of precipitation itself is more important than the actual amount and that nighttime precipitation shows a stronger relationship with the UHI intensity than 24-hour precipitation (Fig 2F).
With the variables mentioned above, MLR modeling was conducted. The results showed that global solar radiation, winds from South (SWI), relative humidity and the Boolean Land use. For this analysis, the mean UHI intensity per station over the 276 days is calculated (n = 61). MLR analyses were conducted for the four models incorporating different CAD variables. However, all models include the same land use variables and the performance is very similar, with R 2 values between 0.86 and 0.87 and RMSEs between 0.31 and 0.32 K (Table 3, Models FA, FAB, TPI, REH).

Combining meteorology and land use
The combined models reach R 2 values from 0.65 to 0.71 and RMSEs from 0.69 to 0.76 K ( Table 4). The performance of the models increase with their complexity, while the differences between the CAD variables remain small (Tables 4 and 5). Irrespective of the model complexity, strongly negative UHI intensities seem to be poorly estimated by all models (Fig 3).

Mapping the urban air temperature distribution of two different nights
The predicted nighttime air temperatures of two meteorologically diverse nights are compared in order to investigate differences in model performance related to different model structures (Figs 4 and S2-S4) and CAD variables (Figs 5 and S5-S10). Favorable conditions to form large UHI intensities were experienced during the night of the 5 th of August 2018: After a sunny day (G = 304.9 Wm -2 ), a dry and calm (FF = 0.73 ms -1 ) night followed, with winds from northeast. While the solar radiation was similar during the day of the 9 th August of 2019 (G = 271.6 Wm - Table 3. Structure of the independent meteorological (MET) and the land use models with the four different CAD proxies tested (FA, FAB, TPI and REH; first step of the modeling process). Indicated are the number of observation of each model (n) and the corresponding R 2 and RMSE values. SWI stands for winds from south and NR for the Boolean nighttime precipitation variable. YES means it rained and the value is 1, otherwise it is 0. While the air temperature distribution and the performance of the models are similar during the night with favorable conditions (Fig 4A, 4C and 4E), differences appear between the model structures for the night with unfavorable conditions to form large UHI intensities ( Fig  4B, 4D, and 4F). A strongly pronounced UHI is predicted by the ADD model (Fig 4B), while almost no temperature variability is estimated by the MULT model (Fig 4D). A weak UHI is predicted by the INT model, which results in a better model performance for that night (r: 0.90 compared to 0.80 (ADD) and 0.73 (MULT); Fig 4B, 4D, and 4F).

Model
The predicted nighttime air temperatures of the INT models look very similar for all CAD variables with differences being mostly smaller than 0.5 K (Figs 4, 5 and S2-S10). Larger differences are observed during the night of the 5 th of August 2018 close to the hills and valleys outside of the city, such as Gurten hill and its adjacent valley in the south (Köniztal; Figs 1B and 5). While the slopes of the hills are estimated to be warmer by the REH model (Fig 5A-5C), the temperatures of Table 5. Structure of the combined TPI models (second step of the modeling process). The models of the other CAD variables can be found in the supplementary material (S1-S3 Tables).   the valleys (as well as the river) are modeled lower with the TPI model (Fig 5A, 5D and 5E). Several cold air flows are modeled by the flow accumulation based models, resulting in these locations appearing cooler in those models (Fig 5B-5E). However, the performance of all CAD models is virtually the same for the two days mapped and for the calibration period.

Model evaluation
Additional stations in 2018. Applied to the data of the 17 additional stations in 2018, the performance of the models is similar than with the calibration data (Tables 4 and 6). The UHI intensity of these additional stations are in general slightly underestimated (MB: -0.21 to -0.29 Model stations during summer 2021. Applied to data of the 55 remaining stations in 2021, the model performance is slightly worse compared to the calibration data (Table 4) and the UHI intensities of all stations are overestimated by 0.41 to 0.43 K (Table 6). While the CAD variables lie close to each other, the performance of the model structures increases with increasing complexity.

Evaluation over small-scale complex topographic features
Gurten hill 2021 and Egelsee depression 2019. The evaluation on the three additional stations at Gurten hill in 2021 shows that the air temperatures are (slightly) overestimated by all models. While the positive mean biases and the RMSEs of the MULT models are large, the INT models show the highest R 2 , the lowest RMSEs, and only very small positive mean biases ( Table 7). The differences between the CAD proxies are small with REH models showing slightly lower model fits (RMSE 1.12 to 1.88 K compared to 1.02 to 1.74 K).
Applying the LUR models on the additional data at the topographic depression "Egelsee" showed that none of them could adequately represent this small-scale cold air pool, which can reach up to -3 K during favorable conditions. The air temperatures of the three stations are overestimated by more than 1 K on average and the variance could not be explained (R 2 = 0.00). All CAD proxies perform similarly poor (Table 7).

Estimation of nocturnal air temperature variability using LUR models
The present study aims to model nocturnal UHI intensities and map summertime nighttime air temperatures in a city with complex terrain using a LUR approach, whereas three different model structures and four CAD variables are tested. Hereby, one crucial step is the selection of the critical meteorological and land use variables. It has been shown in various studies, that large UHI intensities are found during "calm and clear" (low wind speed and clear sky) nights after radiation intense days [4,[18][19][20][21]. This research leads to similar findings for the city of Bern, but with the addition, that wind direction may also play an important role depending on the location of the rural reference station. Here, the rural reference station is located in the northeast of the city (Fig 1B), and if the winds originate from this direction, air masses close to the reference station were likely to be exchanged faster than in the city, due to blocking by buildings or orography. Hence, the temperature differences between urban sites and the rural reference station in Zollikofen are even larger during nights with northerly winds than during calm nights. Conversely, winds from southeast or southwest lead to strongly declining UHI intensities ( Fig 3E). Additionally, little attention has so far been given to the inclusion of precipitation variables in similar LUR modeling studies, since the focus has mainly been on conditions favoring large UHI intensities [4,12,14]. However, the weakening effect of precipitation on the UHI intensity in general is well documented [19,22] and should be incorporated if longer timespans are investigated. Here, a Boolean precipitation variable instead of a numerical variable is used, which can be justified by the fact that precipitation events during the night likely vanishes a substantial share of the urban heat, independently of the amount of precipitation.
The land use variables that are significant in this study differ slightly from a previous study in the same city, which only focused on land use variables and heatwave episodes [17]. While the main pattern of open and vegetated areas having a cooling effect is similar, the cooling effect of forested areas and dense vegetation is more explicit in this study. Contrariwise, the cooling effect of water seems to be less pronounced if the entire summer is analyzed, which might be due to the water temperatures not significantly being cooler than the air temperatures during the night. Within the densely built areas, unsealed garden areas and vegetation lead to lower local temperatures (Figs 4 and S2-S9). Another important difference compared to the previous study is that altitude is directly incorporated in the model with the variable AD. This allows for a modeling of warm hilltops and slopes during calm and clear nights (Fig 4A, 4C and 4D), which is an important feature since such inversion patterns are a typical characteristic of the local climate of Bern [37]. However, when analyzing such rural, hilly areas, it is important to keep in mind that only a few rural stations were available for model calibration and validation, making uncertainties larger in these locations of the study area.
The evaluation with additional data (Table 6) shows that the model performance (R 2 and RMSE) is similar (summer 2018) or a bit lower (summer 2021) compared to the calibration data ( Table 4). The air temperatures are rather underestimated during 2018 (MB -0.21 to -0.29 K) and overestimated during 2021 (MB +0.41 to +0.43 K; Table 6). One likely reason for these differences can be found in the fact that Switzerland experienced a very hot and record dry summer 2018 and a rather cool and very wet summer 2021 [52,53]. The differences in air temperature and precipitation might thus have led to a very different state of the vegetation during these summers, having larger evaporative cooling capacities during summer 2021, which may have caused lower UHI intensities. In this study, the land use variables including vegetation are treated as static and thus do not account for changes due to inter-annual variability of meteorological background conditions. An inclusion of such an additional variable accounting for the year to year changes in the state of the vegetation (e.g. NDVI) would be an interesting supplement for future studies.
The presented LUR models are subject to additional limitations being worthwhile to be mentioned. As such, shallow depressions may cause air temperature differences of several K during nighttime [25], even within a city. The application of the models on a shallow depression shows that none of the models is able to represent such small-scale air temperature differences in an adequate manner (Table 7). This might be due to the rather large buffer radii of the land use variables (up to 1000 m; Table 1) or the cooling effect of the CAD variables being to weak: A detailed analysis of the CAD variables shows that they are only able to capture about 0.6 K of the cooling in the Egelsee depression, whereas differences of up to 3 K are reached under favorable conditions. Another limitation is that the models tested here are static and do not take into account weather patterns from previous days and resulting variables depending on those, such as soil moisture, which may influence the present day air temperature patterns [20]. Lastly, the exposition of a location to the wind is not incorporated, since no land use variable was found which was able to represent this aspect.

Differences between model structures and CAD variables
Three model structures with different degrees of complexity and four straight forward computable variables to incorporate cold air dynamics induced by topographic features are tested in this study. Regarding the model structures, the MULT models shows inaccuracies if either a station is not considerably warmer (or even cooler) than the reference station, or if a day with meteorological conditions unfavorable for the formation of the UHI is modeled. In both cases, one of the factors in the equation of the MULT model is found to be close to zero or negative (Table 2), and the resulting variability of urban air temperatures is not realistic (Fig 3B). Hence, this model structure is only appropriate if "cool" reference and "warm" urban stations are used and only when days with meteorological conditions favoring large UHI intensities are analyzed. Besides that, enhanced complexity leads to better performance. The inclusion of the interaction terms leads to an increase in R 2 between 0.057 and 0.058 (depending on the CAD variable), which implies that the interaction effect accounts for 5.7 to 5.8% of the variance [50].
In the simple ADD model, the magnitude and the shape of the UHI remains similar for every night (Fig 4A and 4B), whereas the effect of different meteorological conditions on the formation of urban heat islands (and also CAD) can be incorporated with the interaction terms.
The resulting air temperature fields of all CAD variables show more realistic patterns than the simple flow accumulation approach of an earlier study [17] (Figs 4 and S2-S9). The largest differences between the CAD variables can be found in areas where strong CAD is expected. While models based on elevation based CADs model the typical pattern of temperature inversions in a more pronounced way, flow accumulation based CADs are able to model cold air paths along the topography (Fig 5). Due to the limited number of stations in the hilly areas and the similar performance in calibration and evaluation, no statement about which CAD variable to be the most appropriate can be made here. Additionally, it has to be stated that none of the CAD variables is capable to model cooling effects of a small-scale depression within the city ( Table 7), and that potential blocking effects of buildings or trees are not incorporated by three out of four CAD variables (FA, TPI and REH). Additional research and the evaluation of these CAD variables in other cities with complex terrain would hence be beneficial.

Conclusions
The aim of this study was to model nocturnal UHI intensities and to map nighttime air temperature fields in the city of Bern during summers 2018 to 2020 with a LUR approach based on air temperature data of a low cost measurement network and publicly available geospatial and meteorological data.
The results show that, in addition to precipitation, global solar radiation and wind speed, wind direction is an important factor for the distribution of the nocturnal urban air temperatures variability in Bern. Further, we tested three different model structures to combine meteorology and land use and four straight-forward computable variables to predict cold air drainage. The resulting 12 models reached good model performances (R 2 0.65 to 0.71; RMSE 0.69 to 0.76) which could be validated with additional data (R 2 0.59 to 0.73; RMSE 0.77 to 0.85). While multiplicative model structures were not able to realistically predict the air temperatures of cool stations and for days with meteorological conditions unfavorable for the formation of urban heat, additive and integrative model structures performed well, of which the latter is favored due to the possibility of modeling different magnitudes of urban heat island intensity.
With the derived models, it is possible to calculate and map the air temperature distribution of a city in complex terrain for every night of a summer using meteorological data of only one official measurement station and publicly available land use data. With this method, only a very low amount of computational (about 5 min on a basic computer for 92 days) and financial (publicly available data) effort is needed to gain valuable information about the intra-urban temperature distribution in a city which can then be used for various applications in urban planning, heat risk management and future research.