City-wide, high-resolution mapping of evapotranspiration to guide climate-resilient planning

The impacts of global change, including extreme heat and water scarcity, are threatening an ever-growing urban world population. Evapotranspiration (ET) mitigates the urban heat island, reducing the effect of heat waves. It can also be used as a proxy for vegetation water use, making it a crucial tool to plan resilient green cities. To optimize the trade-off between urban greening and water security, reliable and up-to-date maps of ET for cities are urgently needed. Despite its importance


Introduction
Climate change impacts are fundamentally threatening the way of life of residents in many urban areas, where risks of heat waves, flooding, and water scarcity are exacerbated (Kundzewicz et al., 2014;Li and Bou-Zeid, 2013;Manoli et al., 2019). These hazards are leading to an increase in casualties, adverse health effects, and infrastructure damage (Kundzewicz et al., 2014;Mora et al., 2017). As a growing majority of the world population resides in cities (United Nations, 2019), urban sustainability measures may have the greatest impact on human qualityof-life (Acuto et al., 2018;Gaffin et al., 2012). To meet UN Sustainable Development Goal 11 of making cities more sustainable and resilient, urban planners and policy-makers urgently need to adopt mitigation measures grounded in evidence-based understanding of urban (eco) systems (Bayulken et al., 2021).
ET is a key indicator for urban sustainability, as it is associated with thermal comfort (Shashua-Bar et al., 2011;Wang et al., 2019), water storage capacity (Jongen et al., 2022), the irrigation demand of green spaces (Litvak et al., 2017;Saher et al., 2021), and stormwater runoff reduction (Berland et al., 2017;Wan et al., 2021). An improved understanding of urban ET and its spatio-temporal patterns is essential to future-proof cities worldwide against climate risks (Berland et al., 2017;Saher et al., 2021). Sustainability initiatives in urban areas need to be adapted to a neighborhood scale relevant for vulnerable sociodemographic groups (Hölzl et al., 2021). Therefore, for sustainable management of resilient cities, accurate, high-spatial and -temporal resolution maps of ET on a city scale are needed.
However, urban greening programs still require more reliable data and models to quantify the ecosystem services of urban vegetation and their spatio-temporal variability (Pataki, 2013). In recent years, machine learning (ML) has emerged as one of the most promising methods to extract knowledge from big Earth system data, which comprises a variety of remote sensing and in situ sensors (Reichstein et al., 2019). Regression ML algorithms have been especially successful in predicting the spatio-temporal variation of biogeophysical parameters by linking point observations with satellite data (Reichstein et al., 2019). Datadriven methods provide numerous advantages, including the capacity to discover patterns and infer from data without process-based formulation (Kuhn and Johnson, 2013;Reichstein et al., 2019).
Despite the relevance of ET estimation in urban areas, urban ET has rarely been modeled and even more rarely mapped (Nouri et al., 2015;Saher et al., 2021;Vulova et al., 2021). Characterizing ET is always challenging, but in particular in an urban environment due to a high heterogeneity of land cover, plant traits, soil characteristics, and microclimate (Nouri et al., 2013b;Saher et al., 2021). Therefore, process-based models have mainly been used to model urban ET, with the most common types being surface energy balance (SEB) models, hydrological models, and urban land-surface models (ULSMs) (Järvi et al., 2011;Karsisto et al., 2016;Nouri et al., 2015;Rafael et al., 2020;Saher et al., 2021;Ward et al., 2016). As the majority of these models have been developed for agricultural and natural areas, their applications in urban areas have been limited in reliability and transferability (Nouri et al., 2015;Saher et al., 2021). Medium-resolution earth observation data, such as Landsat and MODIS, have been integrated into SEB models to improve urban ET mapping (Abunnasr et al., 2022;Cong et al., 2017;Faridatul et al., 2020;Jiang and Weng, 2017). A SEB model modified for urban areas (SEBU) has recently been implemented in Google Earth Engine, which can characterize ET dynamics in cities at a 100-m spatial resolution and monthly temporal resolution (Abunnasr et al., 2022). However, most SEB models still do not perform well in characterizing urban ET, as they exclude anthropogenic heat, assume a single plant species, and are optimized for a regional scale (Saher et al., 2021).
Anthropogenic heat and moisture emissions are associated with energy consumption in cities, with the main sources being commercial and residential buildings, vehicles, and industry (Sailor, 2011). Anthropogenic heat warms the urban surface and atmosphere, contributing to the UHI effect, with a proportionally larger influence in wintertime (Bohnenstengel et al., 2014) and nighttime (Narumi et al., 2009). Anthropogenic heat is highly spatially and temporally variable, with a strong link to human activity and population density (Narumi et al., 2009;Rapsomanikis et al., 2015). Thus, satellite imagery, such as land use and land cover products, can be used as a proxy of anthropogenic heat when estimating urban ET (Cong et al., 2017;Faridatul et al., 2020). The three main types of approaches for estimating anthropogenic heat and moisture emissions are energy budget closure methods, building energy models, and inventory approaches (Pigeon et al., 2007;Sailor, 2011). More anthropogenic heat generally leads to higher ET values (Rapsomanikis et al., 2015). However, anthropogenic heat is omitted from most physical ET models, leading to an underestimation of ET (Allen et al., 2011a;Cong et al., 2017;Faridatul et al., 2020).
Furthermore, physical modeling approaches to quantifying ET developed for urban areas commonly rely on parameters which are impractical to estimate on the city scale at a high resolution, such as soil water potential, wall emissivity, and heat emissions from vehicles (Meili et al., 2020;Saher et al., 2021). For urban planning, regularly updated city-level ET maps at a hyperlocal resolution (10-30 m) (Venter et al., 2020) are needed to assess high-risk areas related to extreme heat and stormwater runoff and to propose integrated solutions, such as street trees and green roofs (Pataki, 2013;Wan et al., 2021). The few studies mapping urban ET so far have been limited in temporal resolution (Abunnasr et al., 2022), extent (Rocha et al., 2022b;Vulova et al., 2021), spatial resolution (Saher et al., 2021;Wang et al., 2016), or accuracy Grimmond et al., 2010;Rafael et al., 2020;Wan et al., 2021). ML, in combination with open-access remote sensing data and in situ measurements, offers a way forward to accurately mapping spatio-temporally dynamic ET in urban areas (Reichstein et al., 2019;Vulova et al., 2021).
Flux towers are increasingly providing direct measurements of ET in cities worldwide (Jongen et al., 2022), many of which have recently begun to provide harmonized data free of charge (Lipson et al., 2022). The expanding access to high-quality urban flux data opens up the possibility to model urban ET empirically. However, to the best of our knowledge, only one study applied an empirical model to map ET in an urban area so far (Wang et al., 2016). Wang et al. (2016) related a buffer average extracted from two MODIS products (albedo and land surface temperature) around a single flux tower site in Phoenix, USA using ordinary least squares (OLS) regression. This regression model was used to map daily ET at a 500-m spatial resolution (Wang et al., 2016). However, the capacity of mapping ET in a city using a data-driven approach was limited by the lack of ET data from various urban land covers (Wang et al., 2016). To address this issue, Vulova et al. (2021) developed a novel approach to model urban ET using flux footprint modeling to incorporate a variety of land cover into artificial intelligence algorithms. However, this methodology was only applied to model ET at two points in a city (Vulova et al., 2021).
Here, we map urban ET at the city scale at a high spatial (10-m) and temporal (hourly) resolution using a novel data-driven approach. Citywide hourly ET maps were generated for the entire year of 2019, which were aggregated to daily, monthly, and annual ET maps. We integrated the influence of heterogeneous land cover on urban ET into ML models by combining open-access remote sensing data and flux footprint modeling. Hourly (10-m resolution) maps for a 2 km 2 extent around the flux towers for two years (2019-2020) were used for validation with flux tower measurements at two sites in Berlin, Germany. The city-wide ET maps were further analyzed using Local Climate Zones (LCZs), a classification developed to better understand urban heat island (UHI) effects, to explore the implications for urban planning. Our key objectives are to: 1) establish a method for estimating urban ET at a high spatial and temporal resolution and to produce a dataset of ET in Berlin and 2) analyze the spatial and temporal variability of urban ET in Berlin.

