Methodology to separate urban from regional heat advection by use of the Weather Research and Forecasting mesoscale model

Recent studies have identified the significance of urban heat advection (UHA) as the process whereby heat, originally generated through urban modifications to the Earth's surface, is transported downwind of urban areas. Current techniques to separate UHA from local heat signals do not exclude the additional potential impacts of regional heat advection (RHA). For example, large‐scale coastal effects, in addition to latitude and longitude variations, could cause downwind temperature gradients to exist. In this study, the numerical Weather Research and Forecasting (WRF) model coupled with the Building Effect Parameterization (BEP) urban scheme is used to simulate meteorological fields for Birmingham, UK, at a high horizontal resolution (1 km2). The model is run over six case‐studies to provide over 1600 h of simulations (called ‘urban‐case’), and evaluated using a unique high‐resolution dataset from 32 weather stations across Birmingham. The UHA component is decomposed from RHA by conducting a second set of simulations (called ‘rural‐case’), where all urban land‐use is replaced with vegetation. Simulated directional ‘rural‐case’ time‐mean temperature fields, that show RHA, are then subtracted from the equivalent ‘urban‐case’ time‐mean fields. This effectively separates UHA from RHA and shows that a significant portion of heat, previously attributed to UHA in mesoscale modelling, is found to be due to RHA. Using the new methodology, a UHA signal up to 1.9 °C is found largely confined to within, and several kilometres downwind of, the urban areas. These UHA effects highlight the importance of using wind direction segmentation when determining local climate.


Introduction
Anthropogenic warming in towns and cities at night is known as the urban heat island (UHI) effect and is recognised worldwide. Urban areas make profound changes to surface properties that contribute to a nocturnal warming that can exceed 10 • C in UK conurbations compared to surrounding rural areas (Smith et al., 2011). Modifications which affect local temperatures include: (i) reduced albedo, (ii) lack of vegetation, (iii) increased roughness, (iv) larger surface area, (v) heat fluxes from buildings and vehicles, and (vi) building geometries (i.e. radiation trapping in urban canyons). The UHI intensity (difference between urban and background rural temperature) is determined by the urban configuration and local meteorology. UHIs are most pronounced under anticyclonic conditions where calm winds and clear skies emphasise the differential heating and cooling rates between urban and rural areas.
Urban conurbations are now home to over half the world's population, higher still in heavily urbanised countries such as the United Kingdom: 81.5% (ONS (Office for National Statistics), 2013). The combined effects of excess heat and rising global temperatures means there is an ever-increasing risk to health and public services. Prolonged high temperatures can have adverse effects on human health, particularly amongst vulnerable citizens. Indeed, the severe August 2003 European heatwave is thought to be responsible for up to 70 000 excess deaths (Robine et al., 2008). Heat-health effects during heatwaves are likely to be exacerbated in urban areas where temperatures are further warmed by the UHI. It has been estimated that the UHI effect contributed to around 50% of the excess mortality in the West Midlands region of the United Kingdom during the August 2003 heatwave . Nonetheless UHIs also present a winter trade-off, e.g. reduced heating requirements (Mavrogianni et al., 2011;Li et al., 2012).
To quantify Birmingham's UHI, the United Kingdom's second most populous city (1.1 million inhabitants), several approaches have been taken: (i) remote sensing (Tomlinson et al., 2012;Azevedo et al., 2016), (ii) modelling (Heaviside et al., 2015), and (iii) observations (Unwin, 1980;Johnson, 1985;Chapman et al., 2014;Bassett et al., 2016). Urban observations, often a challenge within urban environments, have recently become available in the form of a newly installed high-density network of automatic weather stations -Birmingham Urban Climate Laboratory (BUCL: Warren et al., 2016). In addition to determining the static UHI pattern, these observations have also allowed detection of urban heat advection (UHA), i.e. the heat spread from urban areas to surroundings. UHA, or the advectioninduced UHI component, is thought to occur through either: (i) horizontal heat advection within the urban canopy layer, or (ii) horizontal advection of warm air in the urban boundary layer and then mixed downwards. Whilst UHA has been conceptualised for a number of years (Lowry, 1977) and quantification attempts have been made (Brandsma et al., 2003), it is a common limitation of recent UHI studies to confine results to cities and not consider momentum effects. This is in part due to the majority of UHI studies using a small number of fixed weather stations which are unable to capture how the UHI is advected by wind. However wind direction segmentation has been used to study downwind impacts of UHI on thunderstorm generation (Dou et al., 2015).
A methodology to isolate the UHA signal was first put forward by Heaviside et al. (2015) using mesoscale modelling. UHA was calculated by removing the time-mean UHI field from the UHI field associated with a given wind direction, since the temperature at a given location is a function of background temperature (i.e. regional), local effects (i.e. topography and land-use) and that advected from upwind sources (i.e. UHA). The methodology effectively separates the temperature created locally from that advected. A key feature of this methodology is that it only needs temperature values at the same height across the area (with a reasonable spatial resolution) as the main input; thus it can also be used to analyse observational datasets. Using this approach, Heaviside et al. (2015) found an upwind/downwind temperature difference of approximately 2.5 • C in the region around Birmingham during the period of the heatwave in August 2003. This concept was recently adapted further and expanded to cover 20 months of data from the high-density BUCL network where a mean downwind warming of up to 1.2 • C was found (Bassett et al., 2016). However, neither of these studies considers the potential impacts of any additional regional heat advection (RHA). Whilst Birmingham is the most land-locked city in the United Kingdom, large-scale coastal effects, in addition to latitude and longitude changes, could cause downwind temperature gradients to exist. These effects are not excluded in the current interpretation of the UHA signal.
This study uses the Weather Research and Forecasting (WRF) non-hydrostatic mesoscale model. WRF has a large range of applications and can be coupled with several urban parametrization schemes of varying complexity. WRF has already been extensively tested in the urban environment (Loridian et al., 2013;Heaviside et al., 2015) as well as being used to demonstrate how wind can spread heat horizontally (Chemel and Sokhi, 2012;Takane et al., 2013). However, high-resolution observation networks, necessary for model evaluation, have not always been available. This is largely due to the challenges and associated cost of network maintenance within cities (Chapman et al., 2014). The BUCL observation network addresses this issue in Birmingham and allows model evaluation to be conducted across a broad range of urban land-use types within the urban canopy layer.
By demonstrating the suitability of WRF in reproducing the UHI characteristics in Birmingham, the overall aim of this study is to create a new methodology that removes any RHA effects that were not previously accounted for. This approach, along with wind direction segmentation, will subsequently isolate the mean spatial UHA pattern. By developing this tested UHA modelling approach in Birmingham, it is anticipated that it can be applied in other cities worldwide, especially where high-density urban observations are not available. This insight into UHA will be particularly useful to help city planners combat the effects of excess heat on health and infrastructure. The results for Birmingham and potential application to other cities could also be used to assess the location of long-term climate records and whether UHA bias corrections need to be made.

