On the application and grid-size sensitivity of the urban dispersion model CAIRDIO for the simulation of a real city and realistic meteorology

There is a gap between the need for city-wide air-quality simulations considering the intra-urban variability and mircoscale dispersion features and the computational capacities that conventional urban microscale models require. This gap can be bridged by targeting model applications on the gray zone situated between the mesoscale and large-eddy scale. The urban dispersion model CAIRDIO is a new contribution to the class of computational-fluid dynamics models operating in this scale range. It uses a diffuse-obstacle boundary method to represent buildings as physical obstacles at gray-zone resolutions in 5 the order of tens of meters. The main objective of this approach is to find an acceptable compromise between computationally inexpensive grid sizes for spatially comprehensive applications and the required accuracy in the description of building and boundary-layer effects. For this purpose, CAIRDIO is applied in dispersion simulation of black carbon and particulate matter for an entire mid-size city using an uniform horizontal resolution of 40m in this paper. For evaluation, the simulation results are compared with measurements from 5 operational air monitoring stations, which are representative for the urban background 10 and high-traffic roads, respectively. Moreover, the comparison includes the mesoscale host simulation, which provides the boundary conditions. The temporal variability of the concentration measurements at the background sites was largely influenced only by the characteristics of the mixing layer. As a consequence, the model results were not significantly dependent on spatial resolution, so that the mesoscale simulation also performed reasonably well. At the traffic sites, however, concentrations were in addition markedly influenced by the proximity to road-traffic sources and the surrounding building environment. Here, 15 the mesoscale simulation indiscriminately reproduced almost the same urban-background profiles, which resulted in a large positive model bias. On the other hand, the CAIRDIO simulation was able to respond to the significantly amplified diurnal variability with its pronounced rush-hour peaks. This resulted in a consistent improvement of the model deviation to measurements compared to the mesoscale simulation. Nevertheless, discrepancies to measurements remain in the 40m-CAIRDIO simulation, e.g., an underestimation of peak concentrations at two traffic sites inside narrow street canyons. To further research 20 resolution sensitivity, the horizontal grid spacing of locally nested CAIRDIO domains is refined down to 5m. While for the street canyons the representation of peak concentrations can be improved using horizontal grid spacings of up to 10m, no further improvements beyond this resolution can be observed. This suggests that the too low peak concentrations with the default grid spacing of 40m result from an inadequate representation of the traffic emissions inside narrow street canyons. If the total gain in accuracy due to the grid refinements is put in relation to the remaining model error, the improvements are only modest. 25 1 https://doi.org/10.5194/gmd-2021-358 Preprint. Discussion started: 22 November 2021 c © Author(s) 2021. CC BY 4.0 License.

can be further modified locally in the model by a parameterized surface-heat flux from ground and building surfaces. Inflow conditions are in general not only turbulent but also transient, in order to account for an accurate evolution of the larger-scale 95 meteorology. The complexity of the simulation is further increased by using a comprehensive emission inventory that includes all relevant sectors, which are modulated in time to account for diurnal and weekly changes in activity. While this study aims not at analyzing individual processes in depth, its main objectives are to demonstrate the feasibility and practicability of the approach as a downscaling tool for a more accurate representation of the intra-urban air-pollution variability. Therefore, apart from static inputs, the model solely relies on the output fields of a host simulation conducted at the lower end of the 100 mesoscale, for which the CTM COSMO- MUSCAT (Wolke et al., 2012) is used. For validation, we compare modeled PM 10 and/or BC concentrations with measurements at 5 different operational air-monitoring sites in Leipzig for a total period of two consecutive days in spring 2020. To further estimate the sensitivity to the horizontal resolution, locally-nested sensitivity runs are performed, for which the horizontal resolution is increased from a default ∆h = 40 m in steps down to ∆h = 5 m, enabling conventional building-resolved simulations. 105 The paper is structured as follows: Section 2 describes the methodology, in which all the general and technical aspects of the simulations and measurements are described. This also includes a detailed description of the mesoscale coupling. Section 3 includes the presentation and discussion of model results, which is subdivided into a part describing the modeled PBL evolution, a model evaluation with concentration measurements, including a comparison with results from the CTM COSMO-MUSCAT, and the grid-size sensitivity study. Section 4 summarizes the main findings of the study and highlights the advantages but also 110 limitations of the demonstrated approach and the study itself.

