Evaluating a building-averaged urban surface scheme in an operational mesoscale model for flow and dispersion

Atmospheric Environment 88 (2014) 47e58

Contents lists available at ScienceDirect

Atmospheric Environment journal homepage: www.elsevier.com/locate/atmosenv

Evaluating a building-averaged urban surface scheme in an operational mesoscale model for flow and dispersion Ashok K. Luhar*, Marcus Thatcher, Peter J. Hurley Centre for Australian Weather & Climate Research, CSIRO Marine and Atmospheric Research, PMB 1, Aspendale, Victoria 3195, Australia

h i g h l i g h t s
A building-averaged urban canyon scheme in a mesoscale model is evaluated.
This scheme simulates the observed near-neutral to weakly unstable conditions at night.
In contrast, the original slab scheme predicts weakly stable conditions at night.
A better representation of the observed dispersion by the building-averaged scheme.
Computational efficiency of the canyon scheme is on par with the slab scheme.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 October 2013 Received in revised form 22 January 2014 Accepted 24 January 2014

A recently developed building-averaged urban surface scheme as coupled to an operational mesoscale model, TAPM, is evaluated for both flow and tracer dispersion using data from the 2002 Basel UrBan Boundary Layer Experiment (BUBBLE) conducted in the city of Basel, Switzerland. This scheme is based on the so-called town energy balance (TEB) approach and simulates turbulent fluxes using a generic canyon geometry to resolve energy balances for walls, roads and roofs. Air conditioning to close the building energy budget, in-canyon vegetation, and the effects of recirculation and venting of air within the canyon on turbulent fluxes are included. Comparison is also made with the original urban surface scheme of TAPM based on a simple slab approach with separate urban and vegetationesoil tiles and a specified anthropogenic heat flux. The results show that the new scheme leads to an overall improvement in the prediction of surface fluxes, and is able to reproduce the observed near-neutral to weakly unstable conditions at night, which is a feature of urban meteorology. In contrast, the slab scheme predicts stable conditions at night. The observed concentration fields from the tracer experiments are better simulated using the new scheme, but because there were no nighttime tracer releases, the capability of the new scheme under full diurnal conditions could not be demonstrated. For the applications considered here, the computational efficiency of the new scheme in TAPM is on par with the slab scheme. Crown Copyright Ó 2014 Published by Elsevier Ltd. All rights reserved.

Keywords: Air pollution dispersion BUBBLE data TAPM model Town energy balance Urban boundary layer Turbulent fluxes Mesoscale modelling

1. Introduction The role of a surface scheme in a mesoscale atmospheric computer model is to describe surface-atmosphere exchanges of momentum, heat and water, which influence the structure and evolution of the atmospheric boundary layer, and consequently dispersion of pollutants and tracers. Naturally, the type of surface considered is a critical parameter governing these exchanges. Urban surfaces have been shown to induce thermal and aerodynamic modifications, such as urban heat island and city-induced

* Corresponding author. E-mail address: [email protected] (A.K. Luhar).

circulations (e.g., Oke, 1982; Grimmond and Oke, 1995; Bornstein and Lin, 2000; Martilli, 2003), which have important ramifications for air pollution transport. The formulation of surface schemes for urban surfaces has traditionally been based on a simple slab approach that describes the urban surface as a concrete layer with modified roughness length and thermal properties (e.g., Oke, 1988). The operational mesoscale model TAPM (Hurley et al., 2005) for meteorological and air quality predictions uses this approach. More recently, schemes with canyon-based representations of the urban canopy have been developed that consider the effects of buildings, roads, other artificial materials used for construction and anthropogenic emissions on the surface energy budget (e.g., Masson, 2000; Martilli et al., 2002).

1352-2310/$ e see front matter Crown Copyright Ó 2014 Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.atmosenv.2014.01.059

48

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

Urban canyon schemes can be broadly categorised into singlelayer and multi-layer formulations (Masson, 2006). An example of a single-layer scheme is the so-called town energy budget (TEB) approach of Masson (2000), which assumes a generic canyon geometry and separates the energy budget into roofs, roads and walls. All canyons in the domain under consideration are assumed to have the same height and width, are located along identical roads, and are distributed in all directions with the same probability. The two facing walls are treated identically for all processes, except for the direct solar radiation. At a given time step, the canyon orientation effects with respect to the sun or the wind direction are averaged over 360 for roads and walls, which allows the computation of averaged forcing for the road and wall surfaces (instead of resolving the energy exchange for each individual canyon). The approach includes shadowing effects and parameterises the in-canyon exchange of turbulent heat fluxes. The wind, air temperature and humidity profiles within the canyon are specified. Anthropogenic heat fluxes due to domestic heating and combustion are included, and building shapes and construction materials determine parameter values (e.g., emissivity, heat capacity and albedo). When this scheme is coupled to a mesoscale atmospheric model, the surface in the model is located almost at the roof level, the scheme is forced by only the lowest model level lying almost above the roof level, and the model only sees a constant flux layer as its lower boundary. The averaging over all canyon orientations in this scheme results in a computationally efficient formulation because relatively few individual surface energy-balance solutions are required together with simplified radiation interactions. Masson et al. (2002) and Lemonsu et al. (2004) evaluated the TEB scheme in offline mode by driving it with atmospheric data, and Lemonsu and Masson (2002) used it in a mesoscale model and applied it to the Paris area in France. Multi-layer canyon schemes (e.g., Martilli et al., 2002; Otte et al., 2004; Hamdi and Schayes, 2007) use a drag force approach to account for the vertical effects of buildings. In such a scheme, the lowest level corresponds to the real level of the ground (i.e. the road surface), and the street canyon is vertically partitioned into multiple levels above the ground with a separate energy balance equation solved at each prognostic air level inside the street canyon. Although multi-layer schemes are able to determine profiles of wind, temperature and turbulent statistics within the canyon, their main disadvantage is that the atmospheric model equations for momentum, heat, and turbulent kinetic energy (TKE) are modified and need to be solved for roads, roofs and walls separately, leading to a substantially complex coupling between atmospheric-model levels and the canyon-scheme levels and hence increased computational costs. There have also been attempts to extend the single-layer TEB scheme to account for the vertical effects of buildings in a simplified manner (e.g., Hamdi and Masson, 2008; Masson and Seity, 2009). Such a scheme incorporates a drag force approach similar to that of Martilli et al. (2002), except that the prognostic air levels inside the street canyon are independent of the atmospheric model that is coupled above at a single forcing level and only one surface energy balance per wall is resolved (rather than at each level inside the canyon). The single-layer TEB approach provides an efficient and physically realistic platform to incorporate urban surfaces into operational mesoscale models in a simple, averaged manner. It has been shown to accurately reproduce the surface energy budget, canyon air temperature and surface temperatures in urban areas (e.g., Masson et al., 2002; Lemonsu et al., 2004). In this paper, we evaluate a single-layer TEB scheme as modified by Thatcher and Hurley (2012) who coupled it to TAPM. The evaluation is carried out using both flow and tracer dispersion data from the Intensive Observation