WRF modelling framework
WRF (version 3.8) was configured to run using four oneway nested domains (Figure 1(a)) at 3:1 grid ratios centred over Birmingham. The coarse, outer domain 1 extends across northwest Europe at 27 km grid spacing. The fine, inner domain 4 covers an area of 91 × 91 km at 1 km resolution. Model time step was set to 120 s in the outer domain with a 3:1 ratio for the inner domains and hourly output was taken. The National Centers for Environmental Prediction (NCEP) Final (FNL) Operational Model Global Tropospheric Analyses data at 1 × 1 • horizontal and six-hourly temporal resolution are used for initial and lateral boundary conditions. This NCEP data product is widely used in high-resolution WRF urban studies (Loridian et al., 2013).
The model was set up for long-wave radiation using the Rapid Radiative Transfer Model (RRTM: Mlawer et al., 1997). The Goddard (Chou and Suarez, 1994) scheme was selected for short-wave radiation as it was shown to have the best representation of air temperatures, despite a positive short-wave radiation bias (Loridian et al., 2013). This combination of RRTM and Goddard schemes is commonly used in WRF studies (Flagg and Taylor, 2011;Lee et al., 2011). The Noah land surface model (Tewari et al., 2004) that has four vertical soil levels was used for natural rural surfaces and this was coupled with the multilayer Building Energy Parameterization (BEP: Martilli, 2002) for urban surfaces. Energy fluxes are partitioned between the two models using a tiled approach based on urban fraction. The BEP scheme calculates the horizontal and vertical urban impact on wind, temperatures and turbulent kinetic energy . Additionally, through an array of three-dimensional (3D) buildings, radiation trapping and differential heating caused by shading in urban canyons are calculated. The BEP scheme has been extensively tested in urban environments (Salamanca et al., 2011;Chemel and Sokhi, 2012;Liao et al., 2014;Gutiérrez et al., 2015;Heaviside et al., 2015). However, BEP as a standalone model only considers anthropogenic heat fluxes (AHF) through constant internal building temperatures during simulations. Finally the Bougeault and Lacarrere (1989) planetary boundary-layer parametrization scheme is used because it is designed for use with the BEP urban canopy model.