Study area
ET was mapped for the entire city of Berlin, the largest and capital city in Germany with 3.7 million inhabitants and 891 km 2 (Statistical Office of Berlin-Brandenburg, 2019) (Fig. 3b). Berlin is a mid-latitude city characterized by a climate transitioning between a temperate oceanic and a humid continental climate, with exceptionally hot days in summertime (Beck et al., 2018). The city has a flat topography with an average elevation of 34 m above sea level. Daily air temperature (Tair), shortwave downward radiation, precipitation (P), and NDVI during the study period (2019-2020) are shown in Fig. 7 (Supplement). Both years were warmer than the long-term (1981-2010) mean (9.7 • C), with an average Tair of 11.8 • C and 11.7 • C in 2019 and 2020, respectively (DWD, 2021) and drier than the long-term mean (591 mm) (DWD, 2021), with a P of 520 mm and 479 mm in 2019 and 2020, respectively.
Berlin was selected as a pilot for modeling ET as the city is a hotspot for heat stress-related fatalities (Gabriel and Endlicher, 2011;Krug et al., 2019), with 5% of all deaths between 2001 and 2010 associated with elevated air temperatures (Scherer et al., 2013). Approximately 34% of Berlin's land surface is impervious (Senate Department for Urban Development and Housing, 2017), which is largely concentrated around the city center (Fig. 3b). Heat stress risk is most pronounced in the midrise dense residential areas in the city center and high-rise districts in the eastern part of the city, where higher air temperature coincides with a high proportion of elderly residents (Dugord et al., 2014). At the same time, Berlin has been identified as a pioneer city worldwide in the use of green infrastructure to reduce its heat risk and water footprint (Liu and Jensen, 2018). According to the European Environment Agency, Berlin has one of the highest percentages of 'green infrastructure' (51%), a term encompassing both green (vegetated) and blue (water) spaces, among the 38 European capital cities (EEA, 2022). Surface water accounts for 6.6% of the city's area (SenUVK, 2021). Two flux towers, which are part of the Urban Climate Observatory (UCO) Berlin , were used for training and testing RF models. The surroundings of the towers contain heterogeneous vegetation (grassland, shrub, and trees) and impervious cover representative of land cover in temperate-zone cities (Gillefalk et al., 2021;Jongen et al., 2022).

Evapotranspiration mapping
We mapped the ET of Berlin at a high temporal (hourly) and spatial (10-m) resolution using Random Forest (RF), open data available for most large-and medium-size cities worldwide, including remote sensing imagery, geodata, and meteorological data, and eddy covariance (EC) data (Table 1). An overview of the methodology is provided in Fig. 1. We integrated the variability of urban land cover into RF models by estimating the extents of constantly-changing footprints of two flux towers using flux footprint modeling (Kormann and Meixner, 2001). RF models were trained in one year, while testing and mapping was performed with the other year (e.g., training in 2020 and testing in 2019) and vice versa. Measured ET from the two flux towers was used to estimate the model accuracy. We validated our mapping approach by extracting footprints from hourly ET maps for the entire study period (2019-2020). We then aggregated the city-wide hourly ET maps from 2019 to generate daily, monthly, and annual maps. Finally, as ET varies greatly among different land surfaces, we assessed the variation of ET between LCZs. The codes are publicly available in a GitHub repository (Vulova, 2022).