Study time period
The model case study spans 2 consecutive days from 1 March 2020, 00:00 UTC to 3 March 2020, 00:00 UTC, to address the main objectives of this study. Yet, for a more significant model evaluation with observational data, a substantially longer 115 simulation period needs to be simulated. While principally our approach is computationally much cheaper compared to a wellresolved urban microscale simulation, a compromise still had to be found, and we decided to invest computation resources in a spatially more comprehensive simulation that fits better the aspect of a model case study. The specific simulation period was selected based on suitable properties for an investigation of the intra-urban air-pollution variability. Firstly, quality-assured observational data from all operational air-monitoring sites in Leipzig are available during this period. Secondly, significant 120 impacts of the world-wide spreading Covid-19 pandemic had still not reached the German public by early March 2020, as data from the Google mobility report show a significant decline starting after 10 March (Google-LLC, 2020;Forster et al., 2020). This provided confidence that the traffic emissions had not to be adjusted for a reduced mobility, which is a potential source of additional uncertainty. Thirdly and most importantly, the meteorological conditions were suitable to focus on local air pollution and PBL processes affecting it. The large-scale synoptic pattern during the simulation period from 1 to 3 March 125 2020 was dominated by a large and deep low-pressure system situated over northwestern Europe, from which two troughs 4 https://doi.org/10.5194/gmd-2021-358 Preprint. Discussion started: 22 November 2021 c Author(s) 2021. CC BY 4.0 License. protruded southward and eventually moved across Germany (see Figure 1a-c). The associated unsettled weather conditions in Leipzig resulted in diverse PBL characteristics and effects on local air quality, which are interesting to study. Moreover the influence of low pressure favored low background near-surface PM 10 concentrations over most of Germany, as suggested by results from an air-quality simulation for Europe depicted in Fig. 1d-f. According to this, highest PM 10 concentrations 130 apart from the well-known air-polluted regions, like the Po Valley, occurred over the eastern half of Europe. There were also periods before and after the actual simulation period, when the Siberian high pressure system extended westward and brought a polluted continental air mass to central Europe (not shown). During such periods with elevated background concentrations, the intra-urban air-pollution variability was quite insignificant and not worth to study.

Observations
For model evaluation a range of in-situ air quality observations is used. The south-eastern German state Saxony operates 26 air-quality monitoring sites. Three of them are suitable for our model evaluation, as they are located within the city margins of Leipzig and provide PM 10 concentration measurements. One of the three stations considered is also co-operated by the TROPOS and provides BC retrievals too. Three stations are air quality stations operated by the Saxon State office while two additional stations belong to TROPOS. Three out of the five station measure also equivalent black carbon (eBC) mass con-140 centration, which are also relevant to our evaluation. Tab. 1 lists some basic information of the aforementioned measurement sites, while Fig. 2 shows their locations on horizontal maps of the simulation domains to provide a qualitative overview of the characteristic environment, like the distribution of buildings, parks and major roads highlighted by PM 10 line emissions.
Finally, Fig. 3 gives a detailed mapping of the surrounding building environment in spherical coordinates at each exact measurement position, as well as the corresponding sky-view factors f sky . Based on these information, a brief introduction of 145 each measurement site is given in the following. Station Leipzig Lindenau (LL) measures PM 10 and is located in the western city district Leipzig Lindenau within a closed street canyon running in WNW-ESE direction. The street canyon has average dimensions (height × width) of roughly 20 m × 25 m, which results in a sky-view factor of 0.34. The horizontal position of the measurement site is near the northern side of the canyon at a distance of only 5 m to a traffic-busy road (average daily traffic count (ADTC) of 20400 vehicles). This, in combination with an inlet height of 1.7 m leads to a high exposure of this station 150 to the exhaust gases from nearby traffic, which clearly classifies this station as roadside / street canyon. Station Leipzig Mitte (Leipzig Center, LC) provides PM 10 and eBC measurements. It is located southwest of the Leipzig main railway station at a junction of the inner-city ring, which is a multi-line road (ADTC of 47600 vehicles). Therefore, also this station is classified as roadside. Compared to the site LL, there is more open area around the station LC (f sky = 0.79), with the closest significant building (height times width of 27 m × 50 m) being to the south of the station at a distance of approximately 45 m. Further-155 more, due to a local park adjoining to the east, the influence of traffic emissions at the measurement site can be expected to vary according to the prevailing wind direction. Station Leipzig West (LW), located within the western outskirts of Leipzig inside a park, is a background station for PM 10 , as it is also secluded by lines of trees from a nearby road (ADTC of 8600 vehicles). Station Leipzig Eisenbahnstr. (LE) has a long history as a scientific measurement site and is thus well documented from previous air-quality studies (see e.g. Klose et al. 2009, Wiesner et al. 2021. The measurement equipment is located next to a window on 160 the third floor of an apartment house (inlet height is approx. 6 m above the road) flanking a frequently traffic-congested street canyon (ADTC of 11500 vehicles). The cross section of the street canyon is symmetric (20 m × 20 m). Regularly occurring crossroads divide the street canyon into segments of 70 m−110 m length, with the closest crossroad (ADTC of 11800 vehicles) being to the west at a horizontal distance of about 35 m from the measurement site. While this side is also classified as roadside, its inlet position high above the road makes it more representative to the average concentrations within the street canyon.

165
Depending on the development of the street-canyon vortex, however, it may be also more directly exposed to high pollution concentrations. Finally, the side Leipzig TROPOS (LT) is a background station for eBC, as it is located on the roof top of the TROPOS institute's building at a height of 16 m and at a distance of at least 100 m from any busy roads.
PM 10 concentrations are directly and near-continuously measured using the tapered element oscillating microbalance (TEOM) system (scientific ambient particulate monitor TEOM 1405, Thermo Fisher Scientific Inc.). TEOM derives PM mass concen-170 trations from the frequency-change of an oscillating hollow tube caused by deposited material at one end of the tube (Page et al., 2007). Real-time measurements are averaged to hourly-mean values with a stated precision of ±2.0 µg and an accuracy of 0.75 % (TFS, 2019). eBC is indirectly retrieved from optical principles with multi-angle absorption photometers (MAAP 5012,Thermo Fisher Scientific Inc.). MAAP estimated the absorption coefficient of an aerosol probe from the transmission and back-scattering of light at a wavelength of 637 nm, where eBC is the main absorber (Petzold and Schönlinner, 2004). The 175 eBC mass concentrations calculated with a mass absorption cross section of 6.6 m 2 g −1 are assumed to be directly comparable with modeled BC mass concentrations, and have an uncertainty between 5 % and 12 % according to different sources (Wiesner et al., 2021).  Figure 2. Map of the city area of Leipzig, which is also selected as the coarse-grid CAIRDIO simulation domain L0 introduced in Section 2.4.1. Each of the white boxes contains an operational air monitoring site used for model evaluation. In addition, a magnified view of the area within each box shows the local environment around the corresponding air-monitoring site, which is highlighted by a red circle. These areas also correspond to the CAIRDIO subdomains introduced in Section 2.5. Traffic-PPM10 emissions of major roads are represented by line sources. 1997; Taraborrelli et al., 2009). The remaining fraction of PM 10 is primarily emitted (PPM), and approximated by chemically inert tracers that are only subjected to physical atmospheric removal processes. Figure 4 gives an overview of the bulk PM 10 decomposition in MUSCAT as it is available from the model output. COSMO-MUSCAT is applied on a hierarchy of refined 190 domains, with a one-way nesting technique providing the boundary condition for each consecutive simulation (see Figure 5 for an overview of the domains). This model setup has already been used to provide quasi-operational air-quality forecasts for the citizen-science campaign WTImpact (Heinold et al., 2019;Tõnisson et al., 2021). The outermost domain M0 has a spatial resolution of 14 km and covers entire Europe. This domain is initialized and driven by meteorological data from the global model ICON (Zängl et al., 2015) operationally run by the German weather service (Deutscher Wetterdienst, DWD).