Land-use and urban canopy parameters
The US Geological Survey (USGS) 24-category land-use data were used for the outer three domains. Corine Land Cover 2012 (CLC, 100 m spatial resolution) was used for the inner domain and reclassified for rural land-use categories using equivalencies defined by Pineda et al. (2004). The CLC urban land-use categories were refined into three urban land classes used by the BEP parametrization: (i) low intensity residential (CLC: Discontinuous urban fabric; Road and rail networks and associated land; Construction sites), (ii) high-density residential (CLC: Continuous urban fabric) and (iii) commercial (CLC: Industrial or commercial units). The remaining CLC urban categories (e.g. green urban areas) were assigned to the USGS category Dryland Cropland and Pasture. Satellite imagery checks showed these categories to be largely vegetated. Furthermore, these 'urban' categories do not occupy large surface areas and therefore it would not be significant when viewed at a 1 km WRF model resolution. The resulting land-use using the dominant   category for each 1 × 1 km grid square for the inner domain 4 is shown in Figure 1 A 250 m resolution Normalized Difference Vegetation Index (NDVI) image was used to calculate gridded urban fraction for the inner domain (Figure 1(c)). This use of NDVI as a proxy for urban fraction (normalised to between 0 (rural) and 1 (urban)) uses the same approach as Bassett et al. (2016). An NDVI image from July 2013 was used to represent urban fraction. The default WRF-BEP model is configured using generalised urban parameters. However, these are not specific to the study region and BEP specifically requires 3D building and street geometry distributions for each urban land-use category. To provide this level of detail, 3D Ordnance Survey (OS) Mastermap building data are used. The OS data contain building height and areas for individual buildings across most urban regions in the United Kingdom. For each of the three BEP land-use classes in Birmingham, the mean building widths and building heights are calculated using the OS data (Table 1).
However, street width calculations are less straightforward. The OS building data effectively contain two sets of information: (i) total building area and (ii) total area. The totals and means can be calculated for each urban land-use class. Subtracting the total building area from the total area leaves the total non-building area. However, if the mean street width is calculated from total non-building area (weighted by the total number of buildings in that category), this would overestimate street widths, i.e. widths would be calculated in areas containing green space. Therefore, the urban fraction is used to weight the total area in each urban landuse category to that only covered by urban surfaces, and allows for a more accurate calculation of street widths. Street widths are important because a sub-tiling approach is implemented whereby BEP calculates fluxes in urban areas, and the remaining portion of a given grid cell is calculated by the Noah landsurface model. The remaining thermal and radiative parameters used in BEP are taken from Urban Zones to characterize Energy partitioning (UZE: Loridian and Grimmond, 2012;Loridian et al., 2013). Whilst UZE parameters are created for London it is assumed that the building stock between Birmingham and London, located approximately 100 miles (160 km) apart, does not vary considerably.

WRF simulations
In total, six simulations covering a total of 67 days were run across summer 2013 (30 April-8 May; 30 May-8 June; 2 September-23 September; 25 August-6 September) and 2014 (13 May-20 May; 20 July-27 July). The case-studies were chosen because they contain periods of stable weather conditions that are favourable to UHI development. For each simulation period a second run (called 'rural-case', in contrast to 'urban-case' for the first run) is conducted whereby the urban land-use across all domains is replaced with USGS category Dryland Cropland and Pasture. A similar approach, albeit using urban and pre-urban land-use, was used by Comarazamy et al. (2013). The rural-case simulation is used as part of the calculation to separate UHA from RHA described in section 2.5. The first 12 h of each simulation are disregarded to account for model spin-up. The hourly model output for night-time hours (2100-0500 UTC) and low wind speeds (<5 m s −1 ) are chosen for the analysis. Additionally, using the WRF output variable QCLOUD (column liquid water content), simulation hours where the sky is overcast (8 oktas) are excluded. This filtering (leaving 450 h data) is conducted to limit the analysis to weather conditions that are favourable to UHI and hence UHA development.