Flux measurements
The EC instrumentation, data quality control (QC) and preprocessing are elaborated in Vulova et al. (2021). All half-hourly ET observations (in mm/ h) remaining after QC were used for training, consisting of 16,011 and 16,255 observations (45.7% and 46.3% of the raw EC data) combined from both towers for 2019 and 2020, respectively.
Footprint modeling over a geospatial database, including vegetation and impervious surface fraction, can be used to explain the variation of sensible and latent heat fluxes in urban areas (Demuzere et al., 2008;Kotthaus and Grimmond, 2014;Offerle et al., 2006). We used the Kormann and Meixner (2001) model to estimate the footprints for each half hour during the study period. We selected the Kormann and Meixner (2001) footprint model as it (1) has a low computational cost (Kotthaus and Grimmond, 2014;Neftel et al., 2008;Schmid, 2002), (2) is available through open-source tools (Neftel et al., 2008;Xenakis, 2016), (3) is easy-to-implement as the input parameters are provided by the EC measurements and open-source GIS data (Xenakis, 2016), (4) is one of the most commonly used footprint models in urban areas (e.g., (Abunnasr et al., 2022;Christen et al., 2011;Chrysoulakis et al., 2018;Crawford et al., 2011;Crawford and Christen, 2015;Frey and Parlow, 2012;Kotthaus and Grimmond, 2014;Liu et al., 2012;Salmond et al., 2012;Vulova et al., 2021)), (5) contains realistic profiles of eddy diffusivity and wind speed (Schmid, 2002), and (6) has shown no difference in accuracy compared to complex Lagrangian stochastic dispersion models in field tests (Coates et al., 2021). In short, this footprint model is robust (Coates et al., 2021), commonly applied in urban areas (Abunnasr et al., 2022;Chrysoulakis et al., 2018;Kotthaus and Grimmond, 2014) and fulfills our aim of reproducibility, as non- Table 1 Predictor variables used to model and map evapotranspiration (ET); the data source, provider, and their original spatial and temporal resolutions are listed. In the "Temporal resolution" column, "static" refers to geodata which does not vary temporally (such as impervious surface fraction). 1 Static * In the "Spatial resolution" column, "polygons" refers to datasets originally provided as vector data, which were then rasterized to 1-m resolution. specialists or those with restricted computing capacity (e.g., in resourcelimited countries or institutions) can also apply the presented approach in future studies. We integrated the heterogeneity of the urban environment into footprint models using a sophisticated approach incorporating wind direction, 3D geodata, and vegetation phenology. The aerodynamic parameter zd (representing the height above ground where the wind speed is zero due to obstacles to the flow) used in footprint modeling was estimated in an iterative process for each 30-min interval (Grimmond and Oke, 1999). Morphometric parameters for each site were computed using the Urban Multi-scale Environmental Predictor (UMEP) QGIS plug-in  with a 1000 m radius and a 1 • wind direction search interval. Three different rasters were used as input for UMEP: VH, BH, and a raster combining VH and BH (representing the maximum height of either VH or BH at any given pixel). For this analysis, BH and VH were both acquired from the Berlin Digital Environmental Atlas at 1-m resolution (Berlin Senate Department for Urban Development and Housing, 2014). Zd was calculated with the method of Kanda et al. (2013). Using zd only calculated from BH would underestimate zd, as trees increase the aerodynamic roughness. Therefore, the method from Kent et al. (2017) was used to incorporate vegetation as a roughness element. A different aerodynamic porosity value was used for each season (winter, summer, and intermediate) to account for plant phenology (Kent et al., 2017). The seasonally corrected one-degree zd values were then averaged over a running mean including the zd values from wind direction i ± 5 • moving forward by one degree. These zd estimates were assigned to the EC data based upon wind direction and season and used to calculate flux footprints. The parametrization of footprint models is further described in Vulova et al. (2021).

Remote sensing and GIS data
Predictors were selected based on (1) their relevance as a driver of urban ET based on physical principles and a comprehensive literature review and (2) the availability of open-access data in order to enable reproducibility. Normalized Difference Vegetation Index (NDVI) (Tucker, 1979) is an indicator of vegetation health, density, and biophysical characteristics which is well-suited for modeling urban ET (Nouri et al., 2015;Saher et al., 2021;Vulova et al., 2021). The advantage of using NDVI to model ET is that it captures the temporal variability of plant traits due to phenology, water availability, and disease (Damm et al., 2018;Nouri et al., 2015;Saher et al., 2021), while vegetation fraction is a static dataset. Furthermore, NDVI with a 60 m or higher spatial resolution has been found to be particularly useful for modeling the spatial distribution of energy and heat fluxes in urban areas (Rigo and Parlow, 2007). We used the Google Earth Engine platform (Gorelick et al., 2017) to calculate and download NDVI raster time series from Sentinel-2 Level-2A products at a 10-m resolution (image collection "COPERNICUS/S2_SR"). Atmospheric, terrain, and cirrus correction of Sentinel-2 images was performed using the Sen2Cor atmospheric correction processor. Pixels remaining after cloud masking were interpolated to a daily resolution as NDVI is assumed to be stable on a sub-daily scale. Building height (BH) has a major influence on urban microclimate and therefore ET as a driver of wind flow and shade (Masson et al., 2020;Saher et al., 2021). Urban Atlas 2012 (European Environment Agency, 2018) provided BH at a 10-m resolution. Vegetation height (VH) is an indicator of several factors influencing plant transpiration, including the vegetation type (trees, shrubs, or grasslands), the distance water travels from the roots to the leaves, and the gravitational hydrostatic gradient (Damm et al., 2018;Nobel, 2009). VH was provided by the Berlin Digital Environmental Atlas at 1-m resolution (Berlin Senate Department for Urban Development and Housing, 2014). Vegetation fraction (VF) has been found to be well-linked to heat and energy flux variation in urban areas (Offerle et al., 2006). VF was provided by the Berlin Digital Environmental Atlas (Berlin Senate Department for Urban Development and Housing, 2014). Impervious surface fraction (ISF) is strongly linked to the partitioning of convective energy fluxes (Offerle et al., 2006), indicating thermal emissions from paved surfaces (Saher et al., 2021) and anthropogenic heat (Faridatul et al., 2020). ISF was derived from the Urban Atlas 2018 Land Use and Land Cover (LULC) dataset (European Environment Agency, 2018), as explained in Vulova et al. (2020). Footprint modeling allowed for the remote sensing and GIS datasets (Table 1) to be converted to a halfhourly temporal resolution (matching the EC data).
The GIS data (Table 1) was selected to be as temporally close to the mapping year (2019) as data availability allowed. The most recent data was acquired in 2018 for ISF (European Environment Agency, 2018), 2012 for BH (European Environment Agency, 2018), and 2009-2010 for VH and VF (Berlin Senate Department for Urban Development and Housing, 2014). However, urban form and cover in Berlin is assumed to be fairly stable in this time frame, as urbanization has been found to progress slowly in European cities (Sun et al., 2020). ISF and VF, which were available as polygons, were rasterized to a 1-m resolution. Finally, all geodata were resampled to a common (10-m) resolution before footprint modeling and ET mapping. Using 10-m resolution remote sensing and GIS data is necessary to capture the variability of land cover covered by footprints, as a coarser spatial resolution would decrease the number of unmixed pixels overlapping with footprints, reducing the capacity to identify the patterns in the urban surface that determine ET (averaging effect). A 10-m spatial resolution was selected for mapping ET as (a) urban planners and policy-makers require information at a hyperlocal (10-30 m) resolution to develop neighborhood and household-scale strategies for climate change mitigation (Milojevic-Dupont and Creutzig, 2021;Venter et al., 2020), (b) the highest resolution of available remote sensing and GIS data could be exploited to characterize urban heterogeneity, and (c) a 10-m resolution still incorporates mixed land cover, as is the case with the footprint-derived training data. Water bodies were excluded from footprint modeling and ET mapping, as there were nearly no water bodies present in the training regions (footprints). Weighted averages of the surface cover were then computed by multiplying the footprint grids with the raster layers and summing the pixel values on a half-hourly basis.