195
Initialization and boundary condition for air chemistry are interpolated from operational air-quality forecasts with the model system ECMWF IFS (Copernicus Atmosphere Monitoring Service) (Flemming et al., 2014). The domain M0 is simulated for an extended period in time (at least two weeks) ahead of the actual simulation period of two days. This allows for a proper based on wind-tunnel data also presented in Weger et al. (2021a). For DOB, the same geometric building dataset as already used for deriving the parameters in DCEP is used. At the horizontal scale of 40 m, buildings are effectively represented as diffuse obstacles. As a consequence, the subgrid-scale mixing length is of similar magnitude or even larger than the average space between buildings. This implies that the simulation is mostly non-eddy resolving (similar to a RANS approach) within 230 the urban canopy. For a more accurate representation of vertical gradients and mixing processes near the surface, the vertical dimension is kept much finer resolved, with 5 m grid spacing within the urban canopy. Increased grid stretching is applied above, with a maximum stretching factor of 4 near the domain top. Note that vertical resolution is not as computationally expensive as horizontal resolution, because the Courant-Friedrichs-Lewy (CFL) criterion, which limits the explicit integration time step, is mostly determined by the horizontal wind speed. The horizontal coverage of domain L0 is roughly 18 km × 18 km, 235 which is extensive enough to not only accommodate the complete city of Leipzig but also a sufficient fetch to allow for a relaxation of the lateral boundary conditions. The domain top is located at about 350 m height, which is generally below the vertical extend of a typical convective PBL. However, a suitable boundary condition allows for vertical motions exiting/entering at the domain top, which is further explained in the subsequent paragraph 2.4.2 addressing the mesoscale forcing. While the surface orography within domain L0 is not mountainous, subtle effects from it can still influence meteorology. CAIRDIO can 240 be used with terrain-following coordinates, which in this simulation are inferred from surface elevation data (DGM1) provided by the Staatsbetrieb für Geobasisinformation und Vermessung Sachsen (GeoSN). This dataset with a spatial resolution of 1 m has also been corrected for vegetation and buildings and is thus compatible with the explicit representation of buildings by DOB. Surfaces fluxes of momentum, heat and moisture from vegetation, as well as other types of land cover (lakes, bare soil, subgrid-scale structure of buildings) are parameterized using Monin-Obukhov similarity theory. Therein, each surface type is characterized by a parametric roughness length z 0 . Table 2 lists the z 0 values related to each land-cover class used in the model.
Land-cover is based on a combination of the Pan-European land cover map for 2015 with 30 m spatial resolution (Pflugmacher et al., 2018), and the more detailed land-cover map by Banzhaf and Kollai (2018) (better than 5 m resolution) for most of the urban area. The combined dataset is depicted in Fig. 2 for domain L0, as well as for the finer resolved nested subdomains introduced in the paragraph 2.5 addressing the LES-to-LES nesting.