Period (IOP) of one month of the Basel UrBan Boundary Layer Experiment (BUBBLE) conducted in the city of Basel, Switzerland, in the summer of 2002 (Rotach et al., 2005). The model results are also compared with TAPM’s original slab scheme. Hamdi and Schayes (2007) applied their mesoscale model in a single-column mode to the BUBBLE IOP data, Hamdi and Masson (2008) applied their TEB to these data by running it offline forced by measurements, and Roulet et al. (2005) simulated these data using a multi-layer canyon scheme in a single-column mesoscale model driven by measurements. Rotach et al. (2005) and Batchvarova and Gryning (2006) simulated the BUBBLE tracer dataset using a Lagrangian particle model and the standard Gaussian plume approach, respectively, driven by measurements. To our knowledge, all previous mesoscale model applications to the BUBBLE dataset have not involved canyon schemes coupled to three-dimensional atmospheric models and have not simulated the tracer dispersion experiments. 2. BUBBLE field data We use data from the IOP, 10 Junee10 July 2002, of the BUBBLE experiment (see http://www.mcr.unibas.ch/Projects/BUBBLE). The BUBBLE measurements were made with the objective of studying boundary-layer and surface-exchange processes over different types of surfaces (i.e. urban, sub-urban and rural) and their role in the transport and diffusion of air pollution. Basel is a mid-size town with a built-up area of about 130 km2 (Rotach et al., 2005). The main urban measurement tower, Basel-Sperrstrasse (Ue1), was 32m high and located inside a street canyon in an area with dense, fairly homogeneous, residential building blocks, and a mean building height of 14.6 m above ground (or road) level (AGL). Location 3 in Fig. 1 corresponds to the meteorological tower Ue1. In the vicinity of the tower, the building height was 14 m AGL and the street canyon aspect ratio (i.e. height-to-width ratio) was about unity. The surface roughness length was 2.1 m, and the zero-plane displacement height was 9.5 m (Christen and Rotach, 2004). Sonic anemometers were installed at six levels, namely 3.6, 11.3, 14.7, 17.9, 22.4 and 31.7 m AGL. To study dispersion, sulphur hexafluoride (SF6) tracer was released at near roof-level at two locations, namely R1 (18.6 m AGL) and R2 (21 m AGL), over four separate days (Rotach et al., 2004; Gryning et al., 2005) (see Table 1 for release conditions and Fig. 1 for locations). The two sources were approximately 900 m apart. There were 19 SF6 sampling locations, of which 13 were typically positioned 1.5 m above the roof level and 6 were street-level samplers, the latter in relatively open areas. The distance of the sampler closest to a source was 200 m and that farthest from a source was 2.5 km. Most samplers were located within 1.5 km from a given source. Tethered balloon soundings were also carried out for 24 h starting from the afternoon of 4 July at Basel-Messe (Site Ue3 Fig. 1) in the city, measuring profiles of wind speed, wind direction, temperature and humidity. Hourly-averaged data from the fixed monitoring stations and the tethered balloon data are used for model comparison. The BUBBLE data reveal the distinct influence of the urban surface on flow properties (Rotach et al., 2005). Luhar et al. (2006) previously used the BUBBLE flow data in conjunction with TAPM to evaluate relationships between urban and rural near-surface meteorology. All times given here are in Local Standard Time (LST). 3. Mesoscale model The Air Pollution Model (TAPM, v4.0) developed by CSIRO (Australia) is an operational, inline, coupled prognostic

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

49

Fig. 1. Map of Basel showing the roof-level and street-level tracer samplers (circles numbered 1e19). Circle 3 is also the tower location Ue1. Ue3 (diamond) is the location of the tethered balloon soundings. The filled squares are the tracer release sites. The coordinate system is CH1903 used in Switzerland. (Adapted from http://www.mcr.unibas.ch/Projects/ BUBBLE/images/maps/tracermap.gif.)

meteorological and pollutant dispersion model (Hurley et al., 2005; Hurley and Luhar, 2009; http://www.cmar.csiro.au/research/tapm), which has previously been applied to a variety of regional- and local-scale dispersion problems (e.g., Luhar and Hurley, 2003; Luhar et al., 2008; Zawar-Reza and Sturman, 2008; Luhar and Hurley, 2012). It uses global input databases of terrain height, land use, soil type, leaf area index, monthly sea-surface temperature, and synoptic meteorological analyses, and can be used in one-way nestable mode. The meteorological component of TAPM predicts the local-scale flow against a background of larger-scale meteorology provided by input synoptic analyses. It solves momentum equations for horizontal wind components; the incompressible continuity equation for the vertical velocity in a terrain-following coordinate system; and scalar equations for potential virtual temperature, specific humidity of water vapour, cloud water/ice, rain water and snow. Explicit cloud microphysical processes are included. Pressure is determined from the sum of hydrostatic and optional nonhydrostatic components, and a Poisson equation is solved for the non-hydrostatic component. We use the default hydrostatic option. Measurements of wind speed and direction can optionally be assimilated into the momentum equations as nudging terms. The turbulence closure terms in the mean equations use a gradient diffusion approach, including a counter-gradient term for the heat flux, with eddy diffusivity determined using prognostic equations for the turbulent kinetic energy (E) and its dissipation rate (3 ). A vegetative canopy, soil scheme, and urban scheme (see below) are

used at the surface, while radiative fluxes, both at the surface and at upper levels, are also included. Surface boundary conditions for the turbulent fluxes are determined via the MonineObukhov similarity theory and parameterisations for stomatal resistance. The air pollution component uses the predicted meteorology and turbulence from the meteorological component, and consists of an Eulerian grid-based set of prognostic equations for pollutant concentration and an optional Lagrangian mode (used here); the latter to allow a more detailed accounting of near-field dispersion within the innermost nest for pollution. Instead of using a full particle model, the Lagrangian mode uses Hurley’s (1994) particlepuff (PARTPUFF) approach, whereby the released mass at a given instant is represented as a circular puff (or disc) in the horizontal plane and the vertical motion of this puff is determined via the particle approach. The use of a puff approach in the horizontal makes the computations faster without significant loss of accuracy, whereas retaining particle framework in the vertical allows a proper treatment of important flow characteristics such as the vertical turbulence inhomogeneity and wind shear, which a traditional puff model cannot address satisfactorily. The puff position in the horizontal plane is updated through advection by the modelled ambient wind, and the puff spreads in the streamwise and crosswind directions are calculated using Taylor’s statistical diffusion theory given the modelled turbulence levels. The puff position in the vertical is determined from a nonstationary Langevin equation using the modelled mean vertical velocity and turbulence (Luhar and Hurley, 2003). A number of

50

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

Table 1 Details of the SF6 release conditions and background concentrations. Date

Source

Release period (LST)

Release rate (g s1)

Background concentration (ng m3)

26 June 4 July 7 July 8 July

R1 R2 R1 R1

12:00e16:00 14:40e18:00 13:10e17:00 14:00e18:00

0.0503 0.0499 0.3008 0.1319

33.6 32.2 34.0 53.2

puffs are released sequentially to simulate a continuous release over a specified period of time. Once puffs have travelled a certain length of time (default is 15 min), they are no longer tracked and their mass is converted to grid concentration and put onto the Eulerian grid. From then on the plume is spatially large enough to be satisfactorily described by the Eulerian transport equation. 3.1. Current slab-based urban scheme The current slab scheme of TAPM treats the urban cover as a concrete plate with modified roughness length and thermal properties that are appropriate for such a surface (Oke, 1988; Pielke, 2002; Hurley, 2008). It calculates surface moisture, surface temperature and surface fluxes of momentum, heat and moisture for bare soil, vegetation cover, and urban cover separately, and then uses a weighting scheme according to the fraction of the area covered by the three surfaces in order to derive the effective surface values of these parameters. The evaporative heat flux and specific humidity for the urban fraction are assumed to be zero. A zeroplane displacement height for each surface is included. An anthropogenic heat flux is included in the soil, vegetation and urban surface flux equations. Urban surface-layer scaling variables are calculated using the same approach as for soil and vegetation, incorporating the corresponding urban roughness length. TAPM uses US Geological Survey’s land-use data given at approximately 1 km resolution. The generic urban land-use category contained in the default databases can be thought of as medium density urban conditions. For the BUBBLE experiment, we have used the default parameters as specified by: fraction of the urban cover ¼ 0.5, albedo ¼ 0.15, urban anthropogenic heat flux ¼ 30 W m2, and the urban roughness length ¼ 1.0 m (Oke, 1988; Pielke, 2002). 3.2. New canyon-based urban scheme Thatcher and Hurley (2012) modified the TEB approach of Masson (2000) and coupled it to TAPM. The modifications include an efficient and integrated big-leaf model to represent in-canyon vegetation and a parameterisation for air-conditioners for energy conservation when the air temperature is higher than the building comfort temperature. The big-leaf model is based on the work of Kowalczyk et al. (1994), but replacing the soil component with road surfaces. This approach only requires two additional prognostic variables to constrain the amount of water available for evaporation and transpiration (i.e., the leaf water reservoir and the total soil moisture), whereas the canopy temperature is solved iteratively to close the vegetation energy budget. Canyon turbulent heat fluxes between walls, roads and in-canyon vegetation are parameterised according to a modified version of Harman et al.’s (2004) aerodynamic resistance network, which considers recirculation and venting of air within the canyon. This is shown in Fig. 2, in which the bold lines indicate venting and recirculation regions of the air flow in the canyon. The present approach requires separation of the energy budgets of the two walls into an easterly facing component

Fig. 2. Schematic representation of the aerodynamic resistance (U) network used within the canyon in the new urban scheme.

and a westerly facing component (whereas the original TEB model employs a single wall energy budget after averaging the canyon fluxes over 360 of possible canyon orientations). The two wall energy budgets are derived by averaging the canyon fluxes over 180 of possible canyon orientations. Since walls with an easterly facing component have a different temperature diurnal cycle compared to westerly walls, there is an aggregate representation of a temperature differential across the canyon when considering the recirculation of air in the canyon. The parameterisation is then a more process based representation of turbulent heat transfer within the canyon compared to the more empirical based representation used in the original TEB model. In our simulations, the urban surface characteristics were chosen to be as much consistent as with those reported by Christen and Vogt (2004) for the BUBBLE Ue1 site. Following are the parameter values required and used in the new scheme for the BUBBLE simulation: in-canyon vegetation fraction ¼ 0.16, area fraction occupied by buildings ¼ 0.54, mean building height (zH) ¼ 14.6 m, building height to canyon width ratio ¼ 1.0, ratio of roughness length to building height ¼ 0.05, roughness length of in-canyon surfaces ¼ 0.1 m, industrial sensible heat flux ¼ 0 W m2 (default), and daily averaged traffic sensible heat flux ¼ 1.5 W m2 (default). 3.3. Mesoscale model configuration The input synoptic fields of the horizontal wind components, temperature and moisture required as input boundary conditions for the outermost nest were obtained from the Australian Bureau of Meteorology’s GASP (Global AnalySis and Prediction) system with a resolution of 1 longitude  1 latitude at 6-hourly intervals. TAPM was run with four nested domains of 35  35 horizontal grid points at 20-km, 7.5-km, 2-km and 0.5-km spacing for meteorology, and 41  41 horizontal grid points at 2-km, 0.75-km, 0.2-km and 0.05km spacing for pollution/dispersion, all centred at the location (7 360 E, 47 340 S) which is equivalent to 612.144 km east and 268.452 km north in the CH1903 geodetic datum used in Switzerland and is very close to the location of the Ue1 urban tower (see Fig. 1). The number of vertical levels in the model is the same for meteorology and dispersion. The lowest ten of the 25 vertical levels were 10, 25, 50, 100, 150, 200, 250, 300, 400 and 500 m, with the highest model level at 8000 m. The terrain height within the innermost model domain varies from 240 to 660 m above mean sea level. An initial value of 0.25 m3 m3 for the deep soil volumetric moisture content (at a depth of about 1 m from the surface) was used (which is only relevant for the slab scheme). The land-use category is extracted from the land-use database for each grid cell. For an urban grid cell, TAPM’s urban scheme calculates average fluxes, which are only coupled to TAPM’s lowest

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

level lying above the roof level. Thus, TAPM only sees a ‘flat’ constant flux layer as its lower boundary. There are various time steps in the model (see Hurley et al., 2005, particularly their Fig. 1). A large time step of 300 s is used to update components such as surface temperature and fluxes (via the surface scheme) and radiation. The model advection time step depends on the horizontal grid spacing, and therefore can be different for meteorology and pollution and for different nests. For the innermost-nest spacings used here, the time step was 12.5 s for meteorology and 2.5 s for pollution (these values become progressively larger for the outer nests). The modelled meteorological variables are linearly interpolated to the pollution grid and passed to the pollution component every large time step. The hourlyaveraged predictions on the innermost meteorological nest (500 m resolution) and pollution nest (50 m resolution) were used for comparison with the data.

4. Results and discussion Below, we present the temperature, fluxes, wind fields and tracer dispersion results obtained using the two urban schemes in TAPM. Commonly used model performance statistics used are: the linear correlation coefficient (r), the fraction of predictions within a factor of two of the observations (FAC2), and the normalised mean square error (NMSE), the fractional bias (FB) and the index of agreement (d) (Willmott, 1981) defined as follows:

NMSE ¼  FB ¼ 2 

ðO  PÞ2 O P OP OþP

d ¼ 1

;

(1)

 ;

ðP  OÞ2 ; 2 P  O þ O  O

(2)

(3)

where O and P refer to the observed and predicted concentrations, respectively, and an overbar denotes the mean value. FB varies between 2 and 2, with a positive value implying underprediction and good model performance indicated by a value close to zero. Unlike r, which is a measure of the linear correlation (or dependence) between two variables, d is sensitive to differences between the observed and model means as well as to certain changes in proportionality and varies between 0 (no agreement) and 1 (perfect agreement).

4.1. Meteorology and fluxes The new urban scheme is configured so that the first TAPM atmospheric model level is 10 m above the zero-plane displacement height of the urban area (estimated as 8.8 m AGL), or approximately 5.8 m above the BUBBLE building height and 18.8 m AGL. The first level in TAPM with the default slab scheme corresponds to about 16.7 m AGL. The 17.9-m observational level is the first level clearly above the BUBBLE roof-top level and reasonably matches the first model level for comparison. Fig. 3 compares the predicted time series of temperature at the first model level with the observations at 17.9 m. The overall variation is well simulated by both urban schemes in TAPM. A closer inspection reveals that, when not considering the observed outliers that greatly deviate from the

51

Fig. 3. Time series of the observed hourly-averaged temperature at 17.9 m AGL at the Ue1 urban site. The corresponding time series predicted by TAPM using the current and new urban schemes are shown by dotted and solid lines, respectively.

time series, many of the daily minimum temperature values are slightly better predicted by the new scheme. The urban roughness sublayer roughly extends from the average building height to 2e5 times that height (Rotach et al., 2005), with an overlying inertial sublayer. Both urban schemes assume a constant flux one-layer urban treatment, and do not resolve the roughness sublayer. Thus, the predicted fluxes characterise the inertial sublayer, akin to the surface layer over flat surfaces, where turbulent fluxes reach their maximum values. We compare the modelled sensible heat flux (Hs) with the measurements at the level where the magnitude of Hs is the largest. Hence, the implicit assumption here is that the level at which the flux magnitudes are the largest is representative of the modelled inertial sublayer flux. On average, the observed Hs at 22.4 m is about 15% higher than that at 17.9 m during 0000e0900 h and is lower by the same amount at the other times. The average observed Hs at 31.7 m is about 30% smaller than that at the above two levels, and decreases rapidly within the canyon. We use the Hs data at 22.4 m for comparison with the model simulations. Fig. 4 shows the observed and modelled average diurnal variations of Hs over the experimental period. The vertical bars around the points represent one standard deviation of scatter. The observed Hs is positive (i.e. upwards) at all times. The variation predicted by the current scheme (Fig. 4b) is in reasonably good agreement with the data, but it is apparent that it yields zero or negative Hs during 2000e0500 h whereas the observations suggest values of about 40 W m2. The new urban scheme represents this observed nighttime behaviour of Hs much better. Near-neutral to weakly unstable conditions at night are a feature of urban meteorology, which the new scheme is able to reproduce. However, the new scheme overestimates the measured daytime peak by 45%, as compared to 20% by the current scheme. Fig. 5 shows that the new scheme describes the observed probability (or normalised frequency) distribution of Hs very well, whereas the current scheme does comparatively well only for sensible heat fluxes higher than 200 W m2. (The probability distribution is calculated by counting the hourly-averaged values in bins for the experimental period.) The new urban scheme overestimates Hs greater than 400 W m2, which is less apparent in the current urban scheme. This problem may be related to an overestimation of friction velocity as discussed below. The modelled scatter as represented by the vertical bars is slightly larger around the mid day and smaller at night than the observations for both schemes. The index of agreement (d) value of 0.86 for the current scheme is slightly higher than the value of 0.82 for the new scheme due to the somewhat better prediction of peak Hs levels during the day by

52

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

Fig. 4. Average diurnal variation of sensible heat flux (Hs) over the experimental period 10 Junee10 July 2002: (a) data, (b) modelled with the current urban scheme, and (c) modelled with the new urban scheme. The vertical bars correspond to one standard deviation.

the former. But the FAC2, NMSE and FB values of 38%, 0.64 and 0.28, respectively, for the current scheme, compared to 54%, 0.51 and 0.17, respectively, for the new scheme suggest a better performance by the latter. There are no latent heat flux (QE) data at the 17.9-m level, but two sets of such data are available at the 31.7 m-leveldone taken using a Li-Cor analyser (with 30.1% data missing) and the other using a fast hygrometer (with 40.5% data missing). Based on the height variation of Hs discussed above and of friction velocity discussed below, this height level is probably a little too high for comparison with the modelled values and the amount of missing data is considerable, so the comparison given below is rather tentative. Because of their larger sample size, we only present the Li-Cor data. Fig. 6 shows the observed average diurnal variation of QE and the modelled behaviour. Overall, there is a large overprediction of QE by the current scheme and a relatively small underestimation by the new scheme during daytime. As described in Section 3.2, the new urban scheme estimates the turbulent fluxes using a resistance network that considers recirculation and venting within the canyon, whereas the current urban scheme weights the fluxes from an urban slab tile and a vegetation canopy tile. Both urban schemes employed similar values for unconstrained canopy stomatal resistance (rs ¼ 50) and leaf area index (LAI ¼ 2). The results then suggest that integrating urban vegetation into the

Fig. 5. Probability density function (PDF) of sensible heat flux (Hs) for the data and the two urban schemes.

canyon energy and flux budgets can produce a more realistic latent heat flux than simply treating the vegetation as a separate tile. The new scheme simulates the observed scatter in QE better than the current scheme; the latter considerably overestimates the scatter around mid day and underestimates it at the other times. The modelled average diurnal variation of net radiation (Rn), which is the balance of incoming and outgoing solar and terrestrial radiation, is shown in Fig. 7 together with the observed variation based on measurements made at the 31.7-m level. The current scheme overestimates the observed peak by 28%, whereas the new scheme does that by about 16%. Both schemes perform similarly at night. A time series examination of Rn reveals that one reason the model overestimates the peak is that on some days (e.g. 16, 24 June and 2 July) the model is not predicting the cloud episodes over the region properly. It is also seen that on clear days, the current scheme consistently overestimates the peaks in Rn by about 12% whereas the peaks predicted by the new scheme closely match the observed values. This is evident in a typical time series extract shown in Fig. 8 for 16e19 June (this period was also simulated by Hamdi and Masson (2008)) in which the decrease in the measured values during the morning of 16 June is due to clouds. During 2000e0400 h, the Rn values predicted by the new scheme are generally smaller than the data and the current scheme performs better. After a series of parameter sensitivity tests, we found that the difference in the daily minimum Rn could be attributed to the difference in heat storage by the two urban schemes. Specifically, reducing the thickness of the first roof layer from the default 10 cm to 5 cm increased the value of Rn at night by about 25 W m2. However, this change also reduced the daily maximum Rn by a similar amount, thereby somewhat degrading the performance of the model during the day. This result reflects some of the challenges with the appropriate assignment of thermal parameters for urban models as reflected in the urban intercomparison study described by Grimmond et al. (2011). The values of d, FAC2, NMSE and FB are 0.93, 69%, 0.59 and 0.40, respectively, for the current scheme, and are very similar at 0.93, 66%, 0.71 and 0.20, respectively, for the new scheme. Friction velocity (u* ), which is a measure of the vertical mo2 2 mentum flux, is determined as u* ¼ ½u’w’ þ v’w’ 1=4 , where u’w’ and v’w’ are Reynolds stresses. The observed u* is maximum at the 22.4-m level, and, on average, is lower by 10% at 17.9 m and by 15% at 31.7 m. It decreases further down within the building canyon. Rotach et al. (2004) term the maximum u* from BUBBLE a characteristic velocity, noting that the traditional definition of friction velocity pertains to the inertial sublayer. One may consider that the 22.4-m level (which is about 1.5 times the building height) where the observed fluxes are near their maximum marks the transition

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

53

Fig. 6. Average diurnal variation of latent heat flux (QE) over the experimental period: (a) data, (b) modelled with the current urban scheme, and (c) modelled with the new urban scheme. The vertical bars correspond to one standard deviation.

Fig. 7. Average diurnal variation of net radiation (Rn) over the experimental period: (a) data at 31.7 m AGL, (b) modelled with the current urban scheme, and (c) modelled with the new urban scheme. The vertical bars correspond to one standard deviation.

from the roughness sublayer to the inertial sublayer, but this is somewhat inconsistent with the general understanding of the transition to the inertial sublayer being at 2e5 times the building height (Rotach et al., 2005; Batchvarova and Gryning, 2006). There are no flux observations above 31.7 m, and there are no additional sites with a similar multilevel setup as Ue1 to check for any local flux variations, which make is difficult to ascertain the transition height. We use the u* data taken at the 22.4-m level (as was the case for Hs) to compare with the model results. In Fig. 9, both schemes are able to reproduce the overall observed diurnal behaviour of u* qualitatively, with a peak in the afternoon. Of the two schemes, the new scheme estimates u* better

at night, whereas the current scheme does better during daytime. The probability density function (PDF) plot in Fig. 10 suggests that the new scheme performs noticeably better for u* < 0.2 m s1 and the current scheme for u* > 0.6 m s1. The values of d, FAC2, NMSE and FB are 0.69, 68%, 0.35 and 0.16, respectively, for the current scheme, and are 0.70, 76%, 0.26 and 0.06, respectively, for the new scheme, which suggest slightly better simulations by the new scheme. Given the Hs and u* data at 22.4 m, the variation of the stability parameter L, the Obukhov length, can be examined, where L ¼ u3* r cp =½kðg=qÞHs , q is the mean potential temperature, g is the acceleration due to gravity, k is the von Karman constant

Fig. 8. Time series of the observed hourly-average net radiation (W m2) at 31.7 m AGL for the period 16e19 June 2002. The corresponding time series predicted by TAPM using the current and new urban schemes are also shown.

54

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

Fig. 9. Average diurnal variation of friction velocity (u* ) over the experimental period: (a) data (at 22.4 m), (b) modelled with the current urban scheme, and (c) modelled with the new urban scheme. The vertical bars correspond to one standard deviation.

(¼ 0.4), cp is the specific heat of air at constant pressure, and r is the air density. The modelled L is a property of the inertial sublayer. On the other hand, the observed L is determined at the level of maximum sensible heat flux and shear stress within the influence of the roughness sublayer, and should, therefore, be viewed as a local Obukhov length. Fig. 11 presents the average diurnal variation of z/L, where z ¼ 1 m is taken as a scaling height. 1/L ¼ 0 corresponds to neutral conditions, whereas a negative or positive 1/L suggests unstable or stable conditions, respectively. In Fig. 11a, the observed 1/L is positive, with a maximum magnitude during daytime. In Fig. 11b, TAPM with the current scheme is able to simulate the daytime variation qualitatively, but during 1900e0400 h at night the simulated conditions are stable, in contrast to the observed unstable conditions. The new urban scheme represents the observations better (Fig. 11c) with unstable to near-neutral conditions at night. It is instructive to examine how the observed upper air winds compare with the modelled fields. Fig. 12 shows the wind profiles obtained from the ascents/descents of the tethered balloon measurement system during 4e5 July at Site Ue3 (Fig. 1, see Rotach et al., 2005 for details). The data are averaged into vertical layers with 20-m resolution. The flow is initially westerly, but at around 2000 h it starts to weaken and then at 2100 h the flow turns easte southeast in lower layers, probably as a result of cold air drainage from the Swiss Midlands and the High Rhine Valley. The thickness of the eastesoutheast flow grows, and a jet-like structure is evident

during 2300e800 h. For the hours 0900e1100 h, the flow is almost uniform throughout the vertical domain and remains eastesoutheast. Subsequently, the drainage winds stop, and the flow weakens and becomes rather unsteady, perhaps as a result of boundary-layer convective motions generated by a strong warming of the surface. The corresponding TAPM-predicted hourly wind fields at Ue3 are presented in Fig. 13 for both schemes. There is a good qualitative agreement between the modelled and tethered-balloon fields. Some differences may arise due to the fact that the model profiles are hourly-averaged, calculated for a single surface location, and the values at all levels correspond to the same time, whereas the sonde profiles are near-instantaneous values and may not be exactly above the same surface location. In both Fig. 13a and b, following an initial west-northwest flow, the modelled winds quickly shift at around 2100 h to east-southeast, not only near the ground but at all heights. As a consequence, the observed low wind region caused by two opposing flows at a height of around 700 m during 0100e0500 h shown in Fig. 12 is not present in the modelled fields (where it occurs earlier at 2100 h). At 1100 h, the flow turns south-east and strengthens subsequently. The unsteadiness seen in the observed winds around this time is not present in the modelled fields (this may partly be due to the fact that the data are nearinstantaneous values with inadequate averaging over convective motions). It is interesting that in Fig. 13 both urban schemes give very similar results, but it is clear that there are significant differences in the lowest levels for 0800e1200 h. During 1000e1200 h, the flow below 200 m calculated using the current scheme in TAPM has a pronounced southerly component whereas the new scheme introduces a pronounced northerly component; however, the winds are very light in both cases, and the predicted mean wind direction in these conditions can be very sensitive to how a particular scheme is formulated. 4.2. Tracer concentrations

Fig. 10. Probability density function (PDF) of friction velocity (u* ) for the data and the two urban schemes.

The observed concentrations from a series of tracer experiments minus the background values were used for comparison with the modelled concentrations. Fig. 14a presents a scatter plot of the observed vs. predicted concentrations when the current urban scheme is used, and Fig. 14b corresponds to the new urban scheme. There is considerable scatter in both plots, but the performance statistics (barring correlation coefficient r) given in Table 2 indicate that the new urban scheme results in a moderately better performance for dispersion. The index d allows for sensitivity towards differences in observed and predicted values as well as proportionality changes, and is therefore a better performance statistic than r. Note that the concentration measurements were done in the

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

55

Fig. 11. Diurnal variation of z/L (where z ¼ 1 m) averaged over the experimental period: (a) data, (b) modelled with the current urban scheme, and (c) modelled with the new urban scheme.

daytime when the relative differences between the predicted sensible heat fluxes (which determine ambient stability and turbulence) from the two schemes are not as great as at night. Therefore, we anticipate that concentration differences between the two schemes would be greater at night (discussed later). Prediction of concentration at a location due to a point source can be very sensitive to wind direction. In order to reduce the influence of any predicted inaccuracy in wind direction, we assimilated the hourly-averaged wind speed and wind direction observed at Ue1 at 22.4 m AGL in TAPM assuming a 5-km radius of influence. The assimilation method is based on the approach of Stauffer and Seaman (1994) where a nudging term is added to the horizontal momentum equations with the nudging coefficient equal to the inverse of three times the advection time step. The observed wind speeds were extrapolated to the lowest six model levels using a similarity relationship (see Luhar et al., 2006). The extrapolated values were then assimilated into the model assuming that the wind direction was the same at all six model levels. An improvement in the modelled concentrations is obtained in Fig. 15a compared to Fig. 14a when the wind data are assimilated in TAPM with the current urban scheme, which is also apparent from the performance statistics given in Table 2. Fig. 15b corresponds to the new urban scheme with wind data assimilation. The current scheme tends to systematically overestimate observed concentrations above 50 ng m3. Overall, Fig. 15b is the best model configuration of the four (see Table 2). The above model performance can be compared with other modelling studies that have used the same BUBBLE concentration

data. Rotach et al. (2004) simulated these data using a full 3-D Lagrangian particle model that takes into account characteristics of the urban roughness sublayer. They used measured and parameterised flow and turbulence statistics in their model, and compared hourly-averaged concentrations for 26 June and threehour averaged concentrations for all four tracer experiments. Rotach et al. (2005) used the same model, driven by the measured flow and turbulence fields, and compared hourly-averaged concentrations for all four tracer experiments. Batchvarova and Gryning (2006) applied a simple Gaussian plume model and compared the maximum arcwise concentrations for two arcs for 26 June. Since Rotach et al. (2005) present hourly-averaged

Fig. 12. Temporal variation of the wind profile obtained using the tethered balloon soundings at City Centre (Site Ue3) during 4e5 July.

Fig. 13. Modelled temporal variation of the wind profile obtained using (a) the current urban scheme and (b) the new urban scheme in TAPM at Site Ue3 during 4e5 July.

56

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

Fig. 14. Scatter plot of the observed vs. TAPM simulated concentrations: (a) when the current urban scheme is used and (b) when the new urban scheme is used.

Table 2 Performance statistics for concentration (see start of Section 4 for definition of statistics, sample size ¼ 153). Scheme

d

r

NMSE

FB

FAC2 (%)

Current New Current-assim. New-assim. Rotach et al., 2005

0.73 0.79 0.78 0.82 0.48

0.69 0.64 0.77 0.74 0.26

4.65 2.48 2.40 1.51 4.21

0.52 0.07 0.72 0.07 0.01

23 33 38 44 38

concentrations for all experiments we used their scatter plot Fig. 21 to derive the performance statistics presented in Table 2. Comparing these to the corresponding TAPM statistics suggests that, overall, TAPM’s performance in Fig. 15b is significantly better despite TAPM’s use of a simpler Lagrangian dispersion scheme driven by modelled turbulence levels. We have not investigated what the causes for this performance difference are, but factors such as proper inclusion of horizontal advection in TAPM may be responsible for its better performance. Rotach et al. (2004) state that the likely reasons for the non-ideal performance of their particle dispersion model driven by observed turbulence are their parameterisation of the TKE dissipation rate and the nonaccounting of potential recirculation and trapping of plume particles between the roughness elements. The TKE dissipation rate is modelled in TAPM, whereas TAPM too cannot simulate the recirculation and trapping properly unless the individual buildings are

resolved. In the future, it will be of interest to examine the TAPM predicted turbulence levels vis-à-vis observations and their spatial variability. Quantileequantile (qeq) plots, in which sorted predicted concentrations are plotted against sorted observed values (i.e. independent of time and position), are often used in order to examine any model bias over the concentration distribution. In the qeq plot in Fig. 16a, the current urban scheme (without wind assimilation) overpredicts for observed concentrations above 70 ng m3, whereas the new urban scheme performs noticeably better. There are no significant differences between the two schemes for the observed concentrations below 70 ng m3, with both schemes underestimating for concentrations below 30 ng m3. In Fig. 16b with wind assimilation, the previous model underestimation for the lower-end concentrations has improved for both schemes, but the current scheme overpredicts mid- and higher-end concentrations even more. As indicated earlier, the tracer measurements were taken in daytime when the relative differences between the predicted sensible heat fluxes from the two schemes are not expected to be as large as those at night. This is likely to be due to the large-scale convective mixing resulting from the solar heating of the ground, which dominates dispersion instead of the thermal effects of the canyon geometry. To investigate concentration differences between the two schemes on diurnal basis, tracer was released in TAPM as per Table 1, but with a continuous emission over the full diurnal period (instead of the actual release periods). No wind data

Fig. 15. Scatter plot of the observed concentrations versus the concentrations simulated by TAPM with wind data assimilation: (a) when the current urban scheme is used and (b) when the new urban scheme is used.

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

57

Fig. 16. Quantileequantile plot of the modelled versus observed concentration: (a) without and (b) with wind data assimilation in TAPM.

assimilation was used. The concentrations predicted at all the SF6 sampling locations were averaged for each diurnal hour. These concentration averages scaled by the emission rate (Q) are plotted in Fig. 17 as function of time of day. Considering 13:00e18:00 h, which covers the tracer release periods, it is apparent that the differences between the concentration averages from the two schemes are relatively small for 13:00e15:00 h, but they progressively increase with time and by 18:00 h the current scheme yields a concentration average that is almost 4 times as high as that from the new scheme. Beyond 18:00e21:00 h, this difference is even greater. For 22:00e10:00 h, the difference is in the opposite direction where the new scheme results in higher concentration averages. The above demonstrates diurnal differences between the two schemes. It is difficult to pinpoint the exact parameters that are responsible for these, but generally speaking these would be the stability and turbulence properties, and also the potential for the winds, particularly the wind direction, being affected by the choice of the urban scheme. 5. Conclusions We evaluated a recently developed urban surface scheme in the operational mesoscale model TAPM for flow and tracer dispersion using data from the 2002 BUBBLE experiment. This scheme is based on a building-averaged town energy balance (TEB) approach with a generic canyon geometry, and improves upon the current scheme based on a simple slab approach. It is observed that the new

scheme leads to an overall improvement in the prediction of temperature and surface fluxes: the simulated average sensible heat flux is always positive at night, in agreement with the data, and the modelled magnitude of the average latent heat flux is comparable with the observed one. The new scheme is able to reproduce the observed near-neutral to weakly unstable conditions at night, which is a feature of urban meteorology. In contrast, the current scheme predicts stable conditions at night. These differences in the modelled stability have ramifications for dispersion predictions. TAPM was also used to simulate the tracer releases over four separate days, and it is found that the observed concentration fields are better predicted using the new urban scheme, but because there were no nighttime observations, the full capability of the new scheme under all diurnal conditions could not be demonstrated. For the applications considered here, the computational efficiency of the new scheme in TAPM is on par with the current scheme. The new urban scheme is building-averaged, and therefore does not account for any local variations in the building geometry (e.g. the canyon orientation not being uniformly distributed or a surface that cannot be represented by the generic canyon geometry considered) that may influence the fluxes and other parameters. However, it offers a significant improvement in the representation of the urban canopy over the slab approach without requiring considerable additional computational resources, thus making it feasible to apply it for carrying out long-term simulations (e.g. at climate scale). Coupling surface fluxes to the atmosphere via surface similarity requires more work since the roughness sublayer is not explicitly

Fig. 17. Diurnal variation of the scaled modelled concentrations averaged over all the SF6 sampling locations (filled squares) computed using the (a) current urban scheme, and (b) new urban scheme. The vertical bars correspond to one standard deviation.

58

A.K. Luhar et al. / Atmospheric Environment 88 (2014) 47e58

included. In that regard, additional experiments involving a full resolution of the roughness sublayer and its transition to the inertial sublayer will be very valuable. Flow resolving models based on computational fluid dynamics and large-eddy simulation (e.g., Baik et al., 2003; Nakayama et al., 2011) can prove useful in gaining a better understanding of urban processes and for developing parameterisations for use in operational models.

Acknowledgements We thank the BUBBLE Program, in particular R. Vogt, A. Christen, S.-E. Gryning and M.W. Rotach, for the data used in this work. The two anonymous reviewers provided useful comments.

References Batchvarova, E., Gryning, S.-E., 2006. Progress in urban dispersion studies. Theoretical and Applied Climatology 84, 57e67. Baik, J.-J., Kim, J.-J., Fernando, H.J.S., 2003. A CFD model for simulating urban flow and dispersion. Journal of Applied Meteorology 42, 1636e1648. Bornstein, R., Lin, Q., 2000. Urban heat islands and summertime convective thunderstorms in Atlanta: three case studies. Atmospheric Environment 34, 507e 516. Christen, A., Rotach, M.W., 2004. Estimating wind speed at an urban reference height. In: Fifth AMS Symposium on the Urban Environment, 23e27 August, Vancouver, BC, Canada. Christen, A., Vogt, R., 2004. Energy and radiation balance of a central European city. International Journal of Climatology 24, 1395e1421. Grimmond, C., Oke, T., 1995. Comparison of heat fluxes from summertime observations in the suburbs of four North American cities. Journal of Applied Meteorology 34, 873e889. Grimmond, C., Blackett, M., Best, M., Baik, J., Belcher, S., Beringer, J., Bohnenstengel, S., Calmet, I., Chen, F., Coutts, A., Dandou, A., Fortuniak, K., Gouvea, M., Hamdi, R., Hendry, M., Kanda, M., Kawai, T., Kawamoto, Y., Kondo, H., Krayenhoff, E., Lee, S., Loridan, T., Martilli, A., Masson, V., Miao, S., Oleson, K., Ooka, R., Pigeon, G., Porson, A., Ryu, Y., Salamanca, F., Steeneveld, G., Tombrou, M., Voogt, J., Young, D., Zhang, N., 2011. Initial results from phase 2 of the international urban energy balance model comparison. International Journal of Climatology 31, 244e272. Gryning, S.-E., Batchvarova, E., Rotach, M., Christen, A., Vogt, R., 2005. Roof-Level SF6 Tracer Experiments in the City of Basel. Institute for Atmospheric and Climate Science ETH, Zürich, p. 88. Report 83. Hamdi, R., Masson, V., 2008. Inclusion of a drag approach in the Town Energy Balance (TEB) scheme: offline 1D evaluation in a street canyon. Journal of Applied Meteorology and Climatology 47, 2627e2644. Hamdi, R., Schayes, G., 2007. Validation of Martilli’s urban boundary layer scheme with measurements from two mid-latitude European cities. Atmospheric Chemistry and Physics 7, 4513e4526. Harman, I.N., Barlow, J.F., Belcher, S.E., 2004. Scalar fluxes from urban street canyons. Part II: model. Boundary-Layer Meteorology 113, 387e409. Hurley, P., 1994. PARTPUFF e a Lagrangian particle/puff approach for plume dispersion modelling applications. Journal of Applied Meteorology 33, 285e 294. Hurley, P., 2008. TAPM V4. Part 1: Technical Description. CSIRO Marine and Atmospheric Research Paper No. 25, Australia, p. 59. Available at: http://www. cmar.csiro.au/research/tapm/docs/tapm_v4_technical_paper_part1.pdf. Hurley, P., Luhar, A., 2009. Modelling the meteorology at the Cabauw Tower for 2005. Boundary-Layer Meteorology 132, 43e57. Hurley, P.J., Physick, W.L., Luhar, A.K., 2005. TAPM: a practical approach to prognostic meteorological and air pollution modelling. Environmental Modelling and Software 20, 737e752.

Kowalczyk, E.A., Garrat, J.R., Krummel, P.B., 1994. Implementation of a Soil-canopy Scheme into the CSIRO GCM e Regional Aspects of the Model Response, p. 65. CSIRO Division of Atmospheric Research Technical Paper 32. Lemonsu, A., Masson, V., 2002. Simulation of a summer urban breeze over Paris. Boundary-Layer Meteorology 104, 463e490. Lemonsu, A., Grimmond, C.S.B., Masson, V., 2004. Modeling the surface energy balance of the core of an old Mediterranean city: Marseille. Journal of Applied Meteorology 43, 312e327. Luhar, A.K., Hurley, P., 2003. Evaluation of TAPM, a prognostic meteorological and air pollution model, using urban and rural point source data. Atmospheric Environment 37, 2795e2810. Luhar, A.K., Venkatram, A., Lee, S.-M., 2006. On relationships between urban and rural near-surface meteorology for diffusion applications. Atmospheric Environment 40, 6541e6553. Luhar, A.K., Mitchell, R.M., Meyer, C.P., Qin, Y., Campbell, S., Gras, J.L., Parry, D., 2008. Biomass burning emissions over northern Australia constrained by aerosol measurements: IIeModel validation, and impacts on air quality and radiative forcing. Atmospheric Environment 42, 1647e1664. Luhar, A.K., Hurley, P.J., 2012. Application of a coupled prognostic model to turbulence and dispersion in light-wind stable conditions, with an analytical correction to vertically resolve concentrations near the surface. Atmospheric Environment 51, 56e66. Martilli, A., 2003. A two-dimensional numerical study of the impact of a city on atmospheric circulation and pollutant dispersion in a coastal environment. Boundary-Layer Meteorology 108, 91e119. Martilli, A., Clappier, A., Rotach, M.W., 2002. An urban surface exchange parameterization for mesoscale models. Boundary-Layer Meteorology 104, 261e304. Masson, V., 2000. A physically-based scheme for the urban energy budget in atmospheric models. Boundary-Layer Meteorology 94, 357e397. Masson, V., 2006. Urban surface modeling and the meso-scale impact of cities. Theoretical and Applied Climatology 84, 35e45. Masson, V., Grimmond, C.S.B., Oke, T.R., 2002. Evaluation of the Town Energy Balance (TEB) scheme with direct measurements from dry districts in two cities. Journal of Applied Meteorology 41, 1011e1026. Masson, V., Seity, Y., 2009. Including atmospheric layers in vegetation and urban offline surface schemes. Journal of Applied Meteorology and Climatology 48, 1377e1397. Nakayama, H., Takemi, T., Nagai, H., 2011. LES analysis of the aerodynamic surface properties for turbulent flows over building arrays with various geometries. Journal of Applied Meteorology and Climatology 50 (8), 1692e1712. Oke, T.R., 1982. The energetic basis of the urban heat island. Quarterly Journal of the Royal Meteorological 108, 1e24. Oke, T.R., 1988. The urban energy balance. Progress in Physical Geography 12, 471e 508. Otte, T.L., Lacser, A., Dupont, S., Ching, Jason K.S., 2004. Implementation of an urban canopy parameterization in a mesoscale meteorological model. Journal of Applied Meteorology 43, 1648e1665. Pielke, R.A., 2002. Mesoscale Meteorological Modelling. Academic Press, p. 676. Rotach, M.W., Gryning, S.-E., Batchvarova, E., Christen, A., Vogt, R., 2004. Pollutant dispersion close to an urban surface e the BUBBLE tracer experiment. Meteorology and Atmospheric Physics 87, 39e56. Rotach, M.W., Vogt, R., Bernhofer, C., Batchvarova, E., Christen, A., Clappier, A., Feddersen, B., Gryning, S.-E., Martucci, G., Mayer, H., Mitev, V., Oke, T.R., Parlow, E., Richner, H., Roth, M., Roulet, Y.-A., Ruffieux, D., Salmond, J.A., Schatzmann, M., Voogt, J.A., 2005. BUBBLE e an urban boundary layer meteorology project. Theoretical and Applied Climatology 81, 231e261. Roulet, Y.-A., Martilli, A., Rotach, M.W., Clappier, A., 2005. Validation of an urban surface exchange parameterization for mesoscale modelsd1D case in a street canyon. Journal of Applied Meteorology 44, 1484e1498. Stauffer, D.R., Seaman, N.L., 1994. Multiscale four-dimensional data assimilation. Journal of Applied Meteorology 33, 416e434. Thatcher, M., Hurley, P., 2012. Simulating Australian urban climate in a mesoscale atmospheric numerical model. Boundary-Layer Meteorology 142, 149e175. Willmott, C.J., 1981. On the validation of models. Physical Geography 2, 184e194. Zawar-Reza, P., Sturman, A., 2008. Application of airshed modelling to the implementation of the New Zealand National Environmental Standards for air quality. Atmospheric Environment 42, 8785e8794.