Machine learning approach
To map ET, we used the Random Forest (RF) ML algorithm, a supervised regression model which combines classification and regression trees (CART) with bagging (Breiman, 2001). Detailed descriptions of RF and its implementation using the "caret" R package are provided in Kuhn and Johnson (2013). RF models were trained using half-hourly data (the resolution of EC data and footprint modeling) combined from both towers in order to include a greater variety of combinations of land cover covered by footprints and meteorological conditions. The wide range of land cover integrated into RF models (Supplement; Fig. 11) allowed for ET to be predicted for the entire city of Berlin.

Improvement of existing methods for estimating urban evapotranspiration
While most previous studies used process-based models to estimate urban ET at a plot size or a coarse spatial resolution, in this study urban ET is mapped with a data-driven approach using medium-resolution satellite imagery. Vulova et al. (2021) used artificial intelligence models to estimate urban ET on a point basis. However, sustainability initiatives in cities need to be implemented at a neighborhood scale capturing socio-demographic vulnerability to heat stress (Milojevic-Dupont and Creutzig, 2021). Thus, ET estimation on a point basis is insufficient for urban planners and climatologists preparing for the effects of climate and global change. The major gap in the field of urban ET is that high-resolution, spatially distributed data for an entire city are generally unavailable (Saher et al., 2021).
In order to fill this gap and map the ET of Berlin at a high spatial and temporal resolution, it was necessary to further develop the methodology provided by Vulova et al. (2021), as follows: 1) producing spatially distributed meteorological data, 2) combining the EC data from the two flux towers to create a general model, and 3) validating the result in a spatially explicit manner. Lastly, a major innovation of this study is the ET dataset itself, which allowed for the variation of ET in an entire city to be characterized at a high spatial and temporal resolution. This ET product allowed for the variation within and between LCZs to be analyzed, which can contribute further insights for urban planners and climatologists. The following subsections focus on the main novel aspects of our methodology.