250
The emissions used for domain L0 are on the same basis as the ones already used for mesoscale domain M3, but are locally refined were applicable. Road and railway emissions (SNAP 7 and partly SNAP 8) are available as line sources and can in principle be gridded on arbitrary grids. In our setups, both the building geometries and line emissions are geometrically consistent to each other. Other sources from industry (SNAP 3, SNAP 4, SNAP 6, SNAP 9) and residential combustion (SNAP 2), as well as non-classified mobile sources (SNAP 8) are all originally given on a grid with 500 m spacing. Conveniently, 255 the metadata of the building-geometry dataset provide additional purpose-based classification into residential and commercial buildings. This permits to relocate the corresponding area sources at the building roofs (whereas the building-volumes are taken as weights) and a subsequent gridding of the building emissions on the finer L0 grid. The advantage of this emission downscaling is currently investigated in a separate sensitivity study. The static emission fields are modulated with a time profile also taken from the TNO-MACC2 database. The spatially integrated and temporally modulated emissions are shown in Fig. 6.

260
Accordingly, the traffic sector is by far the most important contributor to BC emissions during the simulation period.

Mesoscale forcings and boundary conditions
Simulation results from COSMO-MUSCAT mesoscale domain M3 are used to drive the meteorological and air-pollution fields of the city-scale domain L0. Initial and boundary conditions for the meteorological prognostic fields, which include the 3-D wind components, potential temperature Θ, specific humidity Q v , and subgrid-scale TKE E sgs , are spatially interpolated using 265 tricubic interpolation. For the air-pollution fields instead, trilinier interpolation is used. The 3-D interpolation procedure is carried out as a sequence of 2-D horizontal interpolation followed by vertical interpolation. For the horizontal interpolation, the Climate Data Operators (CDO) software (Schulzweida, 2019) is used, which is convenient to remap data from rotated lat/lon coordinates of the COSMO-MUSCAT model directly to Lambertian azimuthal equal-area coordinates (epsg:3035) of the L0 grid. Vertical interpolation is based on the 3-D height of half levels (HHL), which coincide with the locations of vertical 270 velocity of the staggered grid. After horizontal remapping of the HHL field of the M3 grid, the vertical interpolation weights are generated by computing Lagrange polynomials of the desired accuracy from the HHL data. As the CAIRDIO simulation employs a finer grid spacing near the ground than COSMO-MUSCAT, the first vertical levels need to be extrapolated, for which a level with zero-height is introduced. At this zero-height level, all wind components are set to zero, and the potential temperature as well as specific humidity fields assume the respective surface values Θ S and Q S v . For the air-pollution fields, for the inertial subrange E(k) ∝ k −5/3 up to the different cut-off wave-numbers k min of the fine and coarse grids. k min can be directly related to the subgrid-mixing scale ∆ sgs , then the following expression follows: (1) ∆ f sgs can be crudely approximated by twice the horizontal grid spacing (corresponding to the Nyquist wavenumber). Note that the horizontal grid spacing in typical PBL simulations is equal or larger than the vertical grid spacing and is thus the 285 dominant cut-off scale. ∆ c sgs , on the other hand, can be related to the master-mixing length of the mesoscale simulation. E f sgs is finally the lateral boundary condition for the prognostic subgrid-scale TKE equation, which determines the eddy diffusivities.
The lateral boundary condition for the 3-D wind vector is a Dirichlet/radiation condition that can flexibly distinguish inflow from outflow regions. For inflow regions, a superposition of the interpolated mesoscale wind field and recycled turbulence is prescribed. The scale-separation applies in a similar way for the turbulence recycling scheme, as the cut-off wavelength of 290 the extraction filter is chosen similar to ∆ c sgs . Consequently, the recycled turbulent fluctuations are scaled with the resolvable energy part E f rec . Convective PBLs that cannot be simulated for their entire vertical extend require a boundary condition that allows for turbulent motions transcending the top-domain boundary, while the mesoscale forcing shall still apply at the same time. A mixed Dirichlet/Neumann condition addressing both scale ranges independently in principle can satisfy these requirements. In practical simulations, however, the small-scale fluctuating Neumann contribution eventually tended to grow 295 exponentially, which required a rescaling with E f rec and a limiting of the maximum amplitude of vertical motions for numerical stabilization.
Besides the initial and boundary conditions for the prognostic fields, another important mesoscale forcing includes the coupling with the land surface, as surface potential temperature Θ S and surface specific humidity Q S v are not computed in CAIRDIO. In this study, a land-cover based fitting method is favored over a more simple polynomial interpolation of the corresponding 300 mesoscale surface fields. In doing so, much of the spatial detail can be preserved by assuming that land-cover is the most significant co-determinant of the small scale variability. This may only be valid for small horizontal domain extents and may require the application of an additional bilinear, or higher-order correction if more significant large-scale variations are present.
In the following, the unknown potential temperature fluctuations Θ S Lc of each considered land-cover class from a spatially filtered state Θ S are obtained by solving a least-square problem using the mesoscale data fields from the COSMO model: with Θ S = Θ S − Θ S . Note that the large scale variations not related to different land-cover classes have to be excluded in Θ S at this point. The matrix L is of shape m × n, and contains the fractional land cover f Lc for the total number of landcover classes m with unknown potential temperature and number of grid cells n of domain M3. The vector b contains the potential temperature contribution from a-priori determined land-cover classes. A simplification leading to more robust fitting as independent classes. Already determined are the surface temperature and specific humidity values of water surfaces, as these values are spatially constant in the mesoscale simulation. Additionally, the urban parameterization DCEP in COSMO provides potential temperature estimates for impervious ground surfaces and roof surfaces at building sites, which both are included in b. After solving Eq. 2, the obtained land-cover dependent fluctuations Θ S Lc are added to the interpolated filtered 315 state Θ S to compose the new surface potential temperature field on domain L0 using the much higher resolved land-cover data. Because CAIRDIO uses a 3-D building structure, potential temperature values from these elevated horizontal and vertical surfaces are still missing. These values can, however, directly be interpolated from the additional output fields of the mean roof and wall temperature computed with the DCEP parameterization. Q S V is fitted in a similar fashion to Θ S , with a further simplification that Q S V is set to zero for the impervious surface fraction. Thus, rain evaporation is neglected. Fig. 7 demonstrates 320 the described downscaling approach based on an exemplary Θ S field for 2 March 2020, 12 : 00 UTC, when sunshine prevailed.
For the resulting reconstructed field of domain L0, a top-down projection is shown. The sun-lit roofs are the warmest surfaces with quite a spatially uniform distribution due to a constant prescribed roof albedo of 0.16 in DCEP. Considerably cooler are the ground surfaces inside the partly shaded street canyons, followed by the forest areas, and finally by the seasonally cold lakes with the lowest surface potential temperature. While in the given example, surface orography is mostly flat and the surface 325 elevation could be neglected as an additional disturbing factor, this approach may be an oversimplification for areas with more mountainous terrain. For the resolved part E f res , which scales the inserted turbulent fluctuations at the inflow boundaries, a slightly modified formula is used in order to consider the missing contribution of numerical diffusion in E c sgs :