Surface observation networks
Observations for the model evaluation are taken from two observation networks: (i) BUCL, installed in 2013 and (ii) Met Office MIDAS surface station network. The BUCL network contains 25 automatic weather stations (Vaisala WXT520, accuracy of ±0.3 • C at 20 • C: Vaisala Weather Transmitter WXT520 User's Guide, 2012) located across Birmingham. The high-density, urban nature of the BUCL network was specifically designed to study urban climate features. Within Birmingham the network has approximately one WXT station per 3 km 2 area with observations taken at 3 m above ground across a range of different land-use types. A full description of the network, including calibration and quality-control checks can be found in Chapman et al. (2014) and Warren et al. (2016). The BUCL network is supplemented by surface observations at seven Met Office stations, accessed through the British Atmospheric Data Centre (BADC). A full list of the stations used in this study can be found in Table 2 and spatial distribution shown in Figure 1 Figure 1(b).

Urban heat advection
Observational studies have demonstrated that horizontal wind flow can transport heat generated by UHIs downwind (Brandsma et al., 2003;Bassett et al., 2016). In order to separate the modelled UHA effect from UHI pattern (Figure 2(b)), a methodology was put forward by Heaviside et al. (2015). Although its first application was for analysis of modelled output, this diagnosis type of methodology, by its nature, has an advantage of being suitable also for analysing measured data (Bassett et al., 2016). The key idea of this method is to subtract the time-mean temperature field (all wind directions, T) from the time-mean temperature field for a given wind direction ( T θ , where θ : northeast (NE), southeast (SE), southwest (SW), northwest (NW)). This effectively separates the advection-induced UHI (T θ UHA ) from the local heating component created by the underlying land-use: In order to further isolate the local UHA from regional effects, we used a modified procedure, detailed as following. Similar to Heaviside et al. (2015), each simulation hour was categorised into one of four wind directions using the cross-domain-mean 10 m wind direction. Within each wind direction group, the timemean 2 m temperature field was calculated for both 'urban-cases' (T θ urban ) and 'rural-cases' (T θ rural ). The time-mean temperature field from all four wind direction groups, shown below in Eq. (2), was then taken to create an overall time-mean temperature field for urban-case (T urban ) and rural-case (T rural ) simulations. This approach was used to account for biases in the number of wind directions in the simulations (see Table 3).
To exclude any non-urban advection components from the time-mean UHI pattern ( T), the rural-case time-mean was subtracted from the urban-case time-mean, indicated in Eq. (3). This effectively removed any local or regional heat patterns that are not caused by urbanization. The same subtraction was then processed for each directional time-mean, shown in Eq. (4). The outputs of Eq. (3) were then subtracted from Eq. (4). The resulting temperature field, T θ UHA in Eq. (5), was considered the deviation from the time-mean due to horizontal wind advecting heat from urban areas (T θ UHA ). This field does not contain the RHA effects as these will have been removed by subtraction of the rural-case time-mean.  resulting UHA (T θ UHA ) therefore represents the departure of temperature from the long-term mean due to advection of the UHI. A negative value indicates cooler than the mean, and a positive value indicates warmer (note these will switch with opposing wind directions).

Model evaluation
In order to demonstrate the suitability of using the WRF-BEP model for urban temperature predictions in Birmingham, modelled 2 m hourly air temperatures are evaluated against observations. Overall model performance (presented in Table 2) is good with RMSE values obtained for the simulations similar to those found in other WRF urban studies (e.g. Liao et al., 2014;Heaviside et al., 2015). Mean RMSE across all stations is 1.68 • C and MBE −0.33 • C. Differences between observed and simulated temperatures could be explained by inherent challenges assessing point observations against model simulations at 1 km resolution. For example, the local land-use surrounding an observation creates its own local climate, whereas the model configuration uses coarser land-use. Additionally, building geometries and distributions in the model are generalised into three categories, far simpler than the multitude of urban configurations found in reality. For example, the model performs noticeably worse at station W001 which, on inspection, is located at the edge of Sutton Park (a large semi-rural park on the periphery of the city) and therefore the model does not fully represent the local land-use at this station.

Urban heat island
UHI intensity is normally calculated by subtracting the temperature at a reference rural weather station from an urban weather station (Stewart, 2011). UHI modelling has the advantage that a rural reference can be effectively provided at every grid point, in this case by use of the 'rural-case' simulations where all urban areas have been replaced with vegetation. Two timemean temperature fields are therefore created using the mean of six urban-case and six rural-case simulations respectively. The temperature difference between the urban-case and rural-case time-mean temperature fields at each grid cell is calculated during night-time hours (2100-0500 UTC). The difference between these two time-mean patterns is interpreted as the time-mean UHI field (Figure 2(a)). However, whilst the same wind categorisations for the urban and rural cases are used, a limitation is noted that for the rural-case, removing all urban components (and associated changes in surface roughness) could change wind flows in the domain.
The largest mean UHI intensities, up to 2.9 • C, are found in the centre of Birmingham. Notable UHIs are also present in smaller conurbations surrounding Birmingham, for example Coventry. Overall the spatial UHI pattern is consistent with other studies in the region (Tomlinson et al., 2013;Heaviside et al., 2015;Azevedo et al., 2016). It should be noted the pattern in Figure 2(a) is the mean nocturnal pattern and therefore higher UHI intensities could be experienced on individual nights. The time-mean UHI field also shows evidence of UHA. This can be seen clearly as the glow in temperatures that extend outwards from Birmingham (Figure 2(a)), which does not show any directional signal.
To illustrate how the spatial time-mean UHI pattern can change with wind direction, the urban-case and rural-case time-means temperature fields are then split by wind direction (θ : NE, SE, SW, NW). This categorisation is based on the cross-domain 10 m wind direction for each simulation hour, following the methodology in Heaviside et al. (2015). Isolating the time-mean UHI by wind direction (Figure 2(b)) shows two notable features: (i) the UHI core shifts downwind of the urban centre for all wind directions, and (ii) the UHI intensity differs between wind directions. Whilst simulations are filtered into night-time hours, wind speeds less  than 5 m s −1 and excluding periods where the sky is completely overcast, differences still exist within the categories (Table 3). For example the mean wind speed is lowest from the SW and highest from the NE. These correspond to the smallest and largest UHI intensities respectively. Secondly, and similarly to the heat glow shown for the time-mean UHI in Figure 2(a), evidence of UHA can be seen at the downwind edge of all plots shown in Figure 2(b). The outwards heat spread from the urban areas appears strongest for the SW and SE cases where the wind speed is lowest, i.e. there is a higher UHI component to be advected. The presence of UHA is further explored in the following section by effectively decomposing the local heating component (UHI) from the advected heat component.

Urban heat advection
The UHA results using the original methodology put forward by Heaviside et al. (2015) are presented in Figure 3(a) here. The results show that there is an approximately diagonal divide across the domain between downwind warming and upwind cooling. The upwind/downwind temperature differences are slightly smaller than found by Heaviside et al. (2015). This could be due to the 67 nights of mild UHI events in 2013 and 2014 used here instead of 10 days of strong UHI event during the 2003 severe heatwave period used in Heaviside et al. (2015). Whilst this methodology quantifies upwind/downwind temperature differences, at larger scales it fails to account for regional-scale heat advection.
In order to address this issue, the original methodology is improved by subtracting the mean hourly temperature output of 'rural-cases' from those of 'urban-cases'. The 'rural-case' results presented for domain 4 in Figure 3(b) show it is clear that the upwind/downwind temperature difference is present in each wind direction, although at a reduced magnitude to the simulations that contain urbanization. As there is no urban land-use in this 'rural-case' simulation, there are still regional advection processes at play. These RHA results show a similar pattern to the UHA results found by Heaviside et al. (2015), and in this study we propose a modification to refine the methodology in order to separate UHA from RHA (see section 2.5).
The UHA results using the new methodology to separate UHA from RHA are presented in Figure 4. In all four wind directions a clear positive UHA (T θ UHA ) signal of up to 1 • C is found downwind of the urban areas. The UHA signal is not only restricted to Birmingham, but can be found downwind from smaller urban areas within the domain. A positive UHA value effectively contains half the advection-induced UHI, whilst negative values are a construct of half the advection-induced UHI from the opposite wind directions. The total UHA component is therefore the difference between these two values and is calculated using domains with opposing wind direction, i.e. NW-SE and NE-SW. To quantify the UHA totals, the 75th percentile and the maximum UHA across the domain are adopted. For the NW-SE direction, these two respective values are 0.5 • C and 1.4 • C, and for the NE-SW direction, these are 0.6 and 1.9 • C, respectively. The UHA signal is shown to extend up to approximately 8 km for the SW case. However, the exact UHA distance is difficult to determine due to the land-use complexity.
Although taking the mean of six case-studies has improved the directional representation of UHA, the resulting UHA pattern is not symmetrical. This may be attributed to unevenly distributed wind conditions across quadrants, as shown in Table 3. Although the UHA pattern is calculated using a minimum of 37 nocturnal simulation hours in each direction (Table 3), this suggests that longer simulations are still required. However, the present results are an improvement over Heaviside et al. (2015) where RHA was not excluded and UHA was not present in all wind directions due to a limited number of cases in certain directions over the time period studied. If wind and stability conditions were identical across all categories, the methodology would produce an opposite spatial pattern in opposing wind directions. Overall, once RHA effects are removed (Figure 4), the UHA magnitude is less extensive than found using the previous methodology (Figure 3(a)). The improved UHA methodology shows UHA to  be more confined to immediate rural surroundings and within the urban boundaries. However, to explore the exact UHA contributions from urban areas of different sizes, it is suggested that idealized cases are needed to exclude complex urban land-use patterns.

Conclusion
Previous attempts to model UHA attempted to separate local and regional effects. Whilst over short distances the RHA effect is not large, when a 91 × 91 km domain is considered it could be easy to confuse this with UHA, especially if the domain is centred over an urban area. Here, an existing methodology is extended to characterise UHA, based on separating the time-mean temperature field across an urban area from temperature fields classified by each of four wind directions. The new methodology detailed here utilizes the temperature output from 'rural-case' as well as 'urban-case' simulations, which allows us to successfully separate UHA from RHA.
The WRF-BEP model configuration is run over several time periods across the Birmingham region, totalling 67 days. This enabled a time-mean UHI field with sufficient hourly data compared with previous UHA modelling studies to be created. The method applied to the modelling results to isolate the UHA signal involves subtracting a duplicate 'rural-case' time-mean simulation in which all urban land-use types are replaced by a vegetation category. UHA effects are shown to be largely confined to within and near the edges of urban areas, but can be extended to approximately 8 km downwind. Overall the total warming caused by UHA, depending upon wind direction, is up to 1.9 • C for the simulated periods. The results show UHA is not confined to Birmingham's surroundings but is also present adjacent to small urban areas within the domain.
Whilst this new methodology is an improvement over previous modelling attempts, there are still limitations. For example, as wind direction biases are accounted for in the methodology, wind speed biases within categories still have an effect on results. Although the data are filtered to speeds less than 5 m s −1 , there are still differences in speed between wind direction categories. However, there are too few hours of data to allow further classification, even though a total of 67 days simulations are available. This is particularly noticeable in the results because the UHA at a given location contains information from opposing wind directions. To overcome this, simulations would need to be conducted for longer time periods. It would be expected that increasing the number of simulation hours would reduce these directional biases. However, this is constrained by computational requirements. Additionally, replacing all urban land-use types directly with a vegetation category may alter wind flows in the domain due to changes in surface roughness. It is also difficult to determine exact UHA contributions with distance from urban areas. This is largely because of the complexity of the urban pattern in the region, and therefore UHA could originate from multiple sources. It is therefore suggested that idealized simulations are needed to simplify the urban pattern. This has the added benefit that wind speeds and directions can be controlled, therefore eliminating any biases.
Despite limitations, this UHA modelling highlights the necessity to consider wind direction when calculating temperatures in or around urban areas. As long as there are an adequate number of simulation hours this methodology can be put into practice using mesoscale modelling on any urban area worldwide. With ever-increasing urban populations there is a compelling need to spatially quantify excess heat and its impacts on health and infrastructure. Indeed, urban heat effects can exceed climate projections over the next century. The results could also be extremely useful for siting or correcting long-term temperature measurements taken near cities that do not previously account for UHA biases.