Spatially distributed meteorological data
Reference ET (RET) was selected as a predictor to synthesize the meteorological conditions regulating ET (Allen et al., 2005). For agricultural applications, RET is corrected to actual ET using crop coefficients (constant), which characterize the influence of crop characteristics on water demand (Allen et al., 1998). With the presented machine learning approach, models automatically scale RET based on the remote sensing and GIS predictors to predict actual ET. We used open-access hourly meteorological data from the German Meteorological Service (DWD) to calculate RET (DWD, 2021). From all the DWD stations within 70 km of Berlin, 10 stations were selected based upon the availability of all the data needed to compute RET (air temperature, relative humidity, wind speed, and sunshine duration). Solar radiation was estimated from sunshine duration using the extraterrestrial equation from the MeTo package (ra function) (Dettmann and Grimma, 2019) due to only one DWD station within 70 km having solar radiation data. RET was computed using the hourly ASCE "Standardized Reference Evapotranspiration Equation" for short crops (Allen et al., 2005). The ASCE equation was selected to estimate RET as it (1) incorporates the most important meteorological parameters to accurately estimate RET in an urban environment in comparison to alternative equations relying on limited data (Shafieiyoun et al., 2020), (2) has been identified as well-suited for urban areas in several studies (DiGiovanni-White et al., 2018;Marasco et al., 2014;Nouri et al., 2013a), (3) is highly recommended for calculating RET for consistency and reproducibility due to its global standardization (Allen et al., 1998(Allen et al., , 2005Allen et al., 2011b;Nouri et al., 2013a), and (4) is implemented at an hourly time step (Allen et al., 2005). For further details on the calculation of RET, refer to Vulova et al. (2021). RET was then linearly interpolated to a half-hourly scale for training to match the temporal resolution of the EC data. As spatially continuous rasters of all predictors are necessary for ET mapping, RET was interpolated with ordinary kriging at a 1 km resolution using the R package "gstat" (Gräler et al., 2016). Ordinary kriging was selected as it has been shown to be a suitable and accurate technique for the spatial interpolation of RET (Ashraf et al., 1997;Dalezios et al., 2002;Hodam et al., 2017;Kamali et al., 2015;Mardikis et al., 2005;Martinez-Cob, 1996;Rocha et al., 2022a). For training, RET was extracted from the kriged rasters at the coordinates of the two towers. The kriged RET was validated with RET calculated with independent meteorological data measured at the flux towers, showing a very high correlation (Pearson's r = 0.97) at both towers (see Supplement, Fig. 14).

Training a general model
In Vulova et al. (2021), EC data from a single tower is used to train models and test their accuracy locally. In this study, in order to fit general models incorporating a wide range of land cover, the EC data of the two flux towers in Berlin were combined. This change was crucial for generalizing the models and predicting ET for the entire city rather than for each tower's footprint for each timestamp.

Spatially explicit validation
Using spatially dynamic footprints of two EC towers is the best available solution for validating urban ET maps, as other sources of in situ data, such as lysimeters and sap flow sensors, are point-based and largely unavailable in cities Nouri et al., 2013b). Even when EC data is available, studies modeling ET on a point basis (e.g. (Vulova et al., 2021)) cannot incorporate the changing source area of EC measurements into the validation process. We mapped ET on an hourly basis based upon the highest temporal resolution of predictors. For testing, ET was mapped on an hourly basis at a 10-m resolution for 2 km × 2 km extents (the likely maximum coverage of flux footprints) centered around the flux towers for 2019-2020. To perform validation on an hourly basis, the quality-controlled EC data was averaged to an hourly resolution, omitting hours where one or no half-hourly observations were available. The output of footprint modeling are grids with a 10-m resolution for the 1000-m surrounding area, where each grid cell holds the probability of contribution of that location and time (Supplement; Fig. 15). The grid cells out of the 90% footprint probability are considered to have zero contribution. As the footprint is weighted and varies in size and direction, the ET product can be validated in a robust and spatially dynamic manner. However, this validation approach still cannot cover all types of urban surfaces. A dense, spatially distributed network of sensors such as lysimeters or flux towers would improve the reliability of validation, but such data is currently unavailable for nearly all cities worldwide (Saher et al., 2021) and have their own limitations and biases (Nouri et al., 2013b). Weighted averages from the 2 km × 2 km ET maps were extracted using flux footprints modeled at a 10-m resolution to compare modeled with measured ET. To evaluate the temporal extrapolation capacity of our models, we used 2020 data for training and 2019 data for testing and vice versa. In total, 6696 and 6771 hourly ET observations combined from both towers were used for testing in 2019 and 2020, respectively. The prediction accuracy was evaluated using root mean square error (RMSE), percent bias (pbias), coefficient of determination (R 2 ), and normalized root mean square error (NRMSE).

Local climate zone analysis
LCZs refer to a classification system initially developed to provide a universal framework for UHI studies, consisting of 17 classes of which 10 are urban-specific (Stewart and Oke, 2012). Each LCZ has a characteristic air temperature regime based on surface structure and land cover (Stewart and Oke, 2012). However, air temperature within LCZs can also be highly variable, especially at nighttime ; an extensive analysis of the intra-LCZ air temperature variability in Berlin is provided by Fenner et al. (2017). For an analysis of the homogeneity of LST within LCZs, refer to the Supplement (Fig. 12). The main advantage of the LCZ scheme is comparability, where a LCZ in Tokyo, Delhi, or Paris refers to the same urban form and cover (Demuzere et al., 2021;Stewart and Oke, 2012). LCZs are also advantageous for guiding urban planners, architects, and engineers to construct climate-adapted cities due to simple and clear definitions of the design elements of each LCZ (e. g., building height, building density, and vegetation cover ratio) (Stewart and Oke, 2012). Many of the LCZ classes are defined as mixed in land cover and form; for instance, the "Open midrise" LCZ contains a mixture of midrise buildings, grass, trees, and impervious surfaces (Stewart and Oke, 2012).
LCZs are not used to model and map ET in this study, but rather to aggregate the modeled ET in a manner useful for urban planners and climatologists. Aggregating modeled ET per LCZ has two essential applications: (1) providing urban planners and policymakers with ET estimates that can be linked to surface structure and land cover, which can guide climate-resilient planning and development (Aslam and Rana, 2022), and (2) supplying modelers with ET input values grouped per LCZ, which can improve the parameterization of process-based models (e.g., urban flood models) (Liu et al., 2020). Deriving the annual or summer ET for each LCZ directly in Berlin using only observational data is not possible as not all LCZs are represented near the towers (Supplement; Fig. 10) and the data contain gaps. To assess and demonstrate the variation of ET for different LCZs across Berlin, we used a LCZ map of Berlin produced by the LCZ Generator (Demuzere et al., 2021) using training data provided by Fenner et al. (2017) (Fig. 3b). We determined that this LCZ map was optimal for Berlin based on visual comparison with high-resolution Google Earth imagery and personal communication with experts on LCZ classification (Bechtel et al., 2019;Demuzere et al., 2021;Fenner et al., 2017). For instance, other available LCZ products reporting higher accuracies were not utilized as they were missing LCZs known to occur in Berlin, such as "Open highrise" (Demuzere et al., 2019) and "Large lowrise" (Muhammad et al., 2022). The LCZ map showed an overall accuracy (OA) of 61% (above the 50% quality control standard suggested by Bechtel et al. (2019)) and a weighted accuracy of 92%, meaning that most confusion can be attributed to similar classes (see Bechtel et al. (2017b)). The LCZ Generator allows users to incorporate their local knowledge to map cities into LCZs in an online platform, with >1000 open-access LCZ maps of cities worldwide already posted by users (Demuzere et al., 2021). The LCZ map is representative for the years 2017-2019 based on the preprocessed earth observation data used to map LCZs (Demuzere et al., 2021). LCZ maps are generally provided at 100-m resolution (Demuzere et al., 2021), as LCZs are defined as "regions of uniform surface cover, structure, material, and human activity that span hundreds of meters to several kilometers in horizontal scale" according to Stewart and Oke (2012). As individual pixels do not comprise an LCZ class, a Gaussianfiltered version of the LCZ map, which reduces granularity, was used (Demuzere et al., 2020(Demuzere et al., , 2021. Although the Gaussian filter procedure greatly reduced the number of small (<0.04 km 2 ) individual LCZs (Demuzere et al., 2020(Demuzere et al., , 2021, some LCZs with the minimal possible extent (100 × 100 m) could still be found in the map. LCZs were converted to polygons to extract summary statistics and produce boxplots from the ET maps. Medians, rather than means, were used to describe differences between LCZs, as the mismatch between the spatial resolution of our ET maps and LCZs (10-m vs. 100-m) introduces more noise in the arithmetic mean.

Validation of the proposed method
The performance metrics when testing for each year and flux tower are shown in Table 2. On average, validation at the greener site (ROTH) showed a higher R 2 (0.76) than at the less vegetated site (TUCC; 0.56). On the other hand, RMSE was higher at the greener site (0.0289 mm) than at the city-center site (0.0171 mm). Two recent studies used a process-based model to estimate ET in Berlin, showing comparable R 2 (~0.8 at ROTH; ~0.5 at TUCC) and RMSE (~0.025 at ROTH; ~0.02 at TUCC) to our study (Rocha et al., 2022a(Rocha et al., , 2022b. The higher R 2 and RMSE when testing at the site surrounded by more vegetation (ROTH) can largely be explained by the higher range and maximum ET measured at ROTH compared to TUCC (Kuhn and Johnson, 2013;Rocha et al., 2022b;Vulova et al., 2021). Normalized RMSE (NRMSE) allows us to compare the accuracy between the two sites with a metric that removes the influence of the range of measured ET. The difference between NRMSE between the two towers was minimal (11.2% vs. 11.9% at ROTH and TUCC, respectively). ET was overestimated by a greater magnitude at the more built-up site (TUCC; pbias = 28.7%) than at the more vegetated site (ROTH; pbias = 3.1%). Training models using data from both towers led to a higher overestimation of ET at the dense city-center site where lower ET values were observed (TUCC).
Overall, averaged across the two sites, accuracy was slightly higher in 2019 (R 2 = 0.67) than in 2020 (R 2 = 0.65). ET was on average more overestimated in 2020 (pbias = 21.2%) than in 2019 (pbias = 10.6%). Training in 2019 and predicting in 2020 likely led to the overestimation as the previous year had higher air temperature (Tair) and precipitation (P) in summertime (Tair: 21.5 and 20.3 • C, P: 194 and 126 mm).
To better understand the temporal variation of the mapping error, the difference between the modeled and measured ET in mm/h (residuals) was calculated and plotted as a time series (Fig. 2). At the more vegetated site, residuals were positive in summertime for both testing years, indicating that models were underestimating ET. In contrast, at the dense city-center site, modeled ET was not underestimated in summertime, with near to zero or negative residuals. In April, RF models overestimated ET for both testing years and both sites, likely due to a time lag between the rising NDVI signal in spring and the contribution of plant leaf out to transpiration. At the same flux sites in Berlin, ET was also overestimated in April using a process-based model due to a sharp increase in vegetation greenness, which did not yet correspond to higher ET (Rocha et al., 2022b).

Spatial and temporal distribution of urban evapotranspiration
Annual estimated ET ranged from 119 mm to 648 mm (mean: 300 mm), with strong spatial patterns dependent on the land cover (Fig. 3a). Annual ET was highest in Tempelhof Field park, the city's largest open green space (Fig. 3d), in the NE of the city, which consists of cropland and pasture (Fig. 3f), and smaller green spaces such as community gardens and parks. Areas of low ET were concentrated in highly impervious areas (Fig. 3c and e), such as streets, in the city-center and the East.
The estimated monthly ET maps captured the phenology of the urban vegetation (Fig. 4). Forests in the SE and SW of the city showed very low ET in wintertime, which increased considerably during leaf out (April). Highest ET occurred in June (mean: 44 mm, max: 80 mm), while lowest ET occurred in January (mean: 11 mm, max: 39 mm).

Evapotranspiration of different urban surfaces
We analyzed the variability of urban ET within different LCZ classes both annually and in summertime (Fig. 5). ET varied highly between both the vegetated and urban LCZs, with the highest and lowest ET estimated for "low plants" and "compact midrise", respectively.
LCZ "low plants" showed the highest median ET both annually (390 mm) and in summertime (50 mm), closely followed by "bush, scrub" for both time periods. "Dense trees" had the lowest median ET (317 mm annually, 46 mm in summertime) of the natural LCZs, other than the non-vegetated "bare soil or sand" LCZ. Our study cannot fully account for interception loss, which may explain the lower estimated ET for trees. EC data during and 4 h after precipitation events is unreliable and is therefore removed, which does not allow for the evaporation from intercepted rainwater to be characterized (Vulova et al., 2021;Ward et al., 2013).
ET also showed high variation between the urban-specific LCZs. Among the urban LCZ classes, "sparsely built" showed the highest median ET both annually and in summertime (335 mm and 47 mm). "Compact midrise" showed the lowest median ET both annually and in summertime (246 mm and 29 mm).
In a few cases, the ranking of LCZs differed between annual and summer ET. For instance, "large lowrise" showed a higher median ET than "open midrise" annually, while the opposite was the case in summertime. In winter months, ET was higher in "large lowrise", a highly paved LCZ, than in "open midrise" (Supplement; Fig. 9). In addition, "bare soil or sand" showed higher ET than "open lowrise" annually, while the relationship was reversed in summertime. This discrepancy Table 2 Performance metrics for validation at two flux towers and two years. The performance metrics are root mean square error (RMSE), percent bias (pbias), coefficient of determination (R 2 ), and normalized root mean square error (NRMSE). The best performance metrics are shown in bold. can also be attributed to the higher ET of "bare soil or sand" than "open lowrise" in wintertime (Supplement; Fig. 9). Our approach can still capture the influence of interception loss occurring >4 h after rainfall. Thus, these differences can be explained by evaporation from interception loss, which is generally higher in wintertime and over impervious surfaces (Miralles et al., 2020;Moriwaki and Kanda, 2004;Ramamurthy and Bou-Zeid, 2014). Understanding the variability of modeled ET for each LCZ can inform urban planners aiming to integrate this information into the planning process. We used the interquartile range (IQR) as an indicator of variability (Fig. 5), with higher IQR indicating more variability. Predominantly vegetated LCZs are fairly homogeneous in modeled ET. For both annual and summer ET, "low plants" and "dense trees" showed the lowest IQR. "Bare soil or sand" was the most variable LCZ both annually and in summertime. Bare soil is likely converted to another land cover (vegetated or impervious) more rapidly than other LCZs, causing a mismatch with the LCZ classification from 2017 . To illustrate this point, we analyzed the variability of ISF and NDVI on a summer date within LCZs (Supplement, Fig. 13). For both ISF and NDVI, "Bare soil or sand" showed the highest variability, while natural LCZs such as "dense trees" showed the lowest variability, similarly to modeled ET. The second-highest IQR for modeled ET occurred for "sparsely built" and "open lowrise" annually and in summertime, respectively. The high heterogeneity of modeled ET in LCZs such as "sparsely built" can be explained by their mixture of vegetated and impervious surfaces (Stewart and Oke, 2012). Some variability may also be attributed to LCZ misclassification (Demuzere et al., 2021;Fenner et al., 2017). However, validating the LCZ classification is not the aim of this study, which only uses LCZs to interpret modeled ET. Several LCZs have been identified to be quite similar in the arrangement of buildings, building height, surface cover, and thermal inertia, such as "Open lowrise" and "Sparsely built" (Bechtel et al., 2017a). Nevertheless, we chose not to merge similar LCZs in this analysis to enable direct comparison to other LCZ-based studies and as many highly similar LCZs do not occur in Berlin .

Discussion
Mapping ET is highly relevant to making cities climate-resilient as more of the world population is exposed to extreme heat, flooding, and drought (Kundzewicz et al., 2014;Li and Bou-Zeid, 2013;Mora et al., 2017). In this study, we demonstrated the capacity of an open-source data-driven approach to accurately map urban ET at a high spatial and temporal resolution.

Advantages of this method
Our modeling approach showed a higher reliability and accuracy than the majority of studies modeling hourly urban ET. In most cases, urban ET modeling approaches can only be validated by flux tower measurements at one to two points in the city Karsisto et al., 2016). Our ET maps were validated at two sites, with the more vegetated site showing a R 2 of 0.76 and a RMSE of 0.0289 mm/h, and the dense city-center site showing a R 2 of 0.56 and a RMSE of 0.0171 mm/h. For urban land surface models, R 2 ranges from 0.017 to 0.79 and RMSE ranges from 0.0295 to 0.1143 mm/h (Järvi et al., 2011;Karsisto et al., 2016;Rafael et al., 2020;Ward et al., 2016). A micrometeorological approach to mapping LE in three European cities showed R 2 ranging from 0.07 to 0.42 and RMSE ranging from 0.0268 to 0.0450 mm/h . As in our study, higher R 2 and RMSE values have been reported for greener sites compared to dense inner-city sites (Rafael et al., 2020;Rocha et al., 2022b;Ward et al., 2016), which is likely due to the higher variation of ET at more vegetated sites.
Furthermore, the urban ET dataset produced in this study had a higher spatial (10-m) and temporal (hourly) resolution than nearly all studies modeling ET at the city-scale. To plan sustainable and climateresilient cities, it is necessary to characterize ET at a high spatiotemporal resolution for the entire city extent (Milojevic-Dupont and Creutzig, 2021; Rocha et al., 2022a;Saher et al., 2021). Many of the studies estimating urban ET only produced maps for a limited area rather than city-wide Gillefalk et al., 2022;Nouri et al., 2020), or not at all (e.g., point-scale) (Gillefalk et al., 2021;Karsisto et al., 2016;Rafael et al., 2020;Rocha et al., 2022b;Vulova et al., 2021;Ward et al., 2016). Wang et al. (2016), the only other study to use a data-driven approach to map urban ET, modeled ET at a 500-m spatial resolution every 8 days. Cong et al. (2017) modeled annual ET at a 500-m resolution using an urban SEB model. Rocha et al. (2022a) uses a Soil-Vegetation-Atmosphere Transfer (SVAT) model to produce hourly ET maps at a 1 km resolution, which are then corrected based on the assumption that impervious surfaces do not evaporate. Recently, Landsat imagery has been used to model urban ET at a higher spatial resolution (30-100 m), but the temporal resolution of these studies remains limited by Landsat's long revisiting cycle (16-day) (Saher et al., 2021;Zhou et al., 2019). For instance, Faridatul et al. (2020) (Abunnasr et al., 2022). In Basel, Switzerland, daily ET was modeled at a 100-m spatial resolution for 22 Landsat scenes using a micrometeorological approach . Using Sentinel-2 or harmonized Landsat and Sentinel-2 time series (Claverie et al., 2018) can enhance the spatial and temporal resolution of urban ET products, as demonstrated in our study, which will be beneficial for urban planning applications.

Features of the proposed method
ET maps can support the management of urban green spaces to mitigate heat risk for urban dwellers. An established trade-off exists between the capacity of UGI to provide ecosystem services, in particular evaporative cooling, and the irrigation needed to maintain this cooling capacity (Cuthbert et al., 2022;McCarthy and Pataki, 2010;Nouri et al., 2019). As urban areas are increasingly faced with the double burden of heat waves and droughts, climate change adaptation with UGI may exacerbate water scarcity (Cuthbert et al., 2022;Miller et al., 2020). This trade-off can be visualized on a map showing the ET anomaly between the hottest day in 2019 and the average daily summertime ET (Fig. 6), with the increase in ET varying with the vegetation health, type, and water availability. In Tiergarten, which is dominated by tree vegetation, ET is ~30% higher than average, enhancing the microclimate where heat risk is highest (Fig. 6b) (Dugord et al., 2014). In Tempelhof Field, a large open green space, the increase in ET is more differentiated, being higher (30-50%) in the south of the meadow but lower in the middle (~20%) (Fig. 6c). The Karl-Marx-Allee (KMA), an inner-city neighborhood (above the river, Fig. 6d), is a hotspot for thermal stress, while housing a high proportion of elderly residents (Knaus and Haase, 2020). Notably, the ET anomaly in the KMA area (up to 100%) is even higher than in conventional green spaces, likely due to irrigation maintaining the vitality of the vegetation. Fig. 6e shows a mixture of agricultural fields and residential buildings, with an industrial area showing the highest ET anomaly (>60%). Our mapping approach reveals how the cooling capacity of green spaces responds to drought stress, which can better inform UGI implementation for climate mitigation.
Our proposed approach can capture the impact of extreme events such as heat waves on ET, as demonstrated with the example of 26 June 2019 (Fig. 6), a capacity missing from existing urban ET models due to constraints in temporal resolution. Due to the long repeat cycle (16-day) of Landsat thermal infrared (TIR) imagery, LST is often not captured during extreme weather events (Zhou et al., 2019), when evaporative cooling is especially important. At the same time, ground-based observations cannot reflect the spatial distribution of urban ET, which is critical for urban planning applications (Nouri et al., 2013b;Saher et al., 2021). Therefore, our study fills a significant gap in the field of urban ET by developing a method which can supply city-wide ET maps for the specific time step when this information is most urgently needed, such as extreme heat waves.
Grasslands and shrublands, which showed the highest ET in our study, may provide a critical evaporative cooling service. On average, grasslands and shrubs converted more of annual precipitation to ET (72% and 68%) than dense trees (61%). Our study also provided insights on the variation of ET between urban LCZs. Among the urban-specific LCZs, "sparsely built" showed the highest ET, followed by "open lowrise" and "open highrise", both annually and in summertime. These findings are in line with previous studies, which found that lower development intensity is associated with higher annual ET (Liu et al., 2010;Wan et al., 2021). "Large lowrise", a mostly paved LCZ, showed relatively higher ET rates in wintertime than more vegetated LCZs (Supplement; Fig. 9), which can be attributed to higher rates of evaporation from interception loss from impervious surfaces (Miralles et al., 2020;Ramamurthy and Bou-Zeid, 2014). Artificial surfaces have been found to contribute substantially to ET in urban areas, especially in wintertime (Moriwaki and Kanda, 2004), as our study confirms. Our results emphasize the need to map ET at different temporal scales, rather than only annually, to better characterize its influence on heat mitigation and the urban hydrological balance.
One advantage of a data-driven mapping approach is its capacity to uncover patterns not properly characterized by physical models (Reichstein et al., 2019), including the rooting depth, which may be more developed in grasslands than previously assumed (Williams et al., 2012) or the higher interception loss over paved surfaces in wintertime. Another benefit of using a data-driven approach is the flexible integration of diverse and novel datasets (Kuhn and Johnson, 2013). For instance, urban microclimate can be integrated into ET mapping in future work using crowdsourced weather data (Venter et al., 2020;Vulova et al., 2020). The presented methodology can be also used to map other fluxes measured by eddy covariance, such as CO 2 (Kotthaus and Grimmond, 2014), at a high spatio-temporal resolution by combining remote sensing and GIS data, flux footprint modeling, and ML.

Limitations and future directions of research
There are a few limitations associated with the presented data-driven method. Our approach therefore allows transferability to similar land cover and meteorological conditions to those used in training. Retraining with local flux data may be necessary to achieve comparably high accuracy, especially in other climate zones. Naïve extrapolation is a well-known risk in machine learning (Meyer and Pebesma, 2022), but process-based models also rely on extensive calibration to local observations to optimize the accuracy that may jeopardize transferability (Abunnasr et al., 2022). As the developed method relies on globally available, open-access satellite data and standard meteorological observations, it can be used to predict ET in other cities. The increasing availability of open-access flux tower data from cities around the world will allow for ML models to be trained with local data from urban areas and for the transferability of models to other cities to be evaluated (Lipson et al., 2022). Though directly comparing our model output and performance to process-based ET models is beyond the scope of this study, evaluating multiple models estimating ET in urban areas against EC observations is an important nascent area of research (Best et al., 2015;Lipson et al., 2020). Although eddy covariance (EC) is one of the best-suited methods to directly measure urban ET (Nouri et al., 2013b), the noise of EC data may introduce inaccuracies into our ET maps. Currently, in situ ET data in cities is exceedingly scarce (Nouri et al., 2015;Saher et al., 2021). In this study, we relied on two flux towers for validation based on data availability. In a heterogeneous urban environment, it would be preferable to have a dense network of in situ sensors covering different surfaces for validation of a high-resolution ET product. Therefore, in future work, alternative sources of ground-truth data, such as lysimeter and soil moisture data (Nouri et al., 2013b), should be collected and used to independently validate ET models in urban areas. Collecting more in situ urban ET observations will allow for a better quantification of the modeling uncertainty over different LCZs. As some uncertainty is derived from footprint modeling (Neftel et al., 2008;Schmid, 2002;Vesala et al., 2008), future studies would benefit from comparing the influence of using different footprint models. Lastly, further investigation into which RET or potential ET (PET) equations are most appropriate for urban areas is needed (Shafieiyoun et al., 2020). However, none of the current PET and RET equations can fully capture the dynamics of ET in an urban environment, as impervious surfaces completely alter their assumptions, which are based on soil-vegetationatmosphere interactions (Allen et al., 2005;Priestley and Taylor, 1972).

Conclusion
Our study presented an innovative approach to mapping urban ET at a high spatial and temporal resolution for an entire city. In order to model the influence of heterogeneous urban land cover on ET, Sentinel-2 imagery and open geodata were integrated into Random Forest models using flux footprint modeling. Hourly ET maps for two years were validated at two flux towers in Berlin, Germany, showing higher accuracy than most studies modeling urban ET. Furthermore, urban ET was modeled at high spatial (10-m) and temporal (hourly) resolution, which very few studies have achieved for an entire city so far. The approach provides detailed urban ET maps, a crucial source of data to guide urban planning aiming to mitigate heat-related risks and manage scarce water resources. We illustrated the advantage of the high spatio-temporal resolution of our methodology by mapping the effect of an extreme heat wave on urban ET, when evaporative cooling is most critical for Fig. 6. Map of the anomaly between the daily evapotranspiration (ET) (mm/ day) of a heatwave date (26 June 2019) and the average daily ET in summertime in 2019. The ET anomaly is depicted for a) the entire city and zoomed extents around b) a large urban forest (Tiergarten) and c) the largest park, dominated by low vegetation (Tempelhof Field), d) a typical inner-city neighborhood with mid-rise buildings (Karl-Marx-Allee neighborhood), and e) an area with croplands and residential buildings. White colored areas represent water bodies, which were excluded from ET modeling. human health and well-being. To provide further insights for urban planners, hydrologists, and climatologists, we aggregated the annual and summer ET maps for Berlin using Local Climate Zones (LCZs). The LCZ analysis confirmed that modeled ET varied spatially as expected, with fully vegetated LCZs such as grasslands and shrublands showing the highest ET rates while built-up areas were associated with lower ET rates. Our study mapped urban ET using machine learning, which allowed for patterns not currently well-characterized in physical models to be captured, including the higher interception loss over impervious surfaces in wintertime.
We illustrated the potential of an open-source data-driven approach to accurately map urban ET in a large mid-latitude city. In future studies, the transferability of this methodology to other cities should be tested using the increasingly available flux tower data from cities worldwide. Considering the current scarcity of ground-truth ET data in cities, future work should use alternative sources of ET data, such as lysimeters, to further quantify the modeling uncertainty over different land covers. To facilitate transferability, we rely on open-access remote sensing and meteorological data and provide our code publicly in a GitHub repository (https://github.com/stenkavulova/Map_UrbanET). The methodology and ET maps presented in this study can support decisionmakers in identifying hotspots with high heat risk and targeting mitigation interventions, such as the incorporation of green infrastructure. To support technical users, the ET maps generated in this study are freely available on Zenodo at https://doi.org/10.5281/zenodo.7561125 (Vulova et al., 2023). Our approach can also be used to simulate the effects of climate change and land cover conversion scenarios on urban ET. High-resolution ET data can allow for the world's growing cities to mitigate health risks and water scarcity in the face of accelerating climate and global change.

Data and code statement
The code used for this study is publicly and permanently available on Zenodo at https://doi.org/10.5281/zenodo.6521095. The monthly and annual ET maps are provided open source on Zenodo at https://doi. org/10.5281/zenodo.7561125. The output data from footprint and machine learning modeling are available on request from the first author. The raw eddy covariance data belongs to the Urban Climate Observatory (UCO) Berlin and are available on request from Dr. Fred Meier (fred.meier@tu-berlin.de). The meteorological (DWD) data are freely available at https://opendata.dwd.de/climate_environment/ CDC/observations_germany/climate/ and can also be downloaded using our open code (https://doi.org/10.5281/zenodo.6521095). The source code to compute and download our NDVI time series using Google Earth Engine is also provided in the same Zenodo publication. Urban Atlas products (building height and Land Use and Land Cover (LULC)) are publicly available at https://land.copernicus.eu/local/ur ban-atlas. The code to derive ISF from Urban Atlas LULC is provided in the Zenodo publication. VH and VF are freely available from the Berlin Digital Environmental Atlas (https://fbinter.stadt-berlin.de/fb/ index.jsp). The LCZ map used in this study is freely available at https ://lcz-generator.rub.de/factsheets/6d7f501c212dc888e32a3a9a5740 ce930addd3bc/6d7f501c212dc888e32a3a9a5740ce930addd 3bc_factsheet.html.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The code and data associated with this study are publicly and permanently available on Zenodo. Please see our "Data and code statement" for more information.