M3 L0
where the extraction-filter width is set to λ cut = 140m, which is considerably larger than ∆ c sgs = 80 m. Note that the use of 340 an exponential filter function results in a smooth cut-off range, with λ cut defined as the wavelength that is scaled e −1 -fold.
In Figure 8 the described nesting method with the energy-scale separation is demonstrated by an example with domain LE and 5 m spatial resolution, which features a stably-stratified, shear-driven PBL. The plots of the dominant velocity component v ( Fig. 8a and b) show that well-developed turbulence already exists near the southern inflow boundary, which qualitatively matches the turbulence more distant from the boundary well. This is also quantitatively shown by the derived energy spectra 345 shown for the x-dimension (Fig. 8c), which do not evolve much when moving further away from the inflow boundary. Note that the inertial subrange is followed by the dissipation range, which can be attributed to the combined (dissipative and dispersive) numerical error of the advection scheme at significant convective speeds (Yalla et al., 2021).
The surface fields Θ S and Q S v are again obtained from the corresponding mesoscale fields using the land-cover based method described in Sect. 2.4.2. An additional scaling is applied, such that the computed horizontally averaged values are independent 350 of the spatial resolution and also correspond to the average values of the congruent part of the parent domain.

Synoptic and urban planetary-boundary layer meteorology
In the following, an overview of the meteorological conditions and resulting PBL characteristics during the simulation period from 1 March 2020, 00:00 UTC to 3 March 2020, 00:00 UTC is given. The implications of the variable PBL structure on the 355 modeled concentrations and transport of local air pollution are then qualitatively discussed based on model outputs from the CAIRDIO domain L0. As already briefly mentioned in Section 2, the weather in Leipzig during the simulation period was In Figure 11, qualitative model results for two contrasting PBL states are featured as simulated with the CAIRDIO domain 380 L0. In the first case (Fig. 11a), the dominant horizontal wind component v near the surface is shown for the stably-stratified, shear-driven PBL on 2 March at 00:30 UTC. Turbulence is organized into horizontal streaks near the surface as it is typical for a shear-driven PBL. The effect of significant surface roughness from the city structure and forest patches locally reduced the wind speed. Over these areas, turbulence was also of more intermittent nature due to the reduced vertical wind shear near the surface and limited turbulent energy production. The largest buildings, e.g., from factories can already be resolved with 385 40 m resolution, while the model representation of the air flow through the diffuse city structure consisting of smaller building units is more typical of a porous-media flow (see corresponding insets). Clearly, building-induced turbulence within the diffuse urban canopy cannot be resolved using such a still comparatively coarse grid spacing, but its mixing effects are represented by diffusion from the subgrid-scale model. The second case (Fig. 11b)  also indicates a significant model sensitivity to the mixing parameterization (e.g., the prescription of the static mixing length).

