Urban Climate 39 (2021) 100982 Contents lists available at ScienceDirect Urban Climate journal homepage: www.elsevier.com/locate/uclim Sensitivity study of PBL schemes and soil initialization using the WRF-BEP-BEM model over a Mediterranean coastal city Ricard Segura a, Alba Badia a, Sergi Ventura a, Joan Gilabert a, b, c, Alberto Martilli d, Gara Villalba a, e, * a Sostenipra Research Group (SGR 01412), Institute of Environmental Sciences and Technology (MDM-2015-0552), Z Building, Universitat Autònoma de Barcelona (UAB), Campus UAB, 08193 Bellaterra, Barcelona, Spain b GAMA Team Department of Applied Physics, University of Barcelona (UB), Barcelona, Spain c PCOT Team, Institute Cartographic and Geological of Catalonia (ICGC), Barcelona, Spain d Research Center for Energy, Environment and Technology, CIEMAT, Madrid, Spain e Department of Chemical, Biological and Environmental Engineering, Universitat Autònoma de Barcelona (UAB), Campus UAB, 08193 Bellaterra, Barcelona, Spain A R T I C L E I N F O A B S T R A C T Keywords: Due to increased urbanization and global warming, cities are experiencing more heat wave (HW) Heat wave events that cause extreme heat stress. To mitigate such effects, a better understanding of the Urban parameterization impact of urban morphology on the boundary layer development is needed. This study in­ WRF BEP BEM vestigates the sensitivity of mesoscale simulations using the WRF model coupled with the building Boundary layer Urban morphology effect parameterization and the building energy model (BEP-BEM) at a 1-km resolution to 1) soil moisture initializations; 2) the inclusion of site-specific urban morphology parameters; and 3) the planetary boundary layer (PBL) scheme. A HW episode that occurred in the metropolitan area of Barcelona serves as the case study. We find that the use of a high-resolution land data assimilation system (HRLDAS) to initialize soil properties results in larger temperature diurnal range, but it did not improve the performance of simulated temperatures compared to using low-resolution ERA5 data. The inclusion of site-specific urban parameters improved the representation of urban fractions, reducing the night-time overprediction of 2-m temperatures compared to using default urban parameters. Overall, the Bougeault-Lacarrere (BouLac) scheme represents the PBL- height noontime observations better than the Mellor-Yamada-Janjic (MYJ) scheme. This was related to a better representation of daytime near-surface temperatures by the BouLac scheme compared to the MYJ scheme 1. Introduction More than half of the world's population is currently living in urban areas, and this demographic trend is expected to increase up to 68% in 2050 (83.7% in the case of Europe), according to the United Nations (UN, 2019). With the dramatic global rise in urbanization, the impacts of climate change on urban areas have become a major worldwide concern. Some of these impacts that have already been observed or estimated include increasing temperatures (Hunt et al., 2013; IPCC, 2014), changing precipitation dynamics (IPCC, 2014), * Corresponding author at: Sostenipra Research Group (SGR 01412), Institute of Environmental Sciences and Technology (MDM-2015-0552), Z Building, Universitat Autònoma de Barcelona (UAB), Campus UAB, 08193 Bellaterra, Barcelona, Spain. E-mail address: gara.villalba@uab.cat (G. Villalba). https://doi.org/10.1016/j.uclim.2021.100982 Received 10 April 2021; Received in revised form 4 August 2021; Accepted 10 September 2021 Available online 15 September 2021 2212-0955/© 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). R. Segura et al. Urban Climate 39 (2021) 100982 increasing frequencies and intensities of extreme events such as heat waves (HWs) (Stott et al., 2004; Hunt and Watkiss, 2011; IPCC, 2014), worsening air quality (Jacob and Winner, 2009), and accelerating sea-level rise (Hunt and Watkiss, 2011; IPCC, 2014). The Mediterranean region has been identified as a highly vulnerable area to the effects of climate change (Giorgi, 2006) and is expected to see an increase in the frequency of longer-lasting and more severe HWs (Meehl and Tebaldi, 2004; Ballester et al., 2009). Historically, the most-studied urban micro-climatic feature is the urban heat island (UHI) phenomenon (Oke, 1982; Moreno-Garcia, 1994; Arnfield, 2003), which refers to the tendency of urban areas to be warmer than their peri-urban and rural surroundings, particularly at night. This effect arises due to multiple causes, such as the excessive absorption of heat by building materials, the geometries of urban canopies, the decreased evapotranspiration owing to the lack of vegetation, and the heat released from anthro­ pogenic emissions (Oke, 1982; Grimmond, 2007). Moreover, global warming is expected to exacerbate the differences between urban and rural temperatures (McCarthy et al., 2010; IPCC, 2014), making urban dwellers more vulnerable to heat stress and other heat- related pathologies (Ingole et al., 2020). Li and Bou-Zeid (2013) revealed that HWs interact non-linearly with UHIs to produce extremely high heat stress on urban citizens. The combined effect of UHIs and climate change also implies a significant challenge for the sustainability of energy and water systems in urban environments and the integrity of urban infrastructures and ecosystems (Hunt and Watkiss, 2011). Responding with urgency to these challenges requires a deep understanding of the behaviour of the urban boundary layer (UBL), particularly during extreme heat events. In this context, numerical models play a key role, both as scientific tools to help researchers understand how different physical mechanisms shape UBL evolution, and as practical tools to aid urban planners and policy makers in evaluating the efficiencies of different mitigation strategies in the context of a changing climate (Chen et al., 2011). However, today's atmospheric models have a range of different options in terms of their physical parameterizations, initializations, etc., that often need to be adapted to each case study. The aim of this paper is to test the sensitivity of the mesoscale meteorological Weather Research and Forecasting (WRF) model (Skamarock et al., 2019) with the multilayer urban canopy parameterization (UCP) Building Effect Parameterization and the Building Energy Model (BEP-BEM) (Salamanca and Martilli, 2010) to three configurations that we identified as the most relevant, using the coastal Mediterranean city of Barcelona and its metropolitan area as a case of study. These configurations are described as follows: 1) Soil moisture and temperature initialization. Soil moisture is a key parameter that affects both the heat capacity of the ground and the partitioning between the sensible and latent heat fluxes, while soil temperature affects the amount of heat transferred between the surface and inner soil layers. The initialization of these variables in land surface models (LSMs) represents a challenge in atmo­ spheric modelling due to sparse available measurements (Best and Grimmond, 2014). Global reanalysis datasets, such as ERA5 (C3S, 2017), derive these variables at low resolutions (31 km × 31 km) that cannot resolve the complex heterogeneities of land cover types and soil textures at the urban scale. To minimize the uncertainties related to the initial conditions in WRF simulations, spin-up strategies have been developed, such as the uncoupled LSM high-resolution land data assimilation system (HRLDAS) (Chen et al., 2007), which can be used to simulate seasonal variability in soil properties at a high resolution. The HRLDAS system provides high-resolution soil moisture and temperature fields for WRF initialization, as was reported in Hong et al. (2011) for the simulation of three regions in the United States. 2) Urban morphology parameters. WRF-BEP-BEM requires the initialization of multiple urban canopy parameters that define the geometrical, radiative, energetic, and metabolic properties of urban areas. In this sense, the World Urban Database and Access Portal Tools (WUDAPT) project focused on creating a detailed database of urban description land-use/land-cover (LULC) maps for each city using the local climate zones (LCZ) classification (Stewart and Oke, 2012; Brousse et al., 2016). The LCZ establishes a framework for classifying neighbourhoods within the same city with similar urban cover types, structures, materials, and human activities, thus showing similar tendencies regarding UHI effects and establishing a range of values for the urban parameters of each LCZ that can serve as inputs for the UCP (Stewart and Oke, 2012). Alternatively, the WRF-BEP-BEM system has also incorporated the option to directly initialize urban canopy parameters into WRF grids in a point-by-point based approach (Chen et al., 2011). In this way, each grid has its own unique set of urban morphological parameters, such as canopy geometry and urban fraction (λu, defined as the fraction of the impervious surface area to the total plan area). Hammerberg et al. (2018) compared the WUDAPT methodology with the LCZ-averaged urban morphology to initialize WRF BEP-BEM over the city of Vienna (Austria). Their results showed that using LCZ-averaged urban morphology from specific data of the city of Vienna only provided marginal improvements compared to the WUDAPT methodology due to the high homogeneity of the city and the high horizontal resolution of their WRF simulations (500 m). Additionally, Wong et al. (2019) compared the WUDAPT and the point-by-point methodologies for defining UCP datasets for the city of Hong Kong. Their main findings were that the use of tables with default values instead of local specific data, such as those used in Stewart and Oke (2012), can add a large source of uncertainty in urban WRF-BEP-BEM simulations. 3) Planetary Boundary Layer schemes. In the WRF model, PBL schemes represent the turbulent sub-grid-scale processes by computing the boundary layer fluxes of heat, moisture, and momentum in the lowest part of the troposphere along with resolving vertical diffusion due to turbulence for the entire vertical column. Two available PBL schemes are available in WRF that can be coupled with the BEP-BEM urban canopy model: the Bougeault-Lacarrere (BouLac) (Bougeault and Lacarrere, 1989) and the Mellor- Yamada-Janjic (MYJ) (Janjic, 2002) schemes. Other studies have analysed the performance of the WRF model using these PBL schemes for the city of Barcelona, such as Banks et al. (2015) and Banks and Baldasano (2016). Banks et al. (2015), showed that the MYJ scheme was the least correlated to LiDAR estimates of the PBL height. Banks and Baldasano (2016) showed that the BouLac scheme estimated air and dew point temperature better than MYJ for the cities of Barcelona, Girona and Reus. However, these studies did not use urban canopy parametrization to represent urban atmospheric processes, and they relied on the BULK scheme (Liu et al., 2006), which does not consider urban heterogeneity. 2 R. Segura et al. 3 Fig. 1. (a) Domains of the WRF model for this study with 9, 3 and 1 km grid size for d01, d02 and d03, respectively. (b) Detail of the domain 3 with a terrain height map and the political boundaries of the AMB and Barcelona city (black contours). The red line represents the cross-section shown in Section 3.4. The red crosses show the position of the meteorological stations used in this study for the performance evaluation of the WRF simulations. The white dot shows the radiosonde launching position. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Urban Climate 39 (2021) 100982 R. Segura et al. Urban Climate 39 (2021) 100982 Recent studies have used numerical modelling approaches to investigate the UHI phenomenon of different cities (Arnfield, 2003; Salamanca et al., 2012; Li and Bou-Zeid, 2013; Sharma et al., 2014, 2016; Li et al., 2019). Some observational studies have found a spatial correlation between near-surface temperatures and λu (Kottmeier et al., 2007; Zhang et al., 2010; Ryu and Baik., 2012; Schatz and Kucharik, 2014; Li et al., 2019). Usually, near-surface temperatures increase with λu, presenting a linear trend (Schatz and Kucharik, 2014, 2015; Li et al., 2019). Li et al. (2019), used a linear regression between simulated two-meter temperatures (T2m) and λu using WRF coupled to an UCP for the city of Berlin, and a method for determining the UHI intensity of Berlin based on a set of validated λu linear functions. This comprehensive sensitivity analysis can constitute a reference for future urban modelling exercises by providing the most suitable WRF-BEP-BEM model configuration for a coastal Mediterranean city during a HW event. Moreover, this analysis constitutes a basis to improve our scientific understanding of the physical mechanism responsible for the development of the UBLs of these types of cities, as well as to derive information that can be used to design mitigation strategies. As an example, in the second part of the paper, the reference model configuration is used to investigate the impact of the studied urban area on UBL development and the impact of vegetation on the urban air temperature. This paper is organized as follows: in Section 2, the study region and case study are described, along with the WRF-BEP-BEM setup and configurations and the observation datasets used to evaluate the model performance; Section 3 describes the results and discussion of the various simulations; and finally, in Section 4, some conclusions are remarked. 2. Methods and data 2.1. Region of interest and case study In this study, we analyse urban atmospheric simulations of the metropolitan area of Barcelona (AMB) and its surrounding areas. The AMB is situated in the northwest region of the Mediterranean basin and in the northeast region of the Iberian Peninsula (Fig. 1a). The AMB coast is oriented from southwest to northeast, and the AMB is surrounded by the Mediterranean Sea, the Catalan Coastal Range (near 500 m in altitude) at the NW and two rivers, the Llobregat River (South) and the Besòs River (North) (Fig. 1b). The AMB covers an area of 636 km2 with a population of more than 3 million inhabitants (the 6th largest metropolitan area in Europe). The city of Barcelona has 1.6 million inhabitants and covers a total area of 101 km2, with a population density of 160 individuals/ha. The Mediterranean climate is made up of dry summers with a regime of local breezes and warm temperatures (25–27 ◦ C) influenced by the thermoregulatory effect of the sea, which in turn creates an increasing heat sensation effect due to humidity (Lionello et al., 2006). During the summer months, stagnant synoptic patterns and the absence of synoptic events make mesoscale convective pro­ cesses predominate, especially on the coast, which is influenced by sea/land breezes (Pérez et al., 2006). The day and night transitions of the sea/land breezes together with the UHI effect cause a complex structure in the atmospheric boundary layer, modifying the temperature and the dispersion of contaminants (Soler et al., 2011). Sea breezes start between 8 and 9 UTC (Barcelona's longitude is 2 Degrees E, therefore UTC can be considered equivalent to Local Solar Time), reaching their maximum strength at approximately 12–14 UTC, with higher values recorded near the coast but spanning hundreds of kilometres of the land surface until 18–20 UTC when they gradually disappear (Soler et al., 2011). At night, the coastal regions are characterized by land-to-sea breezes, which have much lower forces than those of the sea breezes, increasing the influence of the possible effects of topography on the winds (Soler et al., 2011). Our study case corresponds to a HW episode that occurred between 4 and 6 July 2015 in southern and central Europe. The maximum daily temperatures surpassed the 98th percentile for three or more consecutive days at numerous meteorological stations, meeting the definition of a HW episode established by the Catalan Meteorological Service (SMC). During this period, the maximum daily temperatures reached 40 ◦ C in the interior of the AMB and 35 ◦ C on the coast (SMC, 2015). The night-time temperatures were also considerably high during the entire HW, exceeding 25 ◦ C for several consecutive days in some neighbourhoods (SMC, 2015). Maximum daily temperatures were considerably high during the previous weeks over the entire Iberian Peninsula and extended for almost a whole month (27 June 2015 to 22 July 2015) with no precipitation. Most of the days during the HW are classified with the Shallow cyclone or Undetermined pressure gradient and Anticyclonic western advection synoptic types according to the synoptic classification methodology described in Miró et al. (2020). These two circulation patterns are the most common types of situations for the summer days in the Iberian Peninsula, making this HW an excellent case study representative of the most frequent synoptic events leading to HWs. Further details of the synoptic situation during the HW can be seen in Section 1 of the Supplementary material and in Section 2 of the Supplementary material of Gilabert et al. (2021a). 2.2. WRF implementation 2.2.1. Model set-up The Weather Research and Forecasting (WRF) v4.0 model (Skamarock et al., 2019) is used to study the urban boundary layer evolution over the AMB. The WRF model is a mesoscale numerical weather prediction system that allows atmospheric simulations to be performed across a wide range of spatial scales (from regional to urban) due to its nesting capabilities and its potential to be coupled to a UCP (Chen et al., 2011). In this study, WRF is coupled with the BEP-BEM scheme (Martilli et al., 2002; Salamanca and Martilli, 2010) to better represent the physical surface energy exchanges inside the urban canopy, the drag effect on the wind and the urban spatial heterogeneities. The WRF configuration used for this study consists of 3 two-way nested domains with horizontal resolutions of 9 km, 3 km, and 1 km (Fig. 1a). The domains were selected to include all geographical features that have synoptic or regional in­ fluences on the meteorology of the AMB. The large domain covers the entire Iberian Peninsula and Northern Africa. The middle domain 4 R. Segura et al. Urban Climate 39 (2021) 100982 contains the entire Catalan territory and extends in the north above the eastern Pyrenees. The small domain covers Barcelona Province, which contains the AMB on its eastern coast (Fig. 1b). All simulations cover the period from 00 UTC on 24 June 2015 to 00 UTC on 10 July 2015, using the first 24 h as the spin-up period and covering the three days considered to be the HW (July 4, 5 and 6, 2015). Only the innermost domain is evaluated against the observed data. The model characteristics and experimental configurations, such as the domain resolution and physical parameteri­ zations, are shown in Table 1, which have been derived from previous sensitivity studies (Ribeiro et al., 2021; Gilabert et al., 2021a). 2.2.2. Land-use classification For the LULC classification, the WRF model was adapted following the indications of Martilli et al. (2016) and the WUDAPT project to use the detailed urban descriptions of urban LCZs (Brousse et al., 2016; Gilabert et al., 2021b). The LCZ map of the AMB and adjoining areas was created following the WUDAPT classification (Bechtel et al., 2015) at a 100-m resolution. Then, it was improved through the reclassification and integration of high-resolution LULC maps such as Urban Atlas 2012, the LULC AMB 2015 and LiDAR data (Gilabert et al., 2021b). The LCZ map for the AMB region and adjoining areas covers a rectangular extension (ROI, Fig. 2a), which was used as the region of study. For this study, 10 urban LCZs and the LCZ E (reformulated as asphalt) are used in combination with the Corine Land Cover 2018 classification (CLC) (Büttner et al., 2017; https://land.copernicus.eu/pan-european/corine-land-cover/ clc2018) at a 100-m resolution for the other areas as the background land-use types. The CLC is remapped from its original 44 classes into 16 classes using the MODIS IGBP land classification and four urban LCZs following the method of Pineda et al. (2004); these maps are then used by the WRF model. Both LULC classifications are interpolated at each model domain resolution using the most common land-use class in each cell. Additionally, the green vegetation fraction (GREENFRAC) was updated using the fraction of green vegetation cover (FCOVER) product of the Copernicus Global Land Service (Smets et al., 2019) with a 1-km resolution. The default FPAR MODIS data did not represent well the vegetation of the highly urbanized coast of the AMB due to the low resolution. After interpolation to the resolution of the innermost domain, the most common LULC class inside the ROI (without considering water) is the evergreen needleleaf forest class (23.2%) (Fig. 2b), representing the natural regions of Collserola (512 m a.s.l) north of Barcelona city and the Garraf Massif (657 m a.s.l.) to the west of the AMB, both of which are located in the Catalan Coastal Range. Urban areas represent 47% of the non-maritime area of the ROI. The largest LCZ is the large low-rise zone (19.6%), followed by the open low-rise zone (12.7%) and the compact mid-rise zone (9.6%), while the other urban LCZs constitute minorities (5.3% combined). 2.2.3. Urban canopy parameters Different approaches are used to define the WRF-BEP-BEM urban parameters. The radiative and thermal properties must be specified for each surface type in the urban canopy (roofs, building walls and ground/roads). The radiative properties of building roofs (albedo and emissivity) were computed using Landsat 5 TM images with 30 m resolution (USGS, 2019) containing roof and ground albedos following the methods of Liang (2001) and Landsat 8 OLI/TIRS images (Vermote et al., 2016) containing roof and ground emissivity values following the methods of Ndossi and Avdan (2016); then, both properties were averaged for each urban LCZ. Geometric and surface cover parameters such as the urban fraction, street width, building width and height distribution were established using the two methodologies described in Section 2.3.2. The rest of the urban parameters (thermal properties, surface roughness, room occupancy and AC systems specifications) were adjusted using different assumptions to be as representative as possible for each urban LCZ and can be seen in Section 3 of the Sup­ plementary material of Gilabert et al. (2021a). The anthropogenic heat released by buildings due to AC systems can have an important contribution to the surface energy balance, especially during HW events. For this reason, the AC option of the BEP-BEM scheme is activated for all urban LCZs (except LCZ E Asphalt) during the daytime hours (6–18 UTC or 8–20 local time, LT), with a target indoor temperature of 24.0 ◦ C and a comfort range of 1.5 ◦ C, based on the recommendations of Burroughs and Hansen (2011). Table 1 Model characteristics and experiment configurations. Resolution and initial conditions Horizontal resolution 9 km × 9 km; 3 km × 3 km; 1 km × 1 km Domain dimensions 150 × 145; 118 × 118; 121 × 121 Vertical layers 57 (16 between the surface and 100 m) Top of the atmosphere 50 hPa Initial conditions ERA5 (C3S, 2017) with 31 km horizontal resolution, 137 vertical levels and 6-h separation Physics parameterizations Microphysics WRF Single-Moment 6-class scheme (Hong et al., 2006) Shortwave and longwave RRTMG scheme (Iacono et al., 2008) radiation Cumulus Kain-Fritsch scheme (Kain and Kain, 2004) (only the outermost domain) Surface / UCP Noah Land Surface Model (pervious areas) (Chen and Dudhia, 2001) / BEP-BEM (impervious areas) (Martilli et al., 2002; Salamanca and Martilli, 2010) Surface layer Monin-Obukhov Eta similarity scheme with Zilitinkevich thermal roughness length 5 R. Segura et al. 6 Fig. 2. a) Land-use classification for the innermost domain d03, from the interpolation of LCZ and CLC with the red box defining the borders of the ROI; b) Histogram of land-use classes inside the ROI (without considering water). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Urban Climate 39 (2021) 100982 R. Segura et al. Urban Climate 39 (2021) 100982 2.3. Sensitivity studies Various simulations were conducted, as summarized in Table 2, to perform the various sensitivity studies that are described in the next three sections. 2.3.1. High-resolution land data assimilation system implementation To properly initialize the soil temperature and moisture variables in rural grid cells and the pervious part of urban grid cells (1 - λu), LCZ-adapted HRLDAS v3.7.1 (Chen et al., 2007) was used. These variables were input to the Noah LSM in the WRF system. In this study, a comparison between a simulation initialized with low-resolution 1-day spin-up soil variables from the ERA5 reanalysis (BL- UM) data and a simulation initialized with high-resolution long spin-up soil variables obtained from HRLDAS outputs (BL-UM-HRL) is presented (see Table 2 for simulation descriptions). The BL-UM simulation was chosen as the reference simulation in this study because it displays the best performance for near-surface temperatures in comparison with the observations (see Section 3). HRLDAS was run independently for each nested domain using analysed meteorological forcing data and high-resolution static data. ERA5 reanalysis hourly 0.25◦ x 0.25◦ data were used to drive HRLDAS during a period of almost one year (1 July 2014–24 June 2015). According to Chen et al. (2007), soil properties require approximately 1 year to reach a quasi-equilibrium state using the HRLDAS model. Since soil variables are only used in the WRF model by the land surface model and not by the UCP, HRLDAS is used without coupling any urban scheme. In this way, HRLDAS considers only the pervious fraction of each grid cell (the vegetated and bare soil surfaces inside the grid cell), and the impervious fraction (all built surfaces including buildings) is ignored. Fig. 3 shows a comparison in the WRF model initialization of the temperature and moisture conditions of first soil layer (0 to 0.1 m depths) using ERA5 data and the HRLDAS outputs (after running the HRLDAS model for one year). These data were used to initialize the WRF model at 00 UTC on 24 June 2015; therefore, we compared the input files for the WRF model after running the pre-processor tool (real.exe). In the first soil layer, the ERA5 data had slightly warmer soil temperatures (up to 0.4 ◦ C warmer) than the HRLDAS data in low-altitude regions (Fig. 3c), while colder temperatures were seen in the ERA5 data in mountainous regions (down to − 1.3 ◦ C lower). In the case of the HRLDAS data, the soil temperature values showed high spatial uniformity, with first-layer temperatures ranging between 21.2 and 24.0 ◦ C. The fact that HRLDAS only considers the vegetated and bare soil surfaces of each grid cell reduces the spatial heterogeneity resulting from the LULC variability; thus, the spatial differences depend more on other factors such as the vegetation type and fraction and orography. In Fig. S2 of the Supplementary material, the temporal series of the HRLDAS output for the soil temperature variables show that the soil temperature difference between rural and urban areas barely changes over the one-year simulation. Fig. 4 shows the WRF model soil initialization variables of the second soil layer (0.1 to 0.4 m depths), which, together with the first soil layer, constitute the layers that have the greatest influence on the temporal evolution of near-surface variables in short simulations. For the second soil layer, the HRLDAS data generally had much higher soil temperatures (up to 4.6 ◦ C higher for the second layer) and higher spatial variability than the ERA5 data (see Fig. 4c). Comparing soil moisture, the WRF pre-processor only applies a simple bilinear interpolation for this variable in the case of the low- resolution ERA5 data, with low moisture values ranging from 0.17 m3 m− 3 (0.18 m3 m− 3) in inland areas to 0.04 m3 m− 3 (0.04 m3 m− 3) in coastal areas for the first (second) soil layer (Fig. 3e and 4e). In contrast, HRLDAS generally showed higher soil moisture (0.19 m3 m− 3 and 0.21 m3 m− 3 on average for the first and second layers, respectively) than the ERA5 data, and the model was able to represent spatial heterogeneities in the studied domain. Higher differences between the two initialization methods are seen near the southwest coast of the AMB (Fig. 3f and 4f). Minimum differences between urban and rural areas can be seen in the temporal series of soil moisture from the HRLDAS output (Fig. S3). 2.3.2. Sensitivity to urban morphology parameters In this study, we evaluated the capability of LCZs to represent the urban morphology of the AMB compared to using specific urban morphology for the region. In the case of the WUDAPT methodology, the default urban morphology parameters typical for each urban LCZ class suggested by Stewart and Oke (2012) and Stewart et al. (2014) were used to set the parameters in the URBPARM.TBL file of the WRF-BEP-BEM model (See Table S2 in the supplementary material). The building heights distribution is obtained from Brousse et al. (2016), based on the Stewart et al. (2014) tables but slightly adapted to the city of Madrid (Spain). In contrast, for the point-by-point approach, the buildings heights were extracted from LiDAR captures of the ROI and a digital surface model with a resolution of 2 m (ICGC, 2019). Then, the Urban Multiscale Environmental Predictor (Lindberg et al., 2018) was used to generate the walls height. The histogram of building heights was generated following the methodology described by Table 2 Simulations design for the different WRF-BEP-BEM model sensitivities performed in this study. The reference simulation (BL-UM) is highlighted in bold. Simulation Soil initialization PBL Scheme Specific urban morphology Comment BL-UM ERA5 BouLac Yes – BL-UM-HRL HRLDAS BouLac Yes – BL-noUM ERA5 BouLac No – MYJ-UM ERA5 MYJ Yes – BL-crop ERA5 BouLac No Urban modified to cropland/natural vegetation mosaic 7 R. Segura et al. Urban Climate 39 (2021) 100982 Fig. 3. Comparison of WRF model initialization at 00 UTC on 24 June 2015 in the first soil layer (0–0.1 m depth) for a) soil temperature with HRLDAS; b) soil temperature with ERA5; c) difference in soil temperature (HRLDAS – ERA5); d) soil moisture with HRLDAS; e) soil moisture with ERA5; f) difference in soil moisture (HRLDAS – ERA5). The black hatched lines represent the urban areas classified with an urban LCZ. Fig. 4. Comparison of WRF model initialization at 00 UTC on 24 June 2015 in the second soil layer (0.1–0.4 m depth) for the a) soil temperature with HRLDAS; b) soil temperature with ERA5; c) difference in soil temperature (HRLDAS – ERA5); d) soil moisture with HRLDAS; e) soil moisture with ERA5; f) difference in soil moisture (HRLDAS – ERA5). The black hatched lines represent the urban areas classified with an urban LCZ. Pappaccogli et al. (2017). In the case of urban fractions, yearly maximum NDVI maps derived from Landsat 8 OLI/TIRS images at a 30 m resolution (USGS, 2019) were reclassified to differentiate urban areas from vegetated areas, as is explained in Gilabert et al. (2021a). Then, these data were interpolated specifically for each urban classified grid cell at a 1-km resolution by the WRF pre-processing tool (WPS) (Chen et al., 2011). In this work, both methodologies were compared to evaluate the model improvement in the refinement of the urban morphology description and the capability of the WUDAPT methodology to represent the heterogeneities in the climate of the city (see Section 2 in the Supplementary material). To do this, one simulation where BEP-BEM used the default morphological urban parameters derived by the WUDAPT methodology and based on the ranges proposed by Stewart and Oke (2012) (BL-noUM) was conducted, and another simulation where BEP-BEM used the grid-cell-specific morphological urban parameters extracted with GIS on a point-by-point-based approach was performed (BL-UM). We studied how the simulations BL-noUM and BL-UM differed depending on which urban morphology parameter were used. In the 8 R. Segura et al. Urban Climate 39 (2021) 100982 case of BL-noUM, urban morphology parameters were defined using the default values for each LCZ at a 100-m resolution, with a relatively high homogeneity of building morphology and low fraction of pervious areas. Therefore, when the land-use classification was extrapolated from the LCZ resolution to the WRF grid cells at a 1-km resolution, the model overestimated the λu values, since coarser resolutions tend to contain rural areas inside the grid cell apart from urban areas. This is reflected in Fig. S4, were a comparison of urban morphology parameters between BL-noUM and BL-UM simulations can be seen. BL-noUM showed higher average values for urban grid cells (0.82) than BL-UM (0.42) (see Fig. S4c), especially for LCZ 6 (open low-rise) areas. The difference between the urban fractions used in the two simulations was lower in the city centre, where denser urban areas are located. In contrast, the mean street widths for the WUDAPT methodology are calculated for each LCZ from the mean building height and the mean canyon aspect ratio described in the look-up tables, obtaining much lower values than the ones extracted from the point-by- point methodology. BL-UM showed higher mean street widths (39.1 m) than BL-noUM (24.3 m) (Fig. S4l), especially for LCZ 10 (heavy industry, +64.2 m), LCZ 3 (compact low-rise, +19.3 m) and LCZ 8 (large low-rise, +16.8 m) areas. The area-weighted mean building heights comparison show similar values between the two simulations (Fig. S4f), with an average of 9.3 m for BL-noUM and 10.0 m for BL-UM, respectively. This is related to the similarity of the average values for the 3 more frequent LCZ in the AMB (LCZ 2, 6 and 8) between the city of Barcelona and the values obtained from Brousse et al. (2016), adapted to the city of Madrid, which has a similar building height distribution than Barcelona. This is also seen on the mean building width (Fig. S4i), for which both simulations show similar average values of 19.6 m for BL-noUM and 19.6 m for BL-UM. 2.3.3. Sensitivity to the PBL scheme In this study, the sensitivity of the model to the PBL scheme was analysed by comparing the BouLac (Bougeault and Lacarrere, 1989) and MYJ (Janjic, 2002) schemes. For this study, two WRF-BEP-BEM simulations were compared using these two PBL schemes: MYJ-UM for the MYJ scheme and BL-UM for the BouLac scheme (Table 2). Both schemes are one-and-a-half-order prognostic turbulent kinetic energy (TKE) schemes with local closure, although the BouLac scheme has a non-local counter gradient term for convective conditions. The MYJ scheme, which is more widely used than the BouLac scheme in atmospheric models (Banks et al., 2015), is a modified version of the old ETA scheme from the MM5 model. The BouLac scheme has been designed for use with the BEP multi-layer model and has been tested most extensively with BEP and BEP-BEM (Salamanca et al., 2011, 2012; Sharma et al., 2014; Huang et al., 2019; Teixeira et al., 2019; Ribeiro et al., 2021). 2.3.4. Sensitivity of the UBL to urbanization To determine the impacts of the city of Barcelona and the surrounding metropolitan areas on the UBL and near-surface temper­ atures, along with sea and land breezes, an additional simulation (BL-crop) was performed and compared with the BL-UM simulation. For BL-crop, the same configuration as that of the BL-UM was used, but the land-use of the urban areas inside the ROI was changed to the cropland/natural vegetation mosaic class of the modified MODIS IGBP classification (Table 2). In this way, we changed the AMB landscape into a mix of croplands, shrublands and grasslands, without any type of vegetation occupying more than 60% of the total landscape, such as the current cropland areas near the Llobregat River in the south of the city. This class is used in the WRF-BEP-BEM model to define the characteristics of the natural pervious areas inside an urban classified grid cell. The vegetation fraction variable (SHDFAC in VEGPARM.TBL of the WRF model) and the leaf area index (LAI), were kept constant between the two simulations, with the only difference being that the λu value was decreased to 0 in each grid cell. For all the simu­ lations in this study, the vegetation fraction and the LAI of the pervious urban grid cell components are established with the definitions of the cropland/natural vegetation mosaic class of the VEGPARM.TBL: 0.80 and 4.29 m2 m− 2, respectively, which are values that do not coincide with the real values in the AMB for this period (SHDFAC~0.26 and LAI ~ 0.67 m2 m− 2, respectively) but are used to not tune the model with data that are not replicable for other study cases. 2.4. Observed data 2.4.1. Near-surface datasets The WRF simulations were evaluated using observational data from the SMC and the Spanish Meteorological Agency (AEMET). Urban and rural stations were used for this evaluation with a total of 19 stations (Fig. 1b), 14 of which are located inside the AMB and 5 of which are located outside. The characteristics of these stations are provided in Section 4 of the Supplementary material. The stations were classified into 5 groups according to their LULC at a 100-m resolution: 1) urban fabric and non-vegetated urban areas (UF&N­ VUA); 2) non-agricultural vegetated urban areas (NAVUA); 3) agricultural areas (AA); 4) forest and seminatural areas (F&SNA); and 5) water bodies (WB). The stations provide hourly point measurements that are compared with spatially averaged 1 km × 1 km hourly model output data, so we should expect more spatial variability in the observations than in the model outputs (Sharma et al., 2016). Moreover, a coarse 1 km × 1 km grid cell can contain multiple LULC types in the region it covers; therefore, misclassifications could occur between the LULC class of a given grid cell and the local surroundings of a station. For example, several stations from the NAVUA and AA groups (urban parks and agricultural lands) were placed on grid cell that were classified as urban LCZs at a 1 km × 1 km resolution due to the proximity of larger industrial and residential areas nearby. It is expected that for these stations, discrepancies would be found between the observations and model outputs due to the significant differences between the characteristics of the areas. These discrepancies are discussed during the model evaluation (Section 3.1). In the case of three stations, where the measurement equipment is placed on rooftops at 40–50 m above ground, the temperature and wind speed observations were compared with the outputs of the nearest vertical level in the WRF simulation, while for the other stations, the observations were compared with diagnostic 2-m temperature 9 R. Segura et al. Urban Climate 39 (2021) 100982 (T2m) and 10-m wind speed (WS10m) from the model outputs. 2.4.2. Vertical profiles The surface evaluation was complemented with an assessment of the vertical structure of potential temperature and wind speed using radio-sounding measurements from radiosonde launches performed by the SMC. The launches were made in the city of Barcelona (41.384495 N, 2.117497E, 98 m a.s.l.) every day at 12 UTC. The radiosonde equipment measured temperature (◦ C), relative humidity (%), wind speed (m s− 1) and direction (◦ ), barometric pressure (hPa) and height (m). From these radio-soundings, the planetary boundary layer height (PBLH) could be inferred for comparison with the simulated values. To estimate the PBLH from the radio- soundings, we use the definition of Stull (2000) and Seibert et al. (2000), which defines the top of the PBL as the minimum height where the virtual potential temperature (θv) is equal to that of the surface. This methodology is also used in the BouLac scheme to diagnose the PBLH. Note that the lowest height measured by the radiosonde is at the rooftop of a building (33 m from the ground), while the first level measured by the WRF model is around 2.5 m above the ground. 3. Results and discussion Statistical analyses of the root-mean-square error (RMSE), mean bias (MB) and correlation factor (R) were used to compare the model simulations with the observations. The equations used for computing these indicators and the methodology to average them between stations is included in Section 4 of the Supplementary material. The standard deviation is included when the differences between averaged fields or averages of the statistical indexes are presented. The simulations of each sensitivity test (BL-UM-HRL, BL-noUM and MYJ-UM), further described in the following subsections, are tested against the reference simulation (BL-UM) with a two-tailed Welch's t-test for each variable, statistical index and LULC group. The Fig. 5. Temperature analysis for the performance of WRF simulations. Root mean square error (RMSE) (top), Mean bias (MB) (middle) and cor­ relation coefficient (R) (bottom) at daytime (left) and night-time (right) for each simulation and LULC group of stations. 10 R. Segura et al. Urban Climate 39 (2021) 100982 differences between modelled results and observation data are assessed through two-sample t-tests with two-tailed distribution for each pair of simulations, where we reject the null hypothesis with a 95% confidence interval if the p-value is lower than 0.05 (please see Section 4 of the Supplementary material for a more detailed description). We only remark in the next sections the differences between statistical indexes of each simulation that are statistically significant. 3.1. Evaluation of the reference simulation The near-surface temperature analysis is displayed in Fig. 5 for each LULC group for the period of study (25 June to 10 July 2015) over the daytime (6–18 UTC) and night-time (19–5 UTC) periods. During the daytime, BL-UM showed an average RMSE of 1.61 ± 0.39 ◦ C. The lowest RMSE for this period was obtained for the AA stations (1.38 ± 0.29 ◦ C) corresponding to overpredictions of near-surface temperatures. The highest RMSE was seen for the F&SNA stations (1.79 ± 0.56 ◦ C). Some reasons for this underestimation of T2m in F&SNA stations could be: 1) an overestimation of the Bowen ratio, defined as the ratio between the sensible and latent heat fluxes, and 2) an overestimation of wind speed, that causes the reduction of air temperatures due the advection of cooler air from sea-breezes. The BL-UM simulation shows a high correlation with the average daytime temperature observations (0.91 ± 0.03), with higher correlations observed for AA stations (0.93 ± 0.02). At night, BL-UM showed a reduction in the average RMSE (1.43 ± 0.56 ◦ C) with respect to the daytime period, although the error increased considerably for the AA (2.31 ± 0.71 ◦ C) stations, corresponding to overestimations of T2m. The model also overestimated night-time T2m for NAVUA stations (MB = 1.03 ◦ C) and the other stations in lower measure (MB < 0.4 ◦ C). The overpredictions of the NAVUA and AA stations could be explained since most of these stations were classified as urban LCZs (see Table S4), although other rural stations also showed positive MBs. Stations with higher night-time RMSEs corresponded to stations that are located on agri­ cultural lands or in urban parks but are classified as urban LCZs. LULC misclassifications can have a considerable impact on the Fig. 6. Wind speed analysis for the performance of WRF simulations. Root mean square error (RMSE) (top), Mean bias (MB) (middle) and cor­ relation coefficient (R) (bottom) at daytime (left) and night-time (right) for each simulation and LULC group of stations. 11 R. Segura et al. Urban Climate 39 (2021) 100982 statistical analysis of results when compared to observations (Ribeiro et al., 2021). Overall, the simulated night-time temperatures were less correlated with the observations (0.86 ± 0.06) than the simulated daytime temperatures, except for at the F&SNA stations (R = 0.90 ± 0.02). Similar to Fig. 5, Fig. 6 displays a statistical analysis of the modelled near-surface wind speed compared with the observed data. The statistical analysis of the wind speed showed an average RMSE value of 1.67 ± 0.54 m s− 1 and 1.21 ± 0.48 m s− 1 for daytime and night- time, respectively, with lower RMSE values obtained for the NAVUA stations (0.82 ± 0.02 m s− 1, 0.65 ± 0.02 m s− 1). The analysis revealed a significant variability of the MB (− 0.21 ± 1.34 m s− 1, − 0.39 ± 1.36 m s− 1 and 0.45 ± 1.18 m s− 1 at daytime; − 0.10 ± 0.75 m s− 1, 0.10 ± 0.60 m s− 1 and -0.08 ± 1.35 m s− 1 at night-time for UF&NVUA, AA and F&SNA stations, respectively), which we mostly attribute to poor simulation of wind flow over these stations. For some stations, the wind speed was clearly underestimated because the model considered higher building density than the actual local surrounding of these stations. The presence of dense building exerts a high drag effect on the wind which reduces its strength near the surface. For other stations placed on rooftops and two rural stations that measure 2-m wind speed, wind speed was overestimated because the model was not able to reproduce the increased surface drag of building rooftops or the ground. Our RMSE results (1.57 ± 0.38 ◦ C and 1.49 ± 0.49 m s− 1) for temperature and wind speed, respectively) and correlation factors (0.94 ± 0.03 and 0.73 ± 0.10) for the entire study period agree with the results of the study of Ribeiro et al. (2021), who also analysed meteorology in the AMB. Their study also uses the WRF-BEP-BEM model, with no specific urban morphology, to study the weather of the AMB for an entire summer month (July 2016). We obtain similar RMSE values and correlation factors for near-surface temper­ atures and wind speed than those obtained in their study (1.9 ◦ C and 0.91 for temperature and 1.5 m s− 1 and 0.73 for wind speed), despite the differences between the two studies. However, they observed an underprediction of daytime near-surface temperatures, while in our case, we observed an overprediction of night-time temperatures. These differences may be related to the differences between the two case studies and to the definitions of the urban canopy parameters. Fig. S6 shows a comparison of the time series of simulated near-surface temperatures by the BL-UM simulation and the observed measurements averaged among all the stations. In general, the model can capture the temporal evolution of near-surface temperatures. The BL-UM simulation underpredicts maximum daytime near-surface temperatures for the 1 and 4 of July and overpredicts the temperatures considerably the day after the HW (July 7), considering this day to be the warmest day of the period of study. Day 7th of July was the hottest day of the HW for some inland stations, although temperatures decreased during this day in coastal areas with respect to the previous 3 days. Moreover, as seen in the analysis (Fig. 5d), the model accounts for an overestimation of night-time temperatures during the study period. The model simulates the evolution of the near-surface wind speed reasonably well during the study case (Fig. S7), except for underpredicting the maximum sea breeze speed for 26 June and 5 and 9 July and overestimating the sea-breeze speed for 7 July. 3.2. Sensitivity to the initialization of soil properties Fig. 7 presents a model comparison of T2m and 10-m wind fields temporally averaged for the period of study over daytime and night-time between the simulation using HRLDAS soil data (BL-UM-HRL) and the reference using the ERA5 soil data (BL-UM). The spatial distribution of the standard deviation of the differences in T2m between the two simulations over all the time steps available at Fig. 7. Model comparison of T2m and 10 m wind vectors averaged at daytime (top row) and night-time (bottom row), for (a, d) BL-UM-HRL; (b, e) BL-UM; (c, f) difference BL-UM-HRL – BL-UM. The white colour in difference plots stands for temperature differences lower than 0.1 ◦ C in absolute value. In plots c and f, green vectors stand for the time-averaged wind vectors for BL-UM-HRL and black vectors for BL-UM. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 12 R. Segura et al. Urban Climate 39 (2021) 100982 daytime and night-time can be seen in Fig. S8. From comparison in Fig. 7, we can see that BL-UM-HRL accounts for lower daily temperature variability, with lower (higher) daytime (night-time) temperatures than those obtained from BL-UM, especially in the south and west rural regions of the AMB (Fig. 7c and f). These temperature differences are significant in rural areas where BL-UM-HRL considers higher soil moisture values than BL- UM (see Fig. 3f and 4f), such as in the west region of the AMB and the south coast. These differences between the two simulations can be seen during the entire simulation period for the T2m time series of the rural station Viladecans (see Fig. S9c). Higher soil moisture values enhance evapotranspiration from pervious areas during the daytime, altering the surface energy balance that increases latent heat fluxes in detriment of lower sensible heat fluxes. High night-time temperatures seen for BL-UM-HRL can be explained by two different reasons: 1) the HRLDAS data accounts for higher inner soil temperatures than ERA5, which implies high ground heat fluxes from the inner soil layers to the surface, and 2) moist soils have increased thermal properties (heat capacity and thermal conductivity) which implies higher heat retention in inner soil layers compared to dry soils during the day (Vahmani and Hogue, 2014). Both implies higher retention of the heat or an accumulated heat in inner soil layers which is released at night from the surface, compared to BL-UM. During the daytime, the temperature differences between the two simulations led to enhanced sea breezes in the BL-UM simulation (differences up to 0.5 ± 0.9 m s− 1) compared to those in BL-UM-HRL, although the maximum differences are below the standard deviation of wind speed. Moreover, at night, similar results were obtained by the two simulations, although the simulation without HRLDAS data accounted for slightly higher land breezes on the south-west coast of the AMB, where the temperature differences between the simulations were higher. The comparison between soil property initializations revealed different model performances for near-surface temperatures at daytime and night-time. During the daytime, the BL-UM-HRL simulation did not show a statistical improvement of the RMSE (1.60 ± 0.57 ◦ C) compared to the BL-UM simulation (1.61 ± 0.39 ◦ C) (Fig. 5a). However, the RMSE differences were significant for NAVUA (BL-UM: 1.61 ± 0.21 ◦ C, BL-UM-HRL: 1.35 ± 0.07 ◦ C) stations due to reductions in the hot bias values. Urban irrigation is a common practice in the urban parks of the AMB in summer (Ajuntament de Barcelona, 2013), which enhances the evapotranspiration while reduces the vegetation water stress. The increased soil moisture of the HRLDAS data compared to the ERA5 dataset can resemble the effect of irrigation on these areas and can explain the improvement of the BL-UM-HRL simulation in front of the BL-UM. Even though, the BL-UM-HRL showed higher RMSE in F&SNA areas (1.98 ± 0.85 ◦ C) than BL-UM at daytime, corresponding to an increase in the temperatures underpredictions. During night-time, the model showed higher RMSEs with BL-UM-HRL (1.51 ± 0.59 ◦ C) than with BL-UM (1.43 ± 0.56 ◦ C) due to an increase in the overprediction of night-time temperatures for all the stations (Fig. 5d). This increment in night-time overestimations with HRLDAS data was especially seen for AA stations, which are rural stations nearer to the south-west coast of the AMB, where higher night-time temperature differences are observed between the two simulations compared with other stations. This is related to a higher overestimation of night-time ground heat fluxes by the BL-UM-HRL simulation, due to an overestimation of inner soil temperatures by the HRLDAS data compared to the ERA5 data (see Fig. 4c), or to an overestimation of the thermal capacity and thermal conductivity of the soil (functions of the soil moisture) by the BL-UM-HRL simulation. Overall, BL-UM shows better performance in the representation of night-time temperatures, with lower RMSE and higher correlation values than those obtained for BL-UM-HRL. Fig. 8. Model comparison of T2m and 10 m wind vectors averaged at daytime (top row) and night-time (bottom row), for (a, d) BL-noUM; (b, e) BL- UM; (c, f) difference BL-noUM – BL-UM. The white colour in difference plots stands for temperature differences lower than 0.1 ◦ C in absolute value. In plots c and f, green vectors stand for the time-averaged wind vectors for BL-noUM and black vectors for BL-UM. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 13 R. Segura et al. Urban Climate 39 (2021) 100982 The statistics of WS10m for all LULC groups of stations were all comparable between the two simulations and are shown in Fig. 6. In contrast with our results, other studies have obtained improvements in near-surface variables when HRLDAS data were used for different regions (Hong et al., 2011; Rajesh et al., 2016; Sharma et al., 2017). Hong et al. (2011) obtained a relative improvement in simulations of surface energy fluxes and soil moisture, including simulations using HRLDAS data, when simulating three regions with different climatic characteristics in the United States. However, they did not see a significant effect on the simulations of surface temperatures. Additionally, Rajesh et al. (2016) simulated a heavy rainfall event in northern India, comparing a mesoscale simulation initialized with HRLDAS data and one initialized with low-resolution reanalysis data. Using HRLDAS data, they obtained better representations of soil moisture when validated with station observations, along with more accurate evolutions of the PBL charac­ teristics. Sharma et al. (2017) obtained clear improvements in simulations of T2m and WS10m when soil state variables were initialized with HRLDAS data for the simulation of a HW over the Chicago Metropolitan Area. Best and Grimmond (2014) remarked on the substantial impact of initial soil moisture on the performance of urban LSMs. To our knowledge, this is the first study of this kind for a Mediterranean climate city. Overall, the statistical analysis done in this study has revealed that the inclusion of high-resolved long spin-up soil variables in the WRF-BEP-BEM model do not imply in a sig­ nificant improvement of the model performance when a HW episode in the AMB is considered. The only significant improvement is seen for stations in urban parks at daytime, due to a better representation of urban soil conditions. Alternatively, the inclusion of HRLDAS data in the WRF model has increased the RMSE of daytime temperatures in natural regions stations and night-time tem­ peratures for all the stations. 3.3. Sensitivity to urban morphology parameters From Fig. 8 and Fig. S10, we see that BL-noUM shows higher daytime and night-time temperatures than the BL-UM simulation in the urban areas in the northern region of the AMB (up to 2.2 ± 0.7 ◦ C higher during the daytime and 2.1 ± 0.4 ◦ C at night), while lower temperatures are seen in the south during the day (up to 0.8 ± 0.5 ◦ C), corresponding to industrial areas (LCZ 10) (Fig. 8c). This difference is related to the changes in λu between the two model configurations, as mentioned in Section 2.3.2. Comparing Fig. 8c and f with Fig. S4c, a correlation between daytime temperatures and λu can be observed. This correlation can be explained by an over­ estimation of urban fractions that increase heat absorption and sensible to latent heat ratios, producing higher near-surface tem­ peratures in both daytime and night-time (Nemunaitis-Berry et al., 2017). Both simulations presented similar wind fields during the daytime (Fig. 8c), although BL-noUM presented slightly higher wind speeds in rural areas (differences up to 1.0 ± 0.4 m s− 1). This can be related to the increased temperature difference between sea and land, which causes an increased pressure difference between land and sea, enhancing the development of the sea breezes in this simulation. Considering that urban areas imply higher wind drag due to the presence of buildings (Barlow, 2014), this effect on wind speed is higher in rural areas. Moreover, at night, BL-UM presented considerably stronger wind speeds than BL-noUM (up to 1.1 ± 0.9 m s− 1). Unlike during the daytime, the increased night-time temperatures in urban areas due to λu overestimations by the BL-noUM simulation oppose the formation of a temperature gradient between land and sea, offsetting the formation of land breezes and mountain-to-valley breezes. Overall, the inclusion of specific urban morphology (BL-UM) implies a slight reduction in the RMSEs of daytime near-surface temperatures compared to BL-noUM (1.74 ± 0.43), especially for NAVUA stations (BL-UM: 1.61 ± 0.21 ◦ C; BL-noUM: 2.17 ± 0.23 ◦ C) where daytime overpredictions were considerably decreased (Fig. 5a). However, at night, the inclusion of specific morphology significantly improved the modelled near-surface temperatures in terms of RMSEs for all LULC groups (BL-UM: 1.43 ± 0.56 ◦ C, BL-noUM: 1.82 ± 0.78 ◦ C) (Fig. 5b). The overestimation of λu by the default morphological parameters used in the WUDAPT methodology produced an increase in night-time temperature overpredictions in urban areas (Fig. 5d). This overestimation especially affected the NAVUA and AA stations that were classified with a LCZ, increasing the RMSE of night-time near-surface temperatures. In the case of near-surface wind speed, the model did not show a better performance when specific urban morphology parameters were used (BL-UM), showing that the WUDAPT methodology is capable of representing the effect of buildings on wind effectively. However, BL-UM obtained considerably higher RMSEs than BL-noUM did at two stations classified as urban LCZs where wind speed was significantly underestimated. These two stations are located in open spaces without the presence of buildings in their sur­ roundings, so wind is not obstructed by the presence of buildings at these locations. Therefore, BL-noUM showed better agreement with the observations for these two stations since it considered higher wind speeds during the daytime (see Fig. 8). However, the statistical analysis shows that the two simulations are similarly correlated to the near-surface wind speed observations, for all the LULC groups during the day and night. This contrasts with the results obtained by Hammerberg et al. (2018), who found only marginal improvements using the LCZ- averaged site-specific urban morphology parameters in comparison with the WUDAPT methodology for the city of Vienna, Austria. However, significant differences arise between their study and our study. The city of Austria is characteristic of the central European cities with uniform neighbourhoods, while the AMB is a highly heterogenous region characteristic of the Mediterranean coastal cities. Moreover, Hammerberg et al. (2018), compared the average the WUDAPT method with a gridded mean value calculation of the detailed site-specific Vienna GIS data set, not the point-by-point methodology used in our study. Therefore, their simulations were not able to represent properly any variation of morphology within grid cells of the same LCZ. Additionally, the differences seen in urban fractions between the Vienna urban morphology database and the WUDAPT database was minimal and compared to the differences seen in our study, evidencing that the WUDAPT database was more fitted for the city of Vienna than for our domain. Wong et al. (2019) also pointed out that the use of tables, such as those used in Stewart and Oke (2012), instead of using local specific data for each LCZ can add a large source of uncertainty in urban simulations. Therefore, as derived from the results of our study, we conclude that it is 14 R. Segura et al. Urban Climate 39 (2021) 100982 important to use site-specific urban morphology parameters in simulations of urban climates. 3.4. Sensitivity to PBL scheme The sensitivity analysis (Fig. 9 and Fig. S11) shows that MYJ-UM outputs higher daytime rural temperatures than BL-UM (up to 1.1 ± 0.9 ◦ C) (Fig. 9c). In contrast, different results are seen over urban areas depending on the urban density. In open built areas (LCZ 6, open low-rise), MYJ-UM outputs higher near-surface temperatures than BL-UM, while in high-density built areas (LCZ 2, compact mid- rise), higher near-surface temperatures were obtained with the BouLac scheme (differences up to 1.3 ± 1.0 ◦ C). During night-time, BL- UM presented higher T2m values (differences up to 1.2 ± 0.7 ◦ C) than MYJ-UM in all the domain. Therefore, MYJ-UM showed a higher daily temperature variability than BL-UM in rural areas. From Fig. 9, we can see that the model shows considerable sensitivity to the PBL scheme in its representation of near-surface winds, especially in the sea and rural areas. Interestingly, although the temperature differences between land and sea were lower for the BL- UM simulation during the daytime, BL-UM presented greater magnitudes of sea breezes in rural areas (especially in mountain ranges, with differences up to 1.3 ± 1.2 m s− 1) and sea areas, where changes in wind direction were identified (Fig. 9c). However, MYJ-UM showed considerable sea-breeze strengths at higher levels inside the UBL (see Fig. 11 e). Over urban areas, we did not observe a significant difference in wind speed since wind speed was considerably low over these areas. In the case of night-time winds, MYJ-UM represented larger wind speeds than BL-UM in rural areas due to enhanced mountain-to-valley breezes. The modelled near-surface temperatures showed clear differences between the BL-UM and MYJ-UM simulations (Fig. 5). The simulation conducted with BouLac scheme showed lower daytime RMSE values at all stations compared to those output by the MYJ- UM (1.78 ± 0.34 ◦ C). In especial, MYJ-UM showed higher errors for UF&NVUA (1.83 ± 0.30 ◦ C) and NAVUA (1.84 ± 0.21 ◦ C) stations compared to BL-UM (1.61 ± 0.30 ◦ C and 1.61 ± 0.21 ◦ C, respectively). In the case of night-time temperatures, different performances were obtained depending on the LULC group. From Fig. 5b, we can see that BL-UM had lower night-time RMSEs for the UF&NVUA and F&SNA stations because MYJ-UM underpredicted T2m for these regions. In contrast, for the NAVUA and AA stations, higher RMSEs were seen for BL-UM than for MYJ-UM; this difference is related to the larger overprediction of night-time temperatures by BL-UM (Fig. 5d). Therefore, there is a relation between the MB observed at night-time temperatures and the PBL scheme, considering the night-time differences between the two simulations (see Fig. 9f). Even though, the MB is still high for AA (1.2 ± 0.85 ◦ C) stations, evidencing that there are other systematic errors related to different aspects of the WRF-BEP-BEM model. The choice of the PBL scheme did not significantly impact the performance of the model in the case of the wind speed. The BL-UM simulation showed lower RMSEs for UF&NVUA stations at day (BL-UM: 1.80 ± 0.54 m s− 1; MYJ-UM: 1.95 ± 0.59 m s− 1) and night (BL- UM: 1.18 ± 0.37 m s− 1; MYJ-UM: 1.40 ± 0.61 m s− 1). These reduction in the RMSE is seen in the stations placed over rooftops where BL-UM considered a lower overprediction of the wind speed compared to MYJ-UM. The wind speed at these stations is taken from the closest vertical level to the height of the meteorological station from the ground. Therefore, we can point out that the differences seen in the performance analysis of the wind speed between the two PBL schemes are related principally to the vertical distribution of the wind speed. Moreover, a difference in the performance of the two simulations is seen at daytime for the station placed over the sea, where MYJ-UM showed closer values to the observed ones than BL-UM. Overall, the BL-UM near-surface wind speeds showed higher Fig. 9. Model comparison of T2m and 10 m wind vectors averaged at daytime (top row) and night-time (bottom row), for (a, d) MYJ-UM; (b, e) BL- UM; (c, f) difference MYJ-UM – BL-UM. The white colour in difference plots stands for temperature differences lower than 0.1 ◦ C in absolute value. In plots c and f, green vectors stand for the time-averaged wind vectors for MYJ-UM and black vectors for BL-UM. (For interpretation of the ref­ erences to colour in this figure legend, the reader is referred to the web version of this article.) 15 R. Segura et al. Urban Climate 39 (2021) 100982 correlations (0.51 ± 0.12) to the observed values at night-time than the simulated by MYJ-UM (0.43 ± 0.15), although similar cor­ relations are obtained at daytime (BL-UM: 0.68 ± 0.09; MYJ-UM: 0.70 ± 0.10). Fig. 10 shows the time series of the modelled diagnosed PBL heights (PBLH) at the radiosonde launching location (nearest grid cell values to the radiosonde launching location) for the BL-UM and MYJ-UM simulations and the PBLH computed from the radio-sounding observations. Generally, both simulations reflected a rapid growth of the PBLH during the morning, reaching a maximum at midday, followed by the sinking of the inversion cape during the afternoon (Sicard et al., 2006). We observed that BL-UM considered higher diurnal PBLHs than MYJ-UM, except for on 9 July, when MYJ predicted a PBLH maximum. This difference between schemes was mostly produced by the discrepancies in the methodology used to determine the PBLH by each PBL scheme. These PBL schemes define the PBLH in different ways: MYJ defines the PBLH as the level at which the turbulent kinetic energy (TKE) decreases under a certain threshold, while BouLac defines the PBLH as the level where θv reaches the same value as that at the surface; the latter is the same method used for the observations. If we derive the PBLH in the MYJ-UM simulation using the same method as that used in the BL-UM simulation (MYJ-UM_θv), we can see that the differences between the two simulations are not as high as before. Both simulations considered similar PBLH maxima except for days 25 and 26 of June and 1 and 3 of July when BL-UM show much higher PBLH; and on 5, 6 and 7 July, when MYJ-UM_θv provided higher mixing layer tops. The differences in PBLH between BL-UM and MYJ-UM_θv can be related to the differences seen in daytime near-surface temperatures between the two simulations for these days. As it can be seen in the temporal series of T2m for the Zona Univ. station (Fig. S9b), at 1.2 km distance to the radiosonde launching location, BL-UM showed higher maximum T2m than MYJ- UM days 25 and 26 of June and 1 and 3 of July. As for days 5, 6 and 7 of July; MYJ-UM presents higher midday T2m than BL-UM. Compared to the observations, BL-UM overpredicted the PBLH for all the simulated days except for 28 June, 4, 7 and 9 July. On the other hand, MYJ-UM underpredicted the PBLH at midday for all the measurements using the TKE profile method for the diagnosis, while it overestimated all the measurements using the virtual potential temperature method (MYJ-UM_θv), except for on 28 June, 1, 2, 4 and 9 July. The overestimations of the PBLH by the BL-UM and MYJ-UM_θv can be related to the overestimations of daytime near- surface temperatures by these two simulations as it can be seen on Fig. S9b. The day with the highest underestimation of the PBLH is the 4th July, which corresponds to the only day that both simulations underpredict midday temperatures at the radiosonde launching location (Fig. S9b). Overall, the BL-UM simulations represented PBLHs close to the observations than the MYJ-UM simulation. However, when the PBLHs of the MYJ-UM simulation are computed using the virtual potential method (MYJ-UM_θv), both simulations show similar performance in the prediction of the PBL height. Banks et al. (2015) found similar results, where the BouLac represents closest PBLH at 12 UTC over the city of Barcelona for the main synoptic event of this HW (Shallow cyclone or Undetermined pressure gradient, regional recirculations). Banks et al. (2015) also used the WRF model, although without the coupling with a UCP for the urban areas. The inclusion of the BEP-BEM model and site- specific urban parameters in the PBL comparison has not implied in a modification of the results obtained by Banks et al. (2015). However, more observations recorded at different times of the day (the only available data for use in this study were recorded in the early morning or late evening) could help in the evaluation of the PBL schemes by increasing the reliability of the analysis. In this sense, different studies have used LiDAR measurements (Banks et al., 2015) to infer the mixing layer height at high temporal resolutions and almost continuously. Unfortunately, there were no available LiDAR measurements for our period of study. To analyse the differences in the representation of the vertical structure of the PBL seen in Fig. 10, the vertical profiles of θv (K) and wind speed (m s− 1) derived from the simulations along with available radiosonde observations are displayed in Fig. 11 for the three days that were considered the HW period (4, 5 and 6 July 2015) at 12 UTC. For the first day of the HW (4 July 2015), both simulations showed considerably shallow PBLHs compared with the radiosonde Fig. 10. Time series of diagnosed PBL height at the radiosounding site by BL-UM and MYJ-UM simulations and the computed PBL height from radiosounding observations (Obs.) represented with black dots. 16 R. Segura et al. Urban Climate 39 (2021) 100982 Fig. 11. Comparison of vertical profiles of WRF model simulated potential temperature (K) (top) and wind speed (m s− 1) (bottom) and radiosonde observations from the launches at 12 UTC at the days 4 July 2015 (a, d), 5 July 2015 (b, e) and 6 July 2015 (c, f). observations (Fig. 11a). We can see that both simulations underpredicted θv below 500–600 m, although MYJ-UM had a lower output (Fig. 11a). From 600 m to 2000 m, the model overpredicted θv. In general, MYJ-UM showed the closest simulated values to the ob­ servations. In the case of wind speed, all the simulations showed underestimations above the PBL (Fig. 11d), which was also seen in the other days of study (Fig. 11e and f). The underestimations seen for the wind speeds at higher altitudes were produced by un­ derestimations of the boundary conditions from the ERA5 wind fields (see Fig. S12 for the comparison of the ERA5 wind speed vertical profiles). At lower altitudes, BL-UM showed better results inside the PBL than MYJ-UM for the first day of the HW. For 5 July 2015, the PBL seemed to be more stratified than that in the previous day, possibly due to the increase in the sea breeze force on this day (see Fig. 11e). The WRF model-simulated θv values were slightly overestimated throughout the PBL (up to 1400 m), with the BL-UM simulation showing the closest values to the observations (Fig. 11b). In the case of wind speed, the MYJ scheme showed better results than the BouLac scheme compared to the radio-sounding data inside the PBL. The third day of the HW period (6 July 2015) showed the lowest observed PBLH, with a value of 850 m (Fig. 11c). This low PBLH under unstable conditions was produced by large mesoscale compensatory subsidence over the sea and the formation along the coast of a thermal internal boundary layer, as reported by Sicard et al. (2006), which is typical for the city of Barcelona in summer months and weak pressure gradients with dominating mesoscale forcing conditions (e.g. sea and land breezes). Both simulations presented vertical profiles similar to that of the observations, although with overestimations of θv below 1400 m, which cause the model to overestimate the PBLH (Fig. 11c). For the wind speed, no large difference was observed between the two simulations with underestimations inside the PBL. In general, the WRF simulations showed a cold bias tendency on the first day of the HW (4 July 2015), which changed throughout the simulation period to a warm bias towards the last day of the HW (6 July 2015). When we compare these two schemes for the entire simulation period, we see that the vertical profile of virtual potential tem­ perature using the MYJ scheme showed closer agreement with the observations for the firsts days of the simulation (25–28 June) (see Fig. S13). This is most likely due to a higher overestimation by the BL-UM simulation of the virtual potential temperature inside the PBL, especially near the surface. Between the 5th and 7th of July, the MYJ-UM overpredicted the virtual potential temperature near the surface, so the BL-UM simulation showed closer vertical profiles to the radio-sounding. For the rest of the days, except days 29 and 30 June and 4 July, both simulations show similar profiles. Therefore, we can conclude that the vertical profiles are closely related to the 17 R. Segura et al. Urban Climate 39 (2021) 100982 bias in the near-surface, and that this bias can extend to the entire PBL. We can not identify a better PBL scheme for the representation of the vertical structure of the PBL due to the high variability of the results obtained and the different performance depending on the near-surface bias. Additionally, we have found that BL-UM-HRL simulation has a better agreement than BL-UM for the vertical profile of the virtual Fig. 12. Vertical cross-section of potential temperature (K) (as shown in Fig. 1 as a red line) overlaid by wind vectors on 4 July 2015, for (a) BL-UM during day (14:00 UTC), (b) BL-crop during day (14,00 UTC), (c) BL-UM during night (22 UTC), (d) BL-crop during night (22 UTC). Overlaid in dashed black line is the PBLH diagnosed by the BouLac PBL scheme. Colls. refers to Collserola Mountain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 18 R. Segura et al. Urban Climate 39 (2021) 100982 potential temperature for almost all the studied days, due to the lower overestimation of UBL virtual potential temperatures (not shown). This agrees with the better performance of the BL-UM-HRL simulation representing midday near-surface temperatures for the Zona Univ. station (see Fig. S9b). 3.5. Sensitivity of the UBL to urbanization Fig. 12 shows a vertical cross-section of the potential temperature overlaid by wind vectors for the reference simulation (BL-UM) and the urban-to-croplands/natural vegetation mosaic-modified simulation (BL-crop) on 4 July 2015 at 14 UTC (12a and 12b) and 22 UTC (12c and 12d), the representative moments of the maximum sea breeze periods over the city (daytime) and the period of maximum UHI intensity (night-time), respectively. Overlaid in the black dashed line is the PBLH diagnosed by the BouLac scheme. During the day, the BL-UM simulation showed higher near-surface temperatures than the BL-crop simulation (up to 3.4 ± 1.1 ◦ C, see Fig. S14). Urban areas also modified the potential temperature inside the UBL, increasing the temperatures of the overlying air by convection and, consequently, increasing the mixing height (Stull, 1988). Interestingly, we can also see that the sea breeze front penetrated further inland in the BL-crop simulation than in BL-UM, even though the temperature differences between land and sea were smaller. Similar results were obtained by Sharma et al. (2017), who compared two similar scenarios but for the Metropolitan Area of Chicago and the interaction between the UHI of the city and the Lake Michigan breeze. They suggested that the reduced surface roughness length of cropland fields compared to buildings exerted a lower drag on the wind, implying higher near-surface wind speeds, as can also be seen in our case (see Fig. S14c). Moreover, Barlow (2014) reviewed different studies relating sea breezes and the UBL. They concluded that UHIs delay the passage of sea breeze fronts when passing through urban areas due to the effect of building drag on wind flow. In the BL-UM simulation, the presence of the city near the coast modified the UBL; therefore, sea breezes advected warmer temperatures inland, forming a heat plume that crossed Collserola Mountain and extended to the Vallès region (Valley). Therefore, we can observe that the spatial extent of the sea breeze thermoregulatory effect is influenced not only by the temperature difference between the land and sea but also by the surface friction of the land-use type (Martilli, 2003; Martilli et al., 2003; Barlow, 2014; Sharma et al., 2017). At night, the release of the heat accumulated during the day by urban fabric materials caused the intensification of the UHI effect, increasing temperatures inside the UBL and especially near the surface (up to 4.3 ± 1.1 ◦ C; Fig. 12c, d and S14f). This impact was higher at night than in the daytime, as would be expected due to the properties of the UHI in Barcelona city (Moreno-Garcia, 1994; Martin-Vide and Moreno-Garcia, 2020). From Fig. 12c, we can see that in the BL-UM simulation, a well-formed mixing layer above the city (PBLH ~500 m a.s.l) was sustained even after sunset. The release of sensible heat from urban areas can maintain the structure of the mixing layer due to the rise of thermals from the surface (Stull, 1988). Moreover, the presence of the UHI above the city prevented Fig. 13. Relation between urban fractions (λu) and the simulated mean 2 m temperatures (T2m) in the daytime (a) and night-time (b) during the period of study. The black dots show the mean values of the simulated T2m and λu within all the 0.05 intervals of the λu in the urban areas inside the ROI. 19 R. Segura et al. Urban Climate 39 (2021) 100982 the development of land breezes. When urban areas were modified to cropland fields, the UHI effect did not develop, and a weak land breeze could be seen (Fig. 12d and Fig. S14e). 3.6. Relation between near-surface temperatures and urban fraction To study the relation between near-surface temperatures in urban areas of the AMB and λu, an analysis similar to the one used in Li et al. (2019) was performed. We used simulated T2m values obtained from the BL-UM simulation for the analysis, as BL-UM accounted for the specific λu of each urban grid cell in its initialization of the urban morphology. To compare different urban areas with het­ erogeneous LULCs, surface properties, building morphologies and terrain heights, a zonal statistical analysis was used. All the urban classified grids in the ROI were divided into 20 groups with a constant λu interval of 0.05. The first group contained grids with λu values higher than 0 and lower or equal to 0.05; the second group contained grids with λu values higher than 0.05 and lower or equal to 0.10; and this pattern continued until the last group contained grids with λu values between 0.95 and 1.00. The average number of grids per bin is 26.6 ± 12.4; and the maximum and minimum values are 55 and 3 grids for the 0.00–0.05 and the 0.95–1.00 bins, respectively. The mean simulated T2m and λu were computed within each group. Then, a linear regression was performed to obtain the λu-linear function for T2m. The coefficient of determination (R2) of the linear regression was used to evaluate the accuracy of the linear functions. Fig. 13 shows the relation between the simulated daytime and night-time T2m and λu by computing the mean values for each 0.05-λu interval for the reference simulation (BL-UM). T2m increased with increasing λu, with a higher correlation at night (R2 = 0.96, P < 0.01) than in the daytime (R2 = 0.88, P < 0.01). The near-surface temperatures showed higher sensitivity to urban fractions during daytime, with a slope of 1.95 ◦ C than during the night-time, when the slope was 1.78 ◦ C. The intercept value should be similar to the mean of the rural T2m since it corresponds to the case where λu is 0. During the daytime, the intercept was 29.3 ◦ C, which coincided with the average simulated T2m over rural areas in the ROI for the BL-UM simulation (28.5 ± 2.2 ◦ C). The night-time interception parameter obtained from the linear regression with a value of 23.4 ◦ C also coincided with the average rural T2m obtained for the BL-UM during the whole period of study (23.9 ± 0.9 ◦ C). The differences between the surface energy budgets of urban and vegetated areas caused this observed increase in T2m with ur­ banization (Oke, 1982; Li et al., 2019). The WRF-BEP-BEM computed the T2m for each grid cell from the first-level potential tem­ perature, which is strongly influenced by the sensible heat flux. Furthermore, the sensible heat flux (SH) was computed using the linear relation between the sensible heat flux from the urban part (SHurb) and the term from the rural part (SHrural), averaged by the urban fraction as shown below: SH = λu ⋅ SH urb + (1 − λu ) ⋅ SH rural To investigate the physical mechanism that controls the relation between T2m and λu, the same methodology was used to reflect the relation between SH and λu. We obtained a linear relation with a determination coefficient of 0.96 (p < 0.01) during the daytime and 0.99 (p < 0.01) at night (Fig. S15). Therefore, even though SHurb may vary among different urban areas depending on the urban morphology parameters, the linear relation was strong since SHrural was practically constant for all the urban grid cells because the model used the same land-use class (cropland/natural vegetation mosaics) for all the pervious areas inside any urban grid cell, with the same SHDFAC and LAI values applied to all these grid cells. This influences the strong linear relation of the sensible heat flux with λu, and a linear relation was also derived for T2m; although this variable is not only dependent on the sensible heat flux, it is strongly affected by it. Additionally, the overprediction of the SHDFAC and LAI of the pervious areas in urban grid cells could explain the strong relation between λu and SH at daytime seen in this study during the daytime, which was not seen in the study of Li et al. (2019). However, the evaluation of the model performance in this study revealed that the WRF-BEP-BEM simulations showed a consid­ erable positive bias for night-time near-surface temperatures in urban grid cells with considerable vegetated areas such as urban parks and peri-urban agricultural areas (λu ranging from 0.21 to 0.65). Therefore, we can expect that the actual linear relation is not as correlated as the one that we obtained in our results, because, in reality, heavily vegetated urban areas are likely cooler than what simulated. 4. Conclusions In this study, we analysed a high-resolution simulation of a heat wave (HW) event in the Metropolitan Area of Barcelona using the WRF model coupled to the BEP-BEM urban canopy parameterization at a 1-km resolution. The overall aim was to investigate the sensitivity of high-resolution mesoscale simulations to the 1) soil moisture initialization; 2) inclusion of site-specific urban morphology parameters; and 3) planetary boundary layer (PBL) scheme. The initialization of soil properties using HRLDAS data, in comparison with low-resolution ERA5 data, significantly changed the simulations of the near-surface temperatures, decreasing the diurnal temperature range. Although the HRLDAS data showed a qualitatively better representation of the heterogeneity of soil properties with land-use in the ROI, the results of this study showed that the performance of the model only significantly improved for stations situated in urban green areas at daytime due to a better rep­ resentation of soil moisture. For all the stations, the inclusion of HRLDAS data supposed in an increase in the overestimation of night- time near-surface temperatures. Unfortunately, there are currently no soil moisture measurements available to quantitatively evaluate these soil initialization datasets. More evaluation studies are required for Mediterranean coastal cities to understand the impact of different initializations of soil properties. As we expected, including site-specific urban morphology data in the BEP-BEM model instead of default urban parameters derived 20 R. Segura et al. Urban Climate 39 (2021) 100982 based on the WUDAPT methodology significantly improved the modelled near-surface temperatures during the day and night-time, especially in urban green areas. The usage of the urban morphology parameters, as detailed in the tables of Stewart and Oke (2012) for each LCZ, overestimated urban fractions when the LCZ map, originally at a 100-m resolution, was resampled to the coarser WRF model resolution (1 km), losing important sub-grid-scale variability. Therefore, we can conclude that it is important to use site- specific urban morphology parameters of the studied region. The larger biases seen in this study correspond to stations placed in urban green areas with considerable vegetation surrounded by large urban landscapes, and the model classified the nearest grid cell to the station as an urban LCZ. Because of this, a considerable discrepancy existed between the local surroundings of these stations and the conditions that the model represented. This aspect can reflect an additional source of uncertainty and can interfere in the comparison of different sensitivity runs. Running the model at a higher resolution to fully capture urban green areas should improve the model performance, especially in green urban areas and agricultural lands close to urbanized areas. Future urban studies should be run with higher horizontal resolutions (below 1 km) to increase their representation of urban heterogeneity in the utilized models. The sensitivity study of the WRF model to the PBL scheme had a significant impact on the modelled state variables inside the PBL. The BouLac PBL scheme performed better than the MYJ scheme for near-surface temperatures for all stations during the daytime and for urban stations and forest regions at night-time. Nevertheless, MYJ showed lower overestimation of night-time temperatures in high vegetated urban stations due to a reduction of temperature overestimations. Similar performances were obtained for the modelled wind speed. Additionally, the vertical profile of the simulated virtual potential temperature was used to derive the PBLH for both simulations at midday. BouLac scheme considered much higher PBLH than MYJ, although the differences were reduced when a similar method for deriving the PBLH was used for the two simulations. Overall, BouLac scheme overestimated PBLH for midday radio­ sounding measurements of the PBL top, while MYJ underestimated all the mixing layer heights. Using the BouLac scheme provided better agreement with the PBLH observations, in line with other studies for the city of Barcelona (Banks et al., 2015; Banks and Baldasano, 2016). We also found that the vertical profiles are closely related to the bias in the near-surface temperatures, and that this bias can extend to the entire PBL. In addition to the above-mentioned sensitivity studies, we also investigated a theoretical scenario without urban areas in the AMB with the aim of studying the impact of urbanization on the UBL temperature and wind speed with a special focus on T2m. Increasing the vegetation in urban areas revealed the considerable impact of a decrease in T2m since the UHI effect does not develop under these conditions, which consequently reduces the PBLH in this scenario. Moreover, the strength of sea-land breezes in the daytime and at night was higher in the scenario without urbanization. Even though daytime UHIs favour the temperature gradient between land and sea that forms sea breezes, the surface roughness induced by urban buildings has a greater impact on wind speed, inducing a drag effect on wind flows. During night-time, the UHI of the city prevents the formation of land breezes, which can be seen in the scenario without urbanization. Moreover, in this study, we observed that urban areas also reduced the penetration distance of sea breezes in the land interior throughout the entire height of the UBL, in accordance with other studies (Barlow, 2014). Finally, we showed that the T2m field followed a linear relation with urban fractions for the urban areas inside the ROI, which is more evident at night than in the daytime. This linear relation arises from the difference in temperature between vegetated and urban areas due to the differences in their surface energy balances. However, this relation was unrealistically intensified in our study by the homogeneity of the pervious areas in the urban grid cells. Adapting the WRF-BEP-BEM model to account for sub-grid-scale LULC variability, such as by applying the MOSAIC option available in the WRF model (Li et al., 2013), could lead to improved simulations of urban areas. Overall, more studies are needed to determine the relation between near-surface urban temperatures and urban fractions accounting for specific urban morphology, sub-grid-scale variability in LULC and increased horizontal resolutions to achieve an optimal representation of landscape heterogeneity within the simulations. Understanding this relation will provide an important tool in studies of future heat-abatement urban planning strategies. Author statement Ricard Segura: Conceptualization, Methodology, Running model simulations, Interpretation of results and Writing- Original Draft and Review Editing. Sergi Ventura and Joan Gilabert: preparation of input data for model simulations, review and editing of manuscript. Alberto Martilli and Alba Badia: Supervision, Conceptualization, Review, Editing. Supervision. Gara Villalba: Conceptualization, Review Editing, Supervision, Project PI. Declaration of Competing Interest None. The research leading to these results has received funding from the European Research Council, European Horizon2020 research and innovation program under grant agreement 818002 URBAG (ERC Consolidator), awarded to Gara Villalba.. Acknowledgements This research has been done under the financial support of the ERC Consolidator Grant to the Integrated System Analysis of Urban Vegetation and Agriculture project (818002-URBAG). The authors acknowledge the financial support of the Spanish Ministry of Science, Innovation and Universities, through the “Maria de Maeztu” programme for Units of Excellence (CEX2019-000940-M). The authors thankfully acknowledge the computational resources at PICASSO and the technical support provided by the Universidad de 21 R. Segura et al. Urban Climate 39 (2021) 100982 Málaga (RES-AECT-2020-2-0004). The authors would like to appreciate the transfer of data and the services of the State Meteoro­ logical Agency (AEMET) and the Catalan Meteorological Service (SMC). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.uclim.2021.100982. References Arnfield, A.J., 2003. Two decades of urban climate research: a review of turbulence, exchanges of energy and water, and the urban heat island. Int. J. Climatol. https://doi.org/10.1002/joc.859. Ballester, J., Douville, H., Chauvin, F., 2009. Present-day climatology and projected changes of warm and cold days in the CNRM-CM3 global climate model. Clim. Dyn. 32 (1), 35–54. https://doi.org/10.1007/s00382-008-0371-0. Banks, R.F., Baldasano, J.M., 2016. Impact of WRF model PBL schemes on air quality simulations over Catalonia, Spain. Sci. Total Environ. 572, 98–113. https://doi. org/10.1016/j.scitotenv.2016.07.167. Banks, R.F., Tiana-Alsina, J., Rocadenbosch, F., Baldasano, J.M., 2015. Performance evaluation of the boundary-layer height from Lidar and the weather research and forecasting model at an urban coastal site in the north-East Iberian Peninsula. Bound.-Layer Meteorol. https://doi.org/10.1007/s10546-015-0056-2. Barlow, J.F., 2014. Progress in observing and modelling the urban boundary layer. Urban Clim. 10 (P2), 216–240. https://doi.org/10.1016/j.uclim.2014.03.011. Bechtel, B., Alexander, P.J., Böhner, J., Ching, J., Conrad, O., Feddema, J., Stewart, I., 2015. Mapping local climate zones for a worldwide database of the form and function of cities. ISPRS Int. J. Geo Inf. 4 (1), 199–219. https://doi.org/10.3390/ijgi4010199. Best, M.J., Grimmond, C.S.B., 2014. Importance of initial state and atmospheric conditions for urban land surface models’ performance. Urban Clim. 10 (P2), 387–406. https://doi.org/10.1016/j.uclim.2013.10.006. Bougeault, P., Lacarrere, P., 1989. Parameterization of orography-induced turbulence in a mesobeta-scale model. Mon. Weather Rev. 117 (8), 1872–1890. https://doi. org/10.1175/1520-0493(1989)117<1872:POOITI>2.0.CO;2. Brousse, O., Martilli, A., Foley, M., Mills, G., Bechtel, B., 2016. WUDAPT, an efficient land use producing data tool for mesoscale models? Integration of urban LCZ in WRF over Madrid. Urban Clim. 17, 116–134. https://doi.org/10.1016/j.uclim.2016.04.001. Burroughs, H.E., Hansen, S.J., 2011. Managing indoor air quality. In: Managing indoor air quality. Fairmont Press, pp. 149–151. ISBN 9780881736618. Büttner, G., Kostztra, B., Soukup, T., Sousa, A., Langanke, T., 2017. CLC2018 Technical Guidelines, 3436, pp. 0–60. Chen, F., Dudhia, J., 2001. Coupling and advanced land surface-hydrology model with the Penn State-NCAR MM5 modeling system. Part I: model implementation and sensitivity. Mon. Weather Rev. 129 (4), 569–585. https://doi.org/10.1175/1520-0493(2001)129<0569:CAALSH>2.0.CO;2. Chen, F., Manning, K.W., Lemone, M.A., Trier, S.B., Alfieri, J.G., Roberts, R., Blanken, P.D., 2007. Description and evaluation of the characteristics of the NCAR high- resolution land data assimilation system. J. Appl. Meteorol. Climatol. 46 (6), 694–713. https://doi.org/10.1175/JAM2463.1. Chen, F., Kusaka, H., Bornstein, R., Ching, J., Grimmond, C.S.B., Grossman-Clarke, S., Zhang, C., 2011. The integrated WRF/urban modelling system: development, evaluation, and applications to urban environmental problems. Int. J. Climatol. 31 (2), 273–288. https://doi.org/10.1002/joc.2158. Copernicus Climate Change Service (C3S), 2017. ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS), 13/07/2020. https://cds.climate.copernicus.eu/cdsapp#!/home. Ajuntament de Barcelona, 2013. Guidelines for Irrigation: A Practical Guide for the Irrigation of Green Areas of Barcelona. Medi Ambient i Serveis Urbans. https:// ajuntament.barcelona.cat/ecologiaurbana/sites/default/files/Manual-Reg-Parcs-Jardins.pdf. Gilabert, J., Ventura, S., Segura, R., Martilli, A., Badia, A., Llasat, C., Villalba, G., 2021a. Abating heat waves in a coastal Mediterranean city: what can cool roofs and vegetation contribute? Urban Clim. 37 https://doi.org/10.1016/j.uclim.2021.100863. Gilabert, J., Deluca, A., Lauwaet, D., Ballester, J., Corbera, J., Llasat, M.C., 2021b. Assessing heat exposure to extreme temperatures in urban areas using the Local Climate Zone classification. Nat. Hazards Earth Syst. Sci. 21, 375–391. https://doi.org/10.5194/nhess-21-375-2021. Giorgi, F., 2006. Climate change hot-spots. Geophys. Res. Lett. 33 (8) https://doi.org/10.1029/2006GL025734. Grimmond, S., 2007. Urbanization and global environmental change: Local effects of urban warming. Geogr. J. 173, 83–88. https://doi.org/10.1111/j.1475- 4959.2007.232_3.x. Hammerberg, K., Brousse, O., Martilli, A., Mahdavi, A., 2018. Implications of employing detailed urban canopy parameters for mesoscale climate modelling: a comparison between WUDAPT and GIS databases over Vienna, Austria. Int. J. Climatol. 38, e1241–e1257. https://doi.org/10.1002/joc.5447. Hong, S.-Y., Kim, J.-H., Lim, J., Dudhia, J., 2006. The WRF single moment microphysics scheme (WSM). J. Korean Meteorol. Soc. 42, 129–151. Hong, S., Lakshmi, V., Small, E.E., Chen, F., 2011. The influence of the land surface on hydrometeorology and ecology: new advances from modeling and satellite remote sensing. Hydrol. Res. 42 (2–3), 95–112. https://doi.org/10.2166/nh.2011.071. Huang, M., Gao, Z., Miao, S., Chen, F., 2019. Sensitivity of urban boundary layer simulation to urban canopy models and PBL schemes in Beijing. Meteorog. Atmos. Phys. 131 (5), 1235–1248. https://doi.org/10.1007/s00703-018-0634-1. Hunt, A., Watkiss, P., 2011. Climate change impacts and adaptation in cities: a review of the literature. Clim. Chang. 104 (1), 13–49. https://doi.org/10.1007/s10584- 010-9975-6. Hunt, J.C., Timoshkina, Y.V., Bohnenstengel, S.I., Belcher, S., 2013. Implications of climate change for expanding cities worldwide. Proceed. Instit. Civil Eng. Urban Des. Plann. 166 (4), 241–254. https://doi.org/10.1680/udap.10.00062. Iacono, M.J., Delamere, J.S., Mlawer, E.J., Shephard, M.W., Clough, S.A., Collins, W.D., 2008. Radiative forcing by long-lived greenhouse gases: calculations with the AER radiative transfer models. J. Geophys. Res.-Atmos. 113 (13) https://doi.org/10.1029/2008JD009944. Ingole, V., Marí-Dell'Olmo, M., Deluca, A., Quijal, M., Borrell, C., Rodríguez-Sanz, M., Ballester, J., 2020. Spatial variability of heat-related mortality in Barcelona from 1992-2015: a case crossover study design. Int. J. Environ. Res. Public Health 17 (7). https://doi.org/10.3390/ijerph17072553. Institute Cartographic and Geological of Catalonia (ICGC), 2019. LiDAR Data. https://www.icgc.cat/es/Descargas/Elevaciones/Datos-lidar (accessed 16 July 2019). IPCC, 2014. In: Core writing team, Pachauri, R.K., Meyer, L.A. (Eds.), Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. IPCC, Geneva, Switzerland, 151 pp. Jacob, D.J., Winner, D.A., 2009. Effect of climate change on air quality. Atmos. Environ. 43 (1), 51–63. https://doi.org/10.1016/j.atmosenv.2008.09.051. Janjic, Z., 2002. Nonsingular implementation of the Mellor-Yamada level 2.5 scheme in the NCEP Meso model. In: NCEP Office Note, 437, p. 61. Retrieved from. http://www.emc.ncep.noaa.gov/officenotes/newernotes/on437.pdf. Kain, J.S., Kain, J., 2004. The Kain - Fritsch convective parameterization: an update. J. Appl. Meteorol. 43 (1), 170–181. https://doi.org/10.1175/1520-0450(2004) 043<0170:TKCPAU>2.0.CO;2. Kottmeier, C., Biegert, C., Corsmeier, U., 2007. Effects of urban land use on surface temperature in Berlin: case study. J. Urban Plan. Dev. 133 (2), 128–137. https:// doi.org/10.1061/(asce)0733-9488(2007)133:2.(128). Li, D., Bou-Zeid, E., 2013. Synergistic interactions between urban heat islands and heat waves: the impact in cities is larger than the sum of its parts. J. Appl. Meteorol. Climatol. 52 (9), 2051–2064. https://doi.org/10.1175/JAMC-D-13-02.1. Li, D., Bou-Zeid, E., Barlage, M., Chen, F., Smith, J.A., 2013. Development and evaluation of a mosaic approach in the WRF-Noah framework. J. Geophys. Res.-Atmos. 118 (21), 11,918–11,935. https://doi.org/10.1002/2013JD020657. 22 R. Segura et al. Urban Climate 39 (2021) 100982 Li, H., Zhou, Y., Wang, X., Zhou, X., Zhang, H., Sodoudi, S., 2019. Quantifying urban heat island intensity and its physical mechanism using WRF/UCM. Sci. Total Environ. 650, 3110–3119. https://doi.org/10.1016/j.scitotenv.2018.10.025. Liang, S., 2001. Narrowband to broadband conversions of land surface albedo I algorithms. Remote Sens. Environ. 76 (2), 213–238. https://doi.org/10.1016/S0034- 4257(00)00205-4. Lindberg, F., Grimmond, C.S.B., Gabey, A., Huang, B., Kent, C.W., Sun, T., Zhang, Z., 2018. Urban multi-scale environmental predictor (UMEP): an integrated tool for city-based climate services. Environ. Model. Softw. 99, 70–87. https://doi.org/10.1016/j.envsoft.2017.09.020. Lionello, P., Malanotte-Rizzoli, P., Boscolo, R., Alpert, P., Artale, V., Li, L., Xoplaki, E., 2006. The Mediterranean Climate: An Overview of the Main Characteristics and Issues. Developments in Earth and Environmental Sciences. Elsevier B.V. https://doi.org/10.1016/S1571-9197(06)80003-0. Liu, Y., Chen, F., Warner, T., Basara, J., 2006. Verification of a mesoscale data-assimilation and forecasting system for the Oklahoma City area during the joint urban 2003 field project. J. Appl. Meteorol. Climatol. 45 (7), 912–929. https://doi.org/10.1175/JAM2383.1. Martilli, A., 2003. A two-dimensional numerical study of the impact of a city on atmospheric circulation and pollutant dispersion in a coastal environment. Bound.- Layer Meteorol. 108 (1), 91–119. https://doi.org/10.1023/A:1023044100064. Martilli, A., Clappier, A., Rotach, M.W., 2002. An urban surface exchange parameterisation for mesoscale models. Bound.-Layer Meteorol. https://doi.org/10.1023/A: 1016099921195. Martilli, A., Roulet, Y.A., Junier, M., Kirchner, F., Rotach, M.W., Clappier, A., 2003. On the impact of urban surface exchange parameterisations on air quality simulations: the Athens case. Atmos. Environ. 37 (30), 4217–4231. https://doi.org/10.1016/S1352-2310(03)00564-8. Martilli, A., Brousse, O., Ching, J., 2016. WUDAPT to WRF (W2W): urbanized WRF modeling using WUDAPT. ResearchGate, (march). https://doi.org/10.13140/ RG.2.1.3405.2724. Martin-Vide, J., Moreno-Garcia, M.C., 2020. Probability values for the intensity of Barcelona’s urban heat island (Spain). Atmos. Res. 240 https://doi.org/10.1016/j. atmosres.2020.104877. McCarthy, M.P., Best, M.J., Betts, R.A., 2010. Climate change in cities due to global warming and urban effects. Geophys. Res. Lett. 37 (9) https://doi.org/10.1029/ 2010GL042845. Meehl, G.A., Tebaldi, C., 2004. More intense, more frequent, and longer lasting heat waves in the 21st century. Science 305 (5686), 994–997. https://doi.org/ 10.1126/science.1098704. Miró, J.R., Pepin, N., Peña, J.C., Martin-Vide, J., 2020. Daily atmospheric circulation patterns for Catalonia (Northeast Iberian Peninsula) using a modified version of Jenkinson and Collison method. Atmos. Res. 231 https://doi.org/10.1016/j.atmosres.2019.104674. Moreno-Garcia, M.C., 1994. Intensity and form of the urban heat island in Barcelona. Int. J. Climatol. 14 (6), 705–710. https://doi.org/10.1002/joc.3370140609. Ndossi, M.I., Avdan, U., 2016. Application of open source coding technologies in the production of land surface temperature (LST) maps from Landsat: a PyQGIS plugin. Remote Sens. 8 (5) https://doi.org/10.3390/rs8050413. Nemunaitis-Berry, K.L., Klein, P.M., Basara, J.B., Fedorovich, E., 2017. Sensitivity of predictions of the urban surface energy balance and heat island to variations of urban canopy parameters in simulations with the WRF model. J. Appl. Meteorol. Climatol. 56 (3), 573–595. https://doi.org/10.1175/JAMC-D-16-0157.1. Oke, T.R., 1982. The energetic basis of the urban heat island. Q. J. R. Meteorol. Soc. 108 (455), 1–24. https://doi.org/10.1002/qj.49710845502. Pappaccogli, G., Giovannini, L., Cappelletti, F., Zardi, D., 2017. Sensitivity of WRF/urban simulations to urban morphology parameters: a case study in the city of Bolzano. In: In Building Simulation Applications. Free University of Bozen Bolzano, pp. 433–440. Vol. 2017-February. Pérez, C., Jiménez, P., Jorba, O., Sicard, M., Baldasano, J.M., 2006. Influence of the PBL scheme on high-resolution photochemical simulations in an urban coastal area over the Western Mediterranean. Atmos. Environ. 40 (27), 5274–5297. https://doi.org/10.1016/j.atmosenv.2006.04.039. Pineda, N., Jorba, O., Jorge, J., Baldasano, J.M., 2004. Using NOAA AVHRR and SPOT VGT data to estimate surface parameters: application to a mesoscale meteorological model. Int. J. Remote Sens. 25 (1), 129–143. https://doi.org/10.1080/0143116031000115201. Rajesh, P.V., Pattnaik, S., Rai, D., Osuri, K.K., Mohanty, U.C., Tripathy, S., 2016. Role of land state in a high resolution mesoscale model for simulating the Uttarakhand heavy rainfall event over India. J. Earth Syst. Sci. 125 (3), 475–498. https://doi.org/10.1007/s12040-016-0678-x. Ribeiro, I., Martilli, A., Falls, M., Zonato, A., Villalba, G., 2021. Highly resolved WRF-BEP/BEM simulations over Barcelona urban area with LCZ. Atmos. Res. 248. https://doi.org/10.1016/j.atmosres.2020.105220. Ryu, Y.H., Baik, J.J., 2012. Quantitative analysis of factors contributing to urban heat island intensity. J. Appl. Meteorol. Climatol. 51 (5), 842–854. https://doi.org/ 10.1175/JAMC-D-11-098.1. Salamanca, F., Martilli, A., 2010. A new building energy model coupled with an urban canopy parameterization for urban climate simulations-part II. Validation with one dimension off-line simulations. Theor. Appl. Climatol. https://doi.org/10.1007/s00704-009-0143-8. Salamanca, F., Martilli, A., Tewari, M., Chen, F., 2011. A study of the urban boundary layer using different urban parameterizations and high-resolution urban canopy parameters with WRF. J. Appl. Meteorol. Climatol. 50 (5), 1107–1128. https://doi.org/10.1175/2010JAMC2538.1. Salamanca, F., Martilli, A., Yagüe, C., 2012. A numerical study of the urban Heat Island over Madrid during the DESIREX (2008) campaign with WRF and an evaluation of simple mitigation strategies. Int. J. Climatol. https://doi.org/10.1002/joc.3398. Schatz, J., Kucharik, C.J., 2014. Seasonality of the urban heat island effect in Madison, Wisconsin. J. Appl. Meteorol. Climatol. 53 (10), 2371–2386. https://doi.org/ 10.1175/JAMC-D-14-0107.1. Schatz, J., Kucharik, C.J., 2015. Urban climate effects on extreme temperatures in Madison, Wisconsin, USA. Environ. Res. Lett. 10 (9) https://doi.org/10.1088/1748- 9326/10/9/094024. Seibert, P., Beyrich, F., Gryning, S.E., Joffre, S., Rasmussen, A., Tercier, P., 2000. Review and intercomparison of operational methods for the determination of the mixing height. Atmos. Environ. https://doi.org/10.1016/S1352-2310(99)00349-0. Elsevier Science Ltd. Servei Meteorològic de Catalunya (SMC), 2015. Butlletí climàtic mensual (juliol del 2015). Departament de Territori i Sostenibilitat. https://static-m.meteo.cat/ wordpressweb/wp-ontent/uploads/2015/06/08111706/Butllet%C3%AD-juny15.pdf. Sharma, A., Fernando, H.J.S., Hellmann, J., Chen, F., 2014. Sensitivity of WRF Model to Urban Parameterizations, with Applications to Chicago Metropolitan Urban Heat Island. American Society of Mechanical Engineers, Fluids Engineering Division (Publication) FEDSM. https://doi.org/10.1115/FEDSM2014-21292. Sharma, A., Conry, P., Fernando, H.J.S., Hamlet, A.F., Hellmann, J.J., Chen, F., 2016. Green and cool roofs to mitigate urban heat island effects in the Chicago metropolitan area: evaluation with a regional climate model. Environ. Res. Lett. 11 (6) https://doi.org/10.1088/1748-9326/11/6/064004. Sharma, A., Fernando, H.J.S., Hamlet, A.F., Hellmann, J.J., Barlage, M., Chen, F., 2017. Urban meteorological modeling using WRF: a sensitivity study. Int. J. Climatol. https://doi.org/10.1002/joc.4819. Sicard, M., Pérez, C., Rocadenbosch, F., Baldasano, J.M., García-Vizcaino, D., 2006. Mixed-layer depth determination in the Barcelona coastal area from regular lidar measurements: methods, results and limitations. Bound.-Layer Meteorol. 119 (1), 135–157. https://doi.org/10.1007/s10546-005-9005-9. Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Liu, Z., Berner, J., Wang, W., Powers, J.G., Duda, M.G., Barker, D.M., Huang, X.-Y., 2019. A description of the advanced research WRF version 4. NCAR tech. In: Note NCAR/TN-556+STR. https://doi.org/10.5065/1dfh-6p97, 145 pp. Smets, B., Verger, A., Camacho, F., Van der Goten, R., Jacobs, J., 2019. Product User Manual, Vegetation and Energy, Collection 1km, Version 2. Copernicus Publication – Product User Manual. Retrieved from. https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_PUM_FCOVER1km-V2_I1.33. pdf. Soler, M.R., Arasa, R., Merino, M., Olid, M., Ortega, S., 2011. Modelling Local Sea-breeze flow and associated dispersion patterns over a coastal area in north-East Spain: a case study. Bound.-Layer Meteorol. 140 (1), 37–56. https://doi.org/10.1007/s10546-011-9599-z. Stewart, I.D., Oke, T.R., 2012. Local climate zones for urban temperature studies. Bull. Am. Meteorol. Soc. 93 (12), 1879–1900. https://doi.org/10.1175/BAMS-D-11- 00019.1. Stewart, I.D., Oke, T.R., Krayenhoff, E.S., 2014. Evaluation of the “local climate zone” scheme using temperature observations and model simulations. Int. J. Climatol. 34 (4), 1062–1080. https://doi.org/10.1002/joc.3746. Stott, P.A., Stone, D.A., Allen, M.R., 2004. Human contribution to the European heatwave of 2003. Nature 432 (7017), 610–614. https://doi.org/10.1038/ nature03089. 23 R. Segura et al. Urban Climate 39 (2021) 100982 Stull, R.B., 1988. An introduction to boundary layer meteorology. In: An Introduction to Boundary Layer Meteorology. https://doi.org/10.1007/978-94-009-3027-8. Stull, R.B., 2000. Meteorology for Scientists and Engineers, 2nd edition. Brooks/Cole Thomson Learning, Pacific Grove. Teixeira, J.C., Fallmann, J., Carvalho, A.C., Rocha, A., 2019. Surface to boundary layer coupling in the urban area of Lisbon comparing different urban canopy models in WRF. Urban Clim. 28 https://doi.org/10.1016/j.uclim.2019.100454. UN: United Nations, Department of Economic and Social Affairs, Population Division, 2019. World Urbanization Prospects: The 2018 Revision (ST/ESA/SER.A/420). United Nations, New York. United States Geological Survey (USGS), 2019. Earth Explorer. https://earthexplorer.usgs.gov/ (accessed 16 July 2019). Vahmani, P., Hogue, T.S., 2014. Incorporating an urban irrigation module into the Noah land surface model coupled with an urban canopy model. J. Hydrometeorol. 15 (4), 1440–1456. https://doi.org/10.1175/jhm-d-13-0121.1. Vermote, E., Justice, C., Claverie, M., Franch, B., 2016. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 185, 46–56. https://doi.org/10.1016/j.rse.2016.04.008. Wong, M.M.F., Fung, J.C.H., Ching, J., Yeung, P.P.S., Tse, J.W.P., Ren, C., Cai, M., 2019. Evaluation of uWRF performance and modeling guidance based on WUDAPT and NUDAPT UCP datasets for Hong Kong. Urban Clim. 28. https://doi.org/10.1016/j.uclim.2019.100460. Zhang, P., Imhoff, M.L., Wolfe, R.E., Bounoua, L., 2010. Characterizing urban heat islands of global settlements using MODIS and nighttime lights products. Can. J. Remote Sens. 36 (3), 185–196. https://doi.org/10.5589/m10-039. 24