Quantitative model comparison with air-monitoring measurements
To quantitatively evaluate the model representation of the intra-urban variability of air pollution, hourly averaged model output evening are underestimated. Correspondingly, the model bias is only slightly positive (FB = 0.14). Finally, the measured PM 10 profile at site LC (Fig. 13f) is more comparable to the measured background at site LW, which is a bit surprising given that the station is classified as roadside. In anyway, measured PM 10 concentrations seem to be more influenced by other not well-known where the brackets <> h denote for the horizontal averaging and c * i is the corrected snapshot. Note that the horizontal average within the canopy layer excludes the inaccessible grid-cell volume by using the volume-scaling field χ as weights. For 485 the velocity components, the cell-face area scaling field η is used instead of χ. Vertical profiles of the computed horizontally and hourly averaged variables u, v, Θ v , c BC , u w , v w , and c BC w are depicted in Figure 14 for the two contrasting PBL states already discussed in Section 3.1. Starting with the first case on 2 March at 00:30 UTC, strong southerly winds along with a weakly stable stratification created a shear-driven, turbulent PBL. In the profiles of the mean horizontal wind components ( Fig. 14a-b), grid sensitivity is mainly restricted to the first 20 m within the canopy layer. Inside there, the run with default 490 grid-spacing of 40 m results in a slightly higher wind speed compared to the runs with a better resolution of buildings. Profiles of Θ v show negligible sensitivity (Fig. 14c), while for scalar c BC concentrations within the urban canopy are slightly higher in the 20 m and 10 m runs compared to the 40 m and 5 m runs (Fig. 14d). Significantly more sensitivity is observed in the vertical turbulent fluxes of momentum ( Fig. 14e-g) and scalar BC (Fig. 14h). Although the subgrid-scale contributions (dotted lines) become larger as the resolution is decreased, they do not seem to compensate for the loss of resolved fluxes. Apparently, this issue is not restricted to the urban canopy, but may be influenced by an underestimation of vertical wind shear just above the roof tops in the coarser runs. Nevertheless, sensitivity in c BC is very low, arguably because transport is mostly horizontal in the shear-driven case. Thus, the profiles of c BC and c BC w are only weakly related to each other.
For the featured convective case around 2 March 12:30 UTC, horizontal wind speeds are by one order of magnitude lower compared to the first case, but a positive vertical heat flux is responsible for turbulence generation this time. The averaged 500 wind profiles (Fig. 14a-b) show a similar sensitivity to the already discussed shear-driven case, which may indeed indicate an increase of the prescribed roughness length of the subgrid-scale building structure for momentum transfer in future simulations with diffuse buildings. Compared to the shear-driven case, a more substantial sensitivity in Θ v of approximately 0.5 K and also in c BC can be observed for the height range within the urban canopy (Fig. 14k-i). In contrast, a negligible grid-sensitivity is observed for the vertical flux c BC w across the full height range (Fig. 14q). The reason for this behavior in the given 505 convective case is that scalar transport is dominated by vertical mixing. Since the emissions are constant in all simulations, also no variations in c BC w are expected when assuming an equilibrium state. This, however, does not imply that the profiles of eddy diffusivity are constant across the simulations, as the vertical gradients in c BC adjust to match the flux profile. In fact, there must be a significant sensitivity of the vertical eddy diffusivity within the height range where the scalar profiles start to diverge, which is not just by coincidence also the area with the largest instability (the largest super-adiabacity in Θ v ). In the 510 currently used subgrid-scale model, the turbulent Prandtl number P t , which relates the eddy diffusivity to the eddy viscosity, is not further decreased below the neutral value of 0.66 for unstable stratification. However, in unstable conditions P t may be actually much lower, implying a larger eddy diffusivity. An adjustment of the stability-dependent P t may be needed in future simulations.

515
For the evaluation of grid sensitivity at the air-monitoring sites, the horizontal grid spacing of the locally nested domains centered at roadside-classified air-monitoring sites (LE, LL, or LC, respectively) is varied between 40 m, 20 m, 10 m, and 5 m.
The finest resolution of 5 m permits conventional building-resolved simulations, but is omitted for the background stations LT and LW due to the expected low model sensitivity there. Model results are again hourly averaged and spatially interpolated to the exact locations of the measurement sites. In Figure 16, the obtained time series are plotted against each other and against 520 the measurements as reference, similar to Fig. 13. In addition, the corresponding FB values in relation to the measurements are listed in Tab. 3. For the BC background station LT (Fig. 16a), little grid-size sensitivity is found as to expect. At best, BC peak concentrations around 1 March 18:00 UTC tend to be slightly higher with increased resolution, which is mainly from the sensitivity of the subgrid-scale parameterization, as the PBL was stably stratified during this time. FB varies only slightly from 0.05 to 0.00 for the 10 m grid spacing. More interesting are the results for the roadside station LE (Fig. 16b). For this 525 site, the 40 m resolution resulted in an underestimation of BC peak concentrations. In fact, decreasing the grid spacing down to lower or equal 10 m results in a significant better representation of both evening rush-hour peaks. However, absolute peak concentrations still cannot be fully recovered. Interestingly, no further improvements can be achieved with the finest resolution of 5 m. Aside from the discussed rush-hour peaks, model sensitivity is considerably lower, and it seems like the 10 m and 20 m grid spacings result in slightly higher concentrations compared to the 40 m and 5 m resolutions, which are close together. The 530 improved peak representation results in a slightly lower model FB ranging from 0.06 to 0.16.
For the station LC (Fig. 16c), a distinction of the evening rush-hour peaks from the remaining time series turns out to be reasonable too. For the two evening peaks, a decrease of peak BC concentrations with increasing spatial resolution can be observed, which is in contrast to the sensitivity at the station LE. A closer inspection of the spatial concentration gradients in Figure 15 reveals the reason for this contrasting sensitivity. The station LC lies next to two traffic lanes to the north in 535 the model. In the 40 m run, the exhaust plumes spread over a comparatively large area, causing a smearing of the gradient in the vicinity of the air-monitoring site. In the 5 m run, spatial gradients near the road are much sharper, which places the measurement site mostly outside the exhaust plumes. It is, however, questionable if the line-source representation of traffic emissions is still adequate in combination with such a fine grid spacing, as in reality emissions can be effectively spread over a larger area by car-induced turbulence (Gross, 2016). Grid-sensitivity at the PM 10 -measurement site LW is negligible (Fig. 16d), while it is again much more significant at the streetcanyon site LL, where bulk PM 10 is influenced to a large degree by traffic emissions. Not surprisingly, decreasing the grid spacing to 20 m results in higher modeled peak concentrations compared to the 40 m resolution. A further decrease of the grid spacing leads to some indecisive changes for the first evening peak and no significant changes for the second peak. As a result 555 of the higher modeled peak concentrations in the higher resolved runs, an initially moderately positive FB = 0.21 of the 40 m is turned into slightly negative FB values ranging from -0.11 to -0.19. Finally at site LC, a quite similar behavior to the already discussed pollutant BC can be observed, albeit with an extenuated amplitude from the higher background influence of PM 10 .
Hence, FB varies only slightly between -0.46 and -0.32 for the set of sensitivity runs.
Having discussed in detail the grid-sensitivity of modeled BC and PM 10 at the different measurement sites, it remains to 560 quantify this sensitivity in proportion to the absolute simulation error in order to give a final conclusion. In addition, also the Pearson correlation coefficient r can serve as a criterion to judge the agreement of model results with observations. The model error is computed by following equation, where c mod and c obs are the modeled and measured concentrations, respectively, and t is the time indexing. Both and r values computed for all sensitivity runs are additionally listed in Tab  In this study, we applied the dispersion model CAIRDIO for the first time on a real mid-sized city to simulate dispersion of the pollutants PM 10 and BC using a realistic meteorological setup, which was interpolated from a hosting mesoscale simula-580 tion. For the simulation period, two consecutive days in early March 2020 were selected. During this time, unsettled weather conditions with changing PBL characteristics and a generally pronounced magnitude of the intra-urban variability due to relatively low background pollution concentrations prevailed. The horizontal model resolution was set uniformly to 40 m, which permits only to resolve the largest building structures, like industrial sites, while the majority of buildings within the city is described as diffuse obstacles. Nevertheless, the LES approach allows for an explicit representation of the most important 585 turbulent PBL processes, which also include effects from a thermal surface forcing essential to the evolution of the PBL. This capability of the dynamical approach can be considered as a major advantage over more idealized models considering such effects, like stratification, only in parametric form (e.g. Gaussian plume or street-canyon models). In fact, the modeled PBL in this study showed turbulent features, which were consistent with the expected qualitative characteristics based on thermal stratification and vertical wind shear alone. Periods of intermittent or absent turbulence occurred when the critical gradient 590 Richardson number was exceeded. The different dominating dispersion pathways during a specific PBL state (horizonal advection for shear-driven PBL vs. vertical turbulent mixing for convective PBL) resulted in visible qualitative differences in modeled near-surface BC concentrations, like a significantly higher background concentration and smoothed-out gradients in the shear-driven case compared to the convective case. The quantitative evaluation of modeled pollutant concentrations at the air-monitoring sites representative to the urban background showed a diurnal variability in modeled BC concentrations that 595 was consistent with the measurements and thus provided further evidence for a realistic model representation of PBL transport processes. The model agreement with the measurements was also better for BC than PM 10 , as BC is more locally influenced, while PM 10 includes not only predominantly regional influences, but also uncertainties in the complex precursor chemistry.
Ultimately, the model representation of the intra-urban variability of BC and PM 10 concentrations was evaluated using the measurements at the road sites. These were distinct from the measurements at the background sites by the significantly ele- need to be carried out in future studies focusing mainly on such aspects. Also, we assumed the validity of the downscaled 625 surface potential-temperature fields prescribed in the simulation, which affect also the heat flux from building surfaces, mainly in lack of a more physically-based alternative. Here, the further comparison with a microscale model equipped with an own radiation and surface scheme could provide confidence, which was however out of the scope of this study. Finally, in spite of the observed and discussed sensitivity, the comparison of the error in modeled concentrations at the measurement sites showed only slight improvements with an increasing grid resolution, if any at all. For a more significant model evaluation, definitely 630 a more prolonged simulation period needs to be investigated. Still, the findings from this study point to the necessity of more accurately representing other non-physical components in the model, in order to benefit from a more accurate representation of model physics with building-resolving grids. Most notably to mention are the traffic emissions with their currently limited accuracy and comprehensiveness, which may be improved in future simulations with the incorporation of real-time traffic-flux data. Nevertheless, with the currently available data the showcased modeling approach performed at urban gray-zone horizon-635 tal resolutions showed to be a very promising tool for application on more targeted research questions that previously relied on mesoscale model outputs, like, e.g., urban population exposure to air pollution.
Code availability. The source code of CAIRDIO model version 2.0, as well as utilities for data pre-processing are accessible in release under the license GPL v3 and later at https://doi.org/10.5281/zenodo.5585383 (Weger et al., 2021b).
Data availability. The data used in this study, which include model results and air-monitoring observations in part provided by LfULG and In CAIRDIO v.1.0, advection used linear 5 th -order reconstructions with additional limiting for positive scalars. It is well known that such upwind-biased odd-order schemes result in numerical diffusion, as the leading error term is diffusive. In LES, numerical diffusion has a detrimental impact on the correct energy transfer, as excessive energy is drained from the smallest scales that feed on the larger eddies, thus affecting the entire energy cascade. The manifestation of excessive damping can be seen in the energy spectra of Figure A1. Nevertheless, some sort of numerical damping of the smallest wavelengths is necessary 650 in order to maintain numerical stability, as these scales are flawed by large dispersion errors. Recognizing that the standard 5 th -order linear upwind formulation carried out in each dimension separately is simply the addition of a high-order diffusion term to a non-diffusive 6 th -order central scheme, directly leads to an opportunity to control numerical diffusion: Here, ∆ h is the grid spacing, u + the positive definite transport velocity component, and ∇ 5th a finite difference operator 655 using 5 th -order reconstructions. The product u + ∆ h is called numerical diffusion coefficient, and mainly acts in the dominant transport direction. In the revised scheme, u + is replaced by a constant parameter d = 0.05. Note that this results in isotropic diffusion, similar to the method proposed in Xue (2000) for high-order damping.
A0.2 Prognostic TKE formulation for subgrid-scale mixing 660 In version 1.0 an algebraic eddy-viscosity formulation was used. Therein, eddy viscosity was diagnosed from the strain-rate tensor S without taking buoyancy effects into account. In order to simulate non-neutral PBLs, we implemented a prognostic subgrid-scale TKE formulation similar to Deardorff 1973 in version 2.0. This scheme not only takes buoyancy effects into account but also avoids the local-equilibrium assumption and thus may provide more accurate results with coarse grid spacings.
The first two terms correspond to the advective-diffusive transport, which also incorporate the pressure-correlation term parameterized by a doubling of the diffusive flux. The shear-production term is parameterized with the squared magnitude of l sgs = min ∆, 0.76 √ e/N c , where ∆ = (∆ x ∆ y ∆ z ) 1/3 is the static (grid) mixing length. The eddy viscosity coefficient k m is parameterized according to: Note that we replaced l sgs with ∆ therein, as the original formulation resulted in a too small eddy diffusivity for stable stratifications within the roughness sublayer of diffuse urban canopies. We argue that in such a case, the mixing length is Lastly, the eddy diffusivity k h for heat and scalar transport is related to k m by the inverse turbulent Prandtl number: In Deardorff (1973), P r −1 t is parameterized according to P r −1 t = 1 + 2 l sgs ∆ .
However, as pointed out in paragraph 3.3.1, the neutral value of P r −1 t may be too low in unstable conditions (cf. Li et al. 2015), and if a negative impact from this is further corroborated, the stability dependency will be entirely replaced by another formulation in a future model version.

A0.3 Dry and wet deposition of particulate matter
For a consistent model description of particulate matter dispersion, dry and wet deposition processes are considered too in the new model version 2.0. For dry deposition, the scheme of Zhang and He (2014) was implemented, which considers the bulk-size categories PM 2.5 , PM 2.5−10 and PM 10+ . Accordingly, the deposition flux of a given particle category on a horizontal surface is given by the particle mass concentration c pm times the parameterized deposition velocity v d , which contains the 690 contributions from gravitational settling v g , aerodynamic resistance r a and surface resistance r s :