Including the Urban Canopy Layer in a Lagrangian Particle Dispersion Model

In this study we introduce a novel extension of an existing Lagrangian particle dispersion model for application over urban areas by explicitly taking into account the urban canopy layer. As commonly done, the original model uses the zero-plane displacement as a lower boundary condition, while the extension reaches to the ground. To achieve this, spatially-averaged parametrizations of flow and turbulence characteristics are created by fitting functions to observational and numerical data. The extended model is verified with respect to basic model assumptions (well-mixed condition) and its behaviour is investigated for unstable/neutral/stable atmospheric stabilities. A sensitivity study shows that the newly introduced model parameters characterizing the canopy turbulence impact the model output less than previously existing model parameters. Comparing concentration predictions to the Basel Urban Boundary Layer Experiment—where concentrations were measured near roof level—shows that the modified model performs slightly better than the original model. More importantly, the extended model can also be used to explicitly treat surface sources (traffic) and assess concentrations within the urban canopy and near the surface (pedestrian level). The small improvement with respect to roof level concentrations suggests that the parametrized canopy profiles for flow and turbulence characteristics realistically represent the dispersion environment on average.

air exchange, often lead to high concentrations of air pollutants that adversely affect human health (e.g., US EPA 2011).
Consequently, many operational and experimental air pollution dispersion models are used in urban areas to forecast air pollution concentrations. Examples include models like AER-MOD (Cimorelli et al. 2005), QUIC-PLUME (Williams et al. 2002), MicroSpray (Tinarelli et al. 2012, and references therein), GRAMM/GRAL (Oettl 2014(Oettl , 2015, and SIRANE (Soulhac et al. 2011). All these models have in common that they do-based on various approaches-include the mean flow through the region between buildings, i.e., the urban canopy layer (UCL, see Fig. 1). How they account for the between-building turbulence differs, however.
It is possible to explicitly resolve all obstacles in a type of large-eddy simulation (LES) (e.g., Auvinen et al. 2020) and derive the turbulence characteristics from there. This, however, is computationally highly expensive and not a suitable approach for operational applications. Some dispersion models (e.g., QUIC, GRAMM/GRAL, MicroSpray) use types of Reynoldsaveraged Navier-Stokes approaches approaches to calculate a mean flow field and derive the second-order moments from turbulence kinetic energy, combined with standard boundarylayer parametrizations. We believe that parametrizations specifically designed for the urban boundary layer could be useful for these types of models. This approach is computationally less expensive than LES, but still associated with considerable costs. Note that the main cost factor of all these approaches is the building-resolving flow simulation, not the dispersion. Furthermore, if there are no building-resolving maps available for a specific city, or the computational effort to explicitly calculate flow around resolved buildings is too high, another solution is required to include urban effects.
In this paper, we describe a method for incorporating the UCL in horizontally homogeneous flow models. Specifically, we propose describing the UCL as a partially porous medium with spatially-averaged vertical parametrizations of mean and turbulent flow. To our best knowledge there are no previously published parametrizations of second-order moments, dissipation rate, and vertical skewness specific to the UCL, even though the method of using specialized parametrizations characteristic for the canopy layer is not new (e.g., for vegetation canopies Baldocchi 1997).
Over a rough, urban surface, the surface layer consists of the roughness sublayer (RSL) and the inertial sublayer (ISL), as shown on the right-hand side of Fig. 1. The RSL extends from the surface up to the blending height z * -the height of the maximum Reynolds stressand includes the urban canopy layer (UCL). Reynolds stress has been found experimentally and through simulation to be strongly height dependent in the RSL, i.e., to decrease to very small values close to the zero-plane displacement d. Figure 1 (left) shows the conceptual sketch of a Reynolds stress profile in the urban boundary layer approaching in its lowest part a 'constant stress' portion-as expected in the surface layer. The dashed line depicts the corresponding decrease in magnitude of Reynolds stress in the RSL, i.e., if large roughness elements are present. This dashed line is based on the Reynolds stress profile of the urban RSL introduced by Rotach (2001). He showed that including this Reynolds stress profile into a Lagrangian particle dispersion model (LPDM) has substantial impact on modelled downwind concentrations and that the model performance improved significantly. However, due to the chosen parametrization, the lower boundary condition had to be set to the height of the zero-plane displacement. Thus the lowest tens of metres are still not included in the model domain. In this paper, we aim at explicitly introducing the UCL in the LPDM by extending the Reynolds stress profile down to the ground, as sketched by the blue line in Fig. 1. We also introduce similar extensions to the other necessary profiles of flow and turbulence. We will call the original model-including its urban RSL parametrization-RSM (roughness The solid line is the parametrization for the outer layer and the inertial sublayer from de Haan and Rotach (1998), from the boundary-layer height z i down to the zero-plane displacement d. The dashed line is the RSL parametrization of Rotach (2001) from the blending height z * down to d. The blue line is the parametrization proposed in Sect. 4 from the mean building height z h down to the ground sublayer model), while the present model with the extension down into the urban canopy layer will be called ULM (urban canopy layer model). When referring to general properties of a Lagrangian particle dispersion model, we will refer to an LPDM. Wang et al. (2018) suggested a similar approach of including the UCL in a backward LPDM to calculate footprint functions, but their parametrizations are base on the assumption of Monin-Obukhov similarity theory inside street canyons, which seems to be hard to defend (e.g., Rotach 1999).
In this work, we test the sensitivity of the concentration output to the prescribed turbulence and wind profiles and validate the new model using field experiment data. Ultimately, we seek to determine whether or not a non-building resolving approach of including building effects is viable for simulating dispersion in urban areas. Then, the model may also be used as the core for an urban footprint model, for which it is necessary, or at least highly desirable, to have a domain extending to the physical surface where many potential sources (e.g., traffic) are present.
In the following, Sect. 2 summarizes the original formulation of the LPDM, Sect. 3 presents the underlying datasets used for the UCL parametrizations, and these are then described in Sect. 4. The ULM is verified and tested on its behaviour compared to the RSM in Sect. 5. A validation against the the Basel Urban Boundary Layer Experiment (BUBBLE) dataset is shown in Sect. 6.

Lagrangian Particle Dispersion Model
This study uses an LPDM initially developed by Rotach et al. (1996), later extended by de Haan and Rotach (1998) (crosswind dispersion), de Haan (1999) (more efficient concentration calculation), and Rotach (2001) (urban RSL parametrization). Rotach et al. (2004) compared this urban model with a specifically designed urban tracer experiment and found the model performance to highly depend on "the exact form of the parametrization for the flow and turbulence structure within the urban roughness sublayer". Gibson and Sailor (2012) suggested several corrections to mathematical formulations in Rotach et al. (1996), but Stöckl et al. (2018) found that most of them had been either implemented in the model code long before or mistakes existed only in the publication, not the code. Finally, a more comprehensive summary of the model and a first attempt to include the effects of street canyon flow in the model can be found in Stöckl (2015).
The present LPDM is a one-dimensional (vertical) model with dispersion in three dimensions. The model is horizontally homogeneous, stationary, and does not include chemistry or deposition. The model's domain ranges from the atmospheric boundary layer height z i to the zero-plane displacement d over urban areas. As mentioned above, in this study we extend the domain further towards the ground surface. Rotach et al. (1996) were the first to enable their model to simulate in both convective and stable to neutral conditions. For this, they used the approach by Thomson (1987) and defined the necessary probability density function (p.d.f.) of the particle velocities as a mixture of two limiting states: (1) The first term describes the joint p.d.f. of a convective atmosphere with no correlation between any of the wind components, Gaussian p.d.f.s P u and P v , and P w,c a skewed p.d.f. for the vertical component w . The second term in Eq. 1 describes a neutral or stable atmosphere with purely Gaussian p.d.f.s for u , v , and w , but also a correlation between u and w . Using this combined p.d.f. with the transition function f enables the model to be valid in different atmospheric conditions, while still fulfilling the well-mixed criterion (Thomson 1987). The transition function f is zero everywhere in stable and neutral conditions, while in convective conditions it ranges from one in the mixed layer to zero near the ground (see Rotach et al. 1996 for details). The model uses an explicit Euler forward scheme to calculate the next position of each particle from the current position t: where i = 1, 2, 3 for the three Cartesian directions, x i is the position and u i the velocity of the particle, and du i is a velocity increment. This velocity increment or acceleration is derived with the Langevin equation: where dt is the timestep, a i the correlated and b i j the random part of the acceleration. Here, ξ j describes a Wiener process with mean zero and variance dt (Rotach et al. 1996). Since this equation cannot be solved analytically, Rotach et al. (1996) followed Thomson (1987) in using the stationary Fokker-Planck equation, to describe the same process as Eq. 4. Here P tot is the total p.d.f. of the particle velocities (Eq. 1), as well as the p.d.f. of the Eulerian fluid velocities under the well-mixed assumption (Thomson 1987). See Rotach et al. (1996) for details on the analytical solution: and The functions i are given in Rotach et al. (1996) and C 0 = 3 is the universal inertial subrange constant, see Rotach et al. (1996) for a discussion. The model is completed by assuming vertical profiles of flow and turbulence characteristics -including the dissipation rate of turbulence kinetic energy ε-and thereby defining all p.d.f.s used in Eq. 1, although it would be possible to use output from a numerical model instead (see, e.g., Weil et al. 2004).

Urban Canopy Layer Data
To achieve the goals outlined in Sect. 1, values of ε, u, u w , u 2 , v 2 , w 2 , and w 3 are required for every possible (vertical) position in the model domain. Our model is horizontally homogeneous and thus does not rely on information on topography or building geometry. Instead we use spatially-averaged, vertical profiles over the widest feasible spread of different urban areas as possible, similar to the work of Raupach et al. (1996) for vegetation canopies. Up to date, vertical profiles for a representative spatial average from full-scale experiments are still scarce. Hence most of the following datasets originate from wind tunnel or numerical studies. Note that we refrain from using the classical <> spatial averaging notation for brevity. Unless noted otherwise, all profiles are spatially averaged.
Concerning spatial averages in the UCL, it is critical to distinguish between two types, which we will call 'intrinsic' and 'superficial' after Schmid et al. (2019). An intrinsic spatial average ignores the volume filled by assumed solid obstacles and averages only over the fluid volume. In contrast, a superficial spatial average takes the volume filled by obstacles into account and therefore also depends on the porosity of the city. In principle, it would be preferable to use superficial averages (as explained by Xie and Fuka 2018), because then the canopy could be treated as a homogeneous, porous medium (Böhm et al. 2013). However, most of the datasets only provide the intrinsic averages and not enough data to use porosity to convert. Instead, we average intrinsically and account for the 'missing' building volumes by adjusting the lower boundary condition at 'the surface'. Usually, particles are 'reflected' at the lower domain boundary z r by inverting the sign of w and u , as well as adjusting the new z-position in such a way that the total distance traveled in this timestep remains the same, but the particle ends up above z r instead of below (Thomson and Montgomery 1994). This sustains the particles in the model domain and is done similarly at the top domain boundary z i . In our approach we introduce an elevated partial reflection level at the mean building height z h , where particles are reflected with a chance of λ p when arriving from above. Here, λ p is the plan area fraction of roughness elements (Grimmond and Oke 1999) and thus mimics roof reflection.
The datasets used in this study are briefly introduced in the following. Coceal et al. (2006) (from now on CTCB06ST: staggered, CTCB06AL: aligned, CTCB06SQ: square array) and  (CDTB07) used direct numerical simulations to predict the flow over uniform cubes in regular grids of different layouts (staggered, aligned, square). We used data from their Figs. 20 and 2, respectively. Auvinen et al. (2020) simulated synoptically forced wind from the sea into the coastal city Helsinki, Finland, with an LES model and shared their spatially-averaged data with us (HEL). Xie and Castro (2009) simulated flow through an intersection in London, UK, (DAPPLE project site) using an LES model and reported the spatially-averaged profiles in their Fig. 9 (XC09). Carpentieri et al. (2009) simulated the same intersection as Xie and Castro (2009) in a wind tunnel (1:200 scale) and shared 14 profiles from different locations indicated in their Fig. 3 with us (CRB09). Harman et al. (2016) studied the flow between a model canopy of regularly spaced, thin 'tombstone' obstacles in a wind tunnel. We used data from their Figs. 2, 3, and 4 (HBFH16). Böhm et al. (2013) showed profiles within and above a canopy of solid tree-shaped obstacles (light bulbs) in their wind tunnel that "share[s] characteristics of both vegetation and urban canopies" in their Fig. 10 (BFRH13). Kastner-Klein and  investigated a model of Nantes, France, in a wind tunnel and measured vertical profiles. We used up to 42 of them to calculate the spatial average, depending on the specific variable (KKR04).
In field studies it is prohibitively difficult and expensive to measure enough vertical profiles for spatial averaging, hence we are unaware of any such measurement campaigns. However, Rotach (1993) argued that averaging over multiple profiles of few towers can be used as a raw estimate for a spatial average, considering the wind direction and therefore the upstream geometry changes.  followed this approach for the BUBBLE project in Basel, Switzerland , and provided his dataset from two towers (BUBBLEU1 and BUBBLEU2). Since spring 2017, highly resolved fluxes of momentum, shear, water vapour and chemical compounds have been measured by the Innsbruck Atmospheric Observatory (IAO) at the University of Innsbruck, Austria (Karl et al. 2020;Ward et al. 2022). Data from two sonic anemometers outside the RSL (42.2 m a.g.l. and 39.6 m a.g.l.) and one streetlevel sonic anemometer (3.0 m a.g.l.) were used herein (ACINN). Although the two upper measurements were outside the RSL, we still chose to use the data of the street level to have more data sources for ε in the UCL. These real-world datasets are, on the one hand, suitable for our purpose because they reflect reality. On the other hand, they do almost certainly not represent true spatial averages, so care has to be taken not to over-interpret them. Giometto et al. (2016) investigated the suitability of a single tower measurement as source of spatial average and concluded that the single tower may be severely biased and is therefore unsuitable. However, they arrived at this conclusion by looking at LESs in comparison to a tower measurement, but only looked at two simulations with differing wind directions. Considering more wind directions may somewhat alleviate these results.
An overview of the data sources can be found in Table 1. Since these data are from different sources, their normalizations, scalings, and rotations differ as well. We rotated the profiles into the prevailing wind direction with a single rotation, where appropriate. If the scaling and normalization corresponded to those used herein (see below), they were left as is. Otherwise the profiles were scaled and normalized, as explained in the following paragraphs.
Traditionally, the atmospheric boundary layer is scaled with z i , at least for the outer layer (e.g., Stull 1988), and so are our model's parametrizations too, as described in Rotach et al. (1996). In the urban RSL the height is scaled with z * instead, which corresponds to defining z * as the height of maximum (magnitude) Reynolds stress, such that the peaks of Reynolds stress collapse. Canopy scaling for vegetation canopies is usually done via the mean canopy height z h (e.g., Raupach et al. 1996;Rannik et al. 2003), but there the top of the canopy often coincides with the peak of Reynolds stress and therefore incorporates the entire RSL. This means that in classical canopy scaling (e.g., Raupach et al. 1996) the peaks of Reynolds stress and velocity variances are situated at z/z h = 1. Urban geometries with uniform height show All sources deliver profiles of u, u w , and u 2 , 'additional profiles' in the last column refer to additional variables provided. WT means wind tunnel, FS full scale experiment, LES large-eddy simulation, DNS direct numerical simulation similar characteristics, but several studies have noted that this is not true for city geometries with non-uniform height (Xie et al. 2008;Xie and Castro 2009;Carpentieri and Robins 2015), where the peaks are positioned further aloft. This the origin of the z * definition. Multiple studies suggest the use of the tallest upstream building as scaling height (Xie et al. 2008;Xie and Castro 2009;Kanda et al. 2013;Inagaki et al. 2017;Sützl et al. 2020) instead. Alternatively, it has been suggested (Martilli et al. 2000, as mentioned in Rotach 2001, to use z h + σ h as scaling height, where σ h is the standard deviation of the building heights. Here we find that z h + 1.5σ h fit our datasets best (not shown). Hence, we use z * = z h + 1.5σ h when RSL scaling is required.
To extend the RSL scaling of Rotach (2001) down to the surface we use canopy scaling (z/z h ) for heights z < z h . The transition from the RSL parametrization to the UCL parametrization thus occurs at z = z h . If z * is expressed in terms of z h , we can present the entire RSL profiles as a function of z/z h .
The variables of interest are normalized with powers of u * ,I SL (e.g., Raupach et al. 1996), where the roughness velocity in the inertial sublayer (ISL) u * ,I SL is obtained via the method of Kastner-Klein and . Two notable exceptions are the mean wind speed, which is scaled by its value at z h (e.g., Raupach et al. 1996), and the skewness of the vertical velocity component, which is not scaled. Figure 2 shows a summary of all datasets mentioned previously, specifically profiles of interest in their intrinsically spatially-averaged scaled and normalized forms. As mentioned, the canopy scaling is only intended to be useful for heights below z/z h = 1, hence the profiles are not expected to collapse for z > z h .
The numbers in the legend are σ h /z h , indicating the uniformity of the height distribution. To avoid visual overload, only a selection of data points are shown with markers In Fig. 2a, normalization and scaling causes all profiles of the mean wind speed to collapse onto (1,1). Below z h most profiles exhibit similar behaviour (reminiscent of an exponential profile) and group together tightly, with the exception of three measurement points from real world datasets. The low scatter implies that the scaling and normalization is successful. However, our dataset is too sparse to evaluate whether or not a um scales with λ p , as sometimes suggested (e.g., Castro 2017 for a discussion). Figure 2b shows the covariances u w , which exhibit their minima at different heights. However, all profiles have their minima at or above z h , with the exception of CRB09, which is the wind-tunnel simulation of a real urban area with a focus of one intersection, where all measured profiles are located within or around this intersection. Therefore the profiles are not ideally placed for a spatial average, because that was never the intention of Carpentieri et al. (2009). We have decided to include these profiles nevertheless, because datasets with non-idealized geometries are rare. The configuration of buildings directly upstream of the intersection may have lead to a non-representative z h for our spatial average, given the fact that most of these buildings are uncharacteristically low, leading to an extremum of u w even below the nominal z h . Furthermore, a thin tower on top of one building directly upstream results in a second local minimum at z/z h = 1.7 to 2.0. Interestingly, all profiles based on uniform-height geometry, which means σ h = 0 (see legend of Fig. 2), indeed exhibit their minimum in Reynolds stress at or slightly above z h , while cities with non-uniform heights have their minimum much higher up (except CRB09).
There are far fewer data available for the dissipation rate of turbulence kinetic energy ε (Fig. 2c), but the dissipation rate seems to decrease in the UCL. More research would be needed.
Contrary to the other profiles, the skewness of the vertical velocity component Sk w (Fig. 2d) is scaled with RSL scaling instead of canopy scaling. For a detailed justification see Sect. 4.4. Near z h , the values are clustered close to zero, reflecting the Gaussianity of the mechanically induced shear turbulence. Below, in the UCL, it is noticeable that all profiles exhibit negative values of Sk w , similar to vegetation canopies (Raupach et al. 1996). Near the ground, the values appear to return to zero. Figure 2e shows profiles of the longitudinal velocity variance u 2 . Similar to u w , they peak-this time a maximum-at or above z h . In the UCL the curvature of the profiles seems to range from positive to negative, most likely depending on the specific city geometry. For example, the contrast between CTCB06ST and CTCB06AL is striking, despite the fact that the only difference between them is the different cube layout in Coceal et al. (2006). However, most real cities would not conform to either end of these extremes and in fact the profiles CRB09, HEL, and KKR04, all more or less reflecting real cities, are somewhat in the middle of the spread. Conversely, XC09 is also a simulation with real city geometry but on the extreme left of the spread. Nevertheless, the general shape of peaking at or above z h and then declining towards the ground can be observed in all profiles.
The lateral velocity variance v 2 in Fig. 2f exhibits similar characteristics as u 2 , except that the curvature of the declining profiles is more or less always the same. The same is true for w 2 in Fig. 2g. Here the profiles seem to collapse best, except that both KKR04 and CRB09 show substantially larger values for w 2 within the UCL.

Urban Canopy Layer Parametrizations
The following sections introduce the turbulence profile parametrizations necessary to run the LPDM. They are all devised by imposing a shape that is inspired by the respective profiles in Fig. 2 and additionally through necessary boundary conditions. With a given function and boundary conditions, each profile is then numerically fitted to the literature data described in Sect. 3.
One common boundary condition is continuity with the existing profiles aloft at height z h , with the exception of w 3 , which is continuous at z * . The 'transition value' (denoted by subscript 'h') of each profile is the value of the individual RSL parametrization at height z h .
Since the profiles (i.e., their parametrizations) in the RSL depend on the (variable) meteorological conditions and on urban geometry (building height), continuity of those profiles towards the UCL parametrization would require a fit of the canopy profiles for each simulation separately. To avoid doing this during each model run, we scaled height with z h and normalized all values of the experimental data with their value at z = z h , thus forcing all the data (and the corresponding parametrizations) through (1, 1). This point (1, 1) is then the point where the profiles transition from RSL parametrization above to UCL parametrization below. To keep the numerical fitting algorithm independent of the amount of data points per profile, the data points are weighted in such a way that each individual profile has equal impact. Since we aim to test the model's sensitivity to the fitted parameters (see below, Sect. 5.3), we need a range of possible values for each parameter. We individually fitted the range of the parameters to encompass the approximate range of the available dataset. This range approximates both the uncertainty introduced by having only a few datasets, measurement and simulation uncertainty, as well as uncertainties introduced by the fitting procedure.
The profiles in Fig. 2 also diverge above z h , which cannot be taken into account using the set-up in this work. Since we normalize all profiles to unity at z h and use another parametrization aloft, it does not matter for our purpose that they diverge further aloft. Furthermore, not all profiles of second-order moments peak at height z h . Generally speaking, canopies with uniform heights have a lower maximum of the Reynolds stress and velocity variances compared to more realistic city geometries. However, at height z h all types of profiles are declining from the peak, which means that the general shape of the profile below z h is conserved (see Fig. 3).
In the following, the details of the parametrizations are outlined separately for each variable. The resulting parametrizations are-together with the datasets from Fig. 2-displayed in Fig. 3.

Mean Wind Speed
Mean wind speed in the UCL is parametrized with the exponential function: following common practice in vegetation canopies (Cionco 1965). Castro (2017) has questioned the adequacy of exponential urban canopy profiles for mean wind-his argumentation largely being based on the assumptions needed for its derivation not being fulfilled, and on the failure of some experimental profiles being exponential. Among the alternatives he discusses, there is the approach by Yang et al. (2016), who propose combining a logarithmic profile above the canopy with an exponential profile, i.e., exactly corresponding to the general approach followed herein. Yang et al. (2016) show that for different values of λ p = λ f (for cubes at a non-oblique approach-flow angle) and σ h , an exponential function is a good fit to their simulated wind profiles in the upper 70-80% of the UCL-and even suggest that the validity of the exponential profile might be extendable closer down to the surface (Yang et al. 2016). λ f is the frontal area density after Grimmond and Oke (1999). Wang et al. (2018) use the same Eq. 8 and a method proposed by Macdonald (2000) to determine a um based on building geometry. However, this approach depends on knowledge of the height profile of the sectional drag coefficient and a mixing length scale, both of which are difficult to determine for non-idealized real urban geometries-and therefore unknown for most of the datasets used here.
Consequently, we determine a um using the boundary condition for the canopy profile, i.e., the requirement for a smooth transition at z h . For this we normalize Eq. 8 with u h and scale it with z h , resulting in:û for z = z z h . Then, we fit Eq. 9 to the UCL data gathered from literature (see Fig. 2) and thus determine a um = 1.97 with an uncertainty range of [0.64, 3.31] (see orange shaded area in Fig. 3a). Based on the range of a um in the literature, Macdonald (2000) suggests a um = 9.6λ f , which would correspond to about 3 to 4 for the datasets used in this study. Ramirez et al. (2018) report values < 2. Note that this profile does not generally satisfy the no-slip condition, but a shallow reflection layer at the bottom (see Sect. 4.7) circumvents this.

Turbulent Fluxes
The LPDM assumes no directional shear (i.e., v w = 0) in the surface layer and uses u * = −u w . Rotach (2001) derived the profile of the local u * ,l = 4 u w 2 + v w 2 and then used this in the RSL as basis for all other turbulence profiles. We define the velocity variances, skewness, etc. separately in this study and therefore the profile of u * ,l has far less influence in the UCL. Furthermore, not all datasets provided v w and thus we chose consistency over completeness and did not include v w in u * UCL . It is noted that neglecting v w is not supported by many of the datasets introduced in Sect. 3 and remains an inherent limitation of the model itself that will be addressed in a future study.
The shape of the data in Fig. 2b suggests a function of the form u * UCL = a 1 z z h 1/a Re +a 3 to be useful, which is subject to the boundary conditions which-normalized, scaled, and converted to u w -yields: Once fitted to data in Fig. 3b, a Re = 0.75, with uncertainty range [0.29, 2.82]. Due to spatial and temporal averaging, the covariances, such as u w contain an additional contribution, i.e., the dispersive stress (Raupach and Shaw 1982). Note that we will use <> to denote spatial averaging in the next three paragraphs, hence u w →< u w >. The wind components can be decomposed into, e.g., u =< u > + u + u , where u is the instantaneous flow quantity, < u > its spatial and temporal mean, u = u − < u > the spatial fluctuation around the time-space mean, and u the turbulent fluctuation in both time and space. Given this, a product of two velocity components averaged over time and space becomes, e.g., < uw >=< u >< w > + < u w > + < u w >, where < u w > is the so-called 'dispersive stress' contribution, which essentially describes spatial deviations from time averaged flow. Similar to Reynolds stress, dispersive stress represents momentum transport, for example by a canyon vortex.
Without immediate effect of obstacles, dispersive stress is zero, but in urban areas it is nonnegligible (Coceal et al. 2006Martilli and Santiago 2007;Xie et al. 2008;Giometto et al. 2016;Simón-Moral et al. 2017;Xie and Fuka 2018). Generally, the value of dispersive stress appears to depend strongly on the city geometry (Coceal et al. 2006. Xie et al. (2008) demonstrate the impact of the standard deviation of canopy height on dispersive stress: for uniform-height canopies it vanishes right above z h , while it only gradually approaches zero for their non-uniform-height canopy. Figure 4 shows this by displaying the ratio of dispersive stress to Reynolds stress taken from various sources approaching zero rapidly above z h for the profiles shown in markers, which stem from canopies of uniform height.
Since the Reynolds stress is always negative for these specific profiles, one can also deduce the sign of the dispersive stress. Dispersive stress is negative for all studies using staggered cubes (stars) and realistic city geometries (coloured lines), but for aligned cubes (dots) it is predominantly positive. This is caused by coherent vortices (Martilli and Santiago 2007) forming in the idealized cube geometries. However, Simón-Moral et al. (2017) show that convective conditions can lead to negative dispersive stresses even for aligned cubes (see Fig. 4), even though their neutral run (not shown) produced mostly positive dispersive stresses. Since we are primarily interested in realistic geometries, we ignore positive dispersive stress.
To investigate the impact of dispersive stress, we have designed a crude description of the ratio between the dispersive stress and the Reynolds stress (black line in Fig. 4). Then, the Reynolds stress parametrization is multiplied by 1+ <u w > <u w > , essentially adding the dispersive stress to the Reynolds stress. Since the dispersive stress appears to strongly depend on the canopy, it is even more unlikely than for other variables that a more general parametrization is possible. The results show hardly any difference between switching dispersive stress on or off (not shown). Firstly, this is due to the rapidly decreasing magnitude of < u w > within the UCL (Fig. 3b) and the therefore small impact of the factor 1 + <u w > <u w > . Secondly, Reynolds stress as a whole has small impact, as shown in Sect. 5.3. Given the large canopyto-canopy variability of the dispersive stress and the relatively minor impact its inclusion has on dispersion within the canopy we decided not to include the dispersive stress contribution into the parametrized profile of Reynolds stress.
For the parametrization of sensible heat flux w θ we use the formulation of  in his Equations 4.41 and 4.42. Note that this parametrization is only valid for convective cases and thus only used in such. This parametrization only affects the convective velocity scale w * , by giving it a local value in the RSL. This w * is used by the LPDM to determine turbulence profiles of RSL, inertial sublayer, and mixed layer in convective situations. Note that the UCL parametrizations do not depend on stability. The dataset in Sect. 3 is too limited to separate it according to atmospheric stability and fit profiles depending on stability. Thus, the profile of w θ does not influence the UCL parametrizations directly, but it does influence the RSL parametrizations and therefore the transition values to the UCL below.

Dissipation Rate of Turbulence Kinetic Energy
Since the data available are scarce (see Sect. 3), the shape of the profile is more uncertain than for other profiles. Nevertheless, the available data suggests a decreasing ε from the RSL downwards, so we choose the same general function as u UCL : where ε h = ε UCL (z h ) = ε RSL (z h ) ensures continuity to the RSL and a ε can be numerically fitted to data after canopy-scaling and normalizing with ε h . The result can be seen in Fig. 3c and leads to a ε = 1.01 with uncertainty range [−0.18, 2.90]. Note that Giometto et al. (2016) shows ε in their Fig. 10 too, but not spatially averaged. Nevertheless, the behaviour in their figure is similar to Fig. 3c, with the exception of a second local minimum close to the ground that is most likely hidden by our reflection layer. Di Bernardino et al. (2020) find that the dissipation rate is not strongly dependent on λ p , above the canopy (see their Fig. 8).

Vertical Skewness
As mentioned earlier, the skewness of the vertical velocity Sk w is negative in the whole RSL.
Since the Rotach (2001) approach of including RSL effects only modifies u * , Sk w remained unchanged from Rotach et al. (1996), because their Sk w parametrization does not depend on u * . Consequently, we choose to have the new parametrization of Sk w affect the whole RSL. Furthermore, it is more common to find profiles of Sk w than w 3 in the literature. Therefore, we show and parametrize Sk w here as well, although the model internally uses w 3 . This is simply addressed by de-normalizing w 3 = Sk w w 2 3 2 in the model. Originally, the p.d.f. P tot was designed to incorporate the effect of convection in the mixing layer, i.e., yielding Sk w > 0. Appendix 1 demonstrates that the same formulation can also be used in the UCL to yield a negative skewness.
Unlike the other parametrizations, it is not reasonably possible to collapse all profiles onto (1, 1), because Sk w ≈ 0 at height z * (see Fig. 2d), which means that both positive and negative values are possible. We choose a parabolic shape Sk wUCL = a 1 (z + a 2 ) 2 + a 3 with boundary conditions (i) Sk wUCL (z * ) = Sk wRSL (z * ) = Sk wt , (ii) Sk wUCL (0) = 0, and (iii) Sk wUCL ( z * 2 ) = −a Sk . The last condition is not necessary to fit the function, but without it the minimum of the fitted function depends not only on the free tuning parameter a Sk , but also on the value of z * , which results in Sk w profiles strongly varying with z * . Condition (iii) also helps to keep the resulting function from being strongly dependent on the value of Sk wt , which is unknown at fit-time and changes from city to city. Since it is not possible to collapse the profiles and the function onto (1, 1) as for the other variables, this was a major concern, because that would mean that the function to be fitted is no longer independent of the specific model run. The third condition alleviates this problem. Taking a Sk as free parameter leads to: which can be scaled via z = z z * . We choose Sk wt = 0.05 as somewhat arbitrary transition value, but-as mentioned before-the function is designed not to be sensitive to this choice. Figure 3d shows the resulting function for a Sk = 0.48 with uncertainty range [0.11, 0.70].

Velocity Variances
Since the profiles in Fig. 2e are of similar shape to the mean wind speed profiles in Fig. 2a, we choose to use the same general function, for the longitudinal velocity variance. With canopy scaling this yieldŝ which is then fitted to similarly scaled and normalized data. The result is a u2 = 1.30 with uncertainty range [0.04, 4.43] and can be seen in Fig. 3e. Using the same arguments as for u 2 UCL , which-when scaled an normalized as in Eq. 15-leads to a v2 = 0.72 with uncertainty range [−0.26, 3.17] [Fig. 3f].
For the vertical velocity variance we note that for all the profiles in Fig. 2g reaching the ground, w 2(z = 0) = 0. Thus we choose w 2 UCL = (a 1 z) 1/a w2 + a 2 with a w2 the free parameter and the boundary conditions (i) w 2 UCL (z h ) = w 2 RSL (z h ) = w 2 h and (ii) w 2 UCL (0) = 0. This leads to: and can be scaled and normalized to:ŵ Fitting Eq. 18 to data results in a w2 = 2.06 with uncertainty range [1.10, 5.89] [Fig. 3g]. Rotach et al. (1996) designed the function f for the transition of the p.d.f. P tot from convective to neutral and stable conditions. Since we have already repurposed the convective form of the p.d.f. for w (P w,c ) to achieve the necessary negative skewness in the UCL, we also have to redesign f to activate this skewed part of the p.d.f. If f would stay near zero close to the ground as proposed in Rotach et al. (1996) for the RSL, the skewness of w would have no effect and the resulting Eulerian p.d.f. of model particles would not be skewed. Consequently, we choose a function that reaches unity quickly, descending from z * . Note that this choice is arbitrary, because the only way to determine it objectively would be to measure the height profile of the w p.d.f. horizontally averaged in the UCL and then find a function f such that the resulting Lagrangian p.d.f. of the model particles agrees with the measured Eulerian p.d.f. of the real-world flow for a given Sk w profile. To our best knowledge no such data can be found in the literature. The boundary conditions are continuity at z * (i) f UCL (z * ) = f RSL (z * ) = f t and (ii) f (0) = 1 for a function of the form f UCL = a 1 + a 2 e a f z , with a f as free parameter, which leads to:

Transition Function
We choose a f = 0.4 and [0.01, 2.0] for the uncertainty range. The resulting functions can be seen in Fig. 3h, although without measurements. Note f UCL being unity at the bottom does not mean that the Eulerian p.d.f. of w is necessarily skewed; this also depends on the value of w 3 , which approaches zero close to the ground in any case. Additionally, f UCL approaching one towards the ground is consistent with u w approaching zero near the ground.

Aspects Not Accounted For
Several of the parametrizations presented are problematic directly at the ground: u UCL does not satisfy the no-slip condition; w 2 UCL is zero at z = 0, which leads to a divide-by-zero error in Eq. 23 and P w (not shown); u * UCL is also zero at z = 0, also leading to a divideby-zero error (not shown). To circumvent these issues, the model does not allow particles to reach all the way to the ground, but reflects them slightly higher up at z 0 . This mimics a shallow surface layer near the ground, where the wind speed decays rapidly to zero (Yang et al. 2016) and where most likely an entire set of new parametrizations would be needed.
While v w is generally small over smooth surfaces or far above rough surfaces-in a wind-following coordinate system-it is not so in the UCL. Unfortunately, our model cannot represent this.
Not taken into account at all in this study is the effect of wind direction and wind direction changes, which are a major factor in plume dispersion (e.g., Michioka et al. 2019).
Ideally, the parametrizations would be able to take different types of urban areas into account, because it is likely that a central business district has different profiles than a sparse suburban area (e.g., Badas et al. 2019). To model this, we attempted to use the mean building height, standard deviation of building height, plan area density λ p and frontal area density λ f (Grimmond and Oke 1999) by trying to find possible scaling relationships that improve the collapse of the profiles. Unfortunately, these efforts did not significantly improve the scatter in Fig. 2 when using the datasets of Sect. 3. We did find an estimate for the blending height z * = z h + 1.5σ h , which is used in scaling the height of the w 3 profile (see Sect. 4.4). Perret et al. (2019) were similarly unsuccessful in collapsing σ 2 u by scaling, although they did not use spatially averaged profiles. On the other hand, the Rotach (2001) RSL parametrizations successfully use the zero-plane displacement d, z * , and z 0 , which depend on the building geometry (Grimmond and Oke 1999). Since the UCL profiles are continuous to the RSL profiles and thus also depend on those factors, they also indirectly affect the UCL parametrizations. It is possible that more data or more realistic data would help in finding appropriate scaling parameters for differing geometries.
An additional limitation of the chosen approach is missing the multi-scale roughness, because the obstacles in wind-tunnel studies or simulations are most often aerodynamically smooth, in contrast to real buildings. However, Vanderwel and Ganapathisubramani (2019) indicated that small-scale roughness elements have a negligibly small effect on drag.

Verification
In this section the model's basic assumptions are verified and its general sensitivities are investigated.

Well-Mixed Criterion
The most fundamental test for an LPDM is the well-mixed criterion after Thomson (1987). If a model fulfills this criterion, it does not unmix once the particles are well-mixed in physical and velocity space. To test whether the modified model fulfills this criterion, we initialize the ULM with particles uniformly distributed in height and given a velocity randomly sampled from the local velocity density function. Then the ULM simulates the dispersion for an extended amount of time (in our case about 1 day of simulated time).
The resulting height distribution of particles-similar to, e.g., Bahlali et al. (2020)-is shown in Fig. 5 for the three test cases in Table 2. The normalized concentration in Fig. 5 is not calculated using the kernel method (de Haan 1999) like in the rest of this work. To avoid any smoothing, we simply calculated a histogram over height-bins and normalized the particle count by the the mean particle count in a perfectly uniform distribution. Note that the figure only shows the lower 7.5% (stable), 1.5% (neutral), and 0.8% (convective) of the model domain. Aloft and up to the boundary-layer height z i the value of the normalized  Table 2, initialized with well-mixed particles. The particle counts are binned (bin size = 1 m) and normalized with their theoretical value for each height range. A perfect result would be 1 (dotted line) everywhere. The figure shows only the lower part of the model domain concentration fluctuates around 1, as expected (not shown). Below the mean building height z h , the normalized concentration of particles is calculated only for the fluid volume, analog to intrinsic spatial averages in Sect. 3. This necessitates normalizing the particles by their mean count (as above) and then multiplying by 1 − λ p . After this, the normalized concentration fluctuates around 1, proving the well-mixed state.

Comparison of Urban-Canopy Layer Model and Roughness Sublayer Model
To investigate the impact of the UCL on the simulated concentration fields, we use three synthetic cases taken from Kljun et al. (2015). The relevant values are listed in Table 2. Note that the source height is 1 m above the mean building height. A source height in the street canyon might possibly produce more pronounced differences between the two models. However, as RSM is restricted to heights z > d, we decided to use a source height representative for domestic heating ('chimney height'). Figure 6 shows selected vertical profiles of the crosswind-integrated concentration CIC, the maximum value along an arc ARCMAX and the standard deviation of particle spread in the lateral direction σ y . Note that we use the term 'arcs', but utilize a Cartesian grid, so that the 'arc' strictly speaking refers to a line perpendicular to the mean wind direction. We use the term 'arc' for historical reasons. The model can calculate the concentration at any point in the model domain independent from other points, so the shape of the grid is not important, as long as the resolution is not too low. Also note that the source distances shown vary between the three stability conditions. In all three conditions-stable, neutral, and convective-the . 6 Vertical profiles of CIC (blue), ARCMAX (orange), and σ y (green) for three different stabilities (rows, see Table 2) and different source distances (x, note that the source distances shown vary between the three stability conditions  Note that the mean building height z h = 15 m, zero-plane displacement d = 10 m, z * z h = 1.5, roughness length z 0 = 1 m and source height z ini = 16 m for all three scenarios peak value of the CIC near the source at approximately roof level is smaller for the ULM (Fig. 6a, d, g). This is not offset by increased vertical dispersion, because the CIC values aloft are highly similar and those below do not compensate either. Since the lateral dispersion measured in σ y is similar for both models at these heights, the reduced dispersion of the ULM must be in the longitudinal direction. This is possible even though the model is horizontally homogeneous, because there is a stochastic effect on u and the mean transport changes between the two models, due to changes in domain size and wind speed profile.
The quantity σ y sometimes displays artifacts when too few particles are available at certain heights. This can be seen in Fig. 6a, b, d, at heights where σ y differs for the RSM and the ULM. CIC and ARCMAX are not affected, because their calculation is not as susceptible to few and strongly different values.
At some distance from the source (Fig. 6b, e, h), the roof level CIC and ARCMAX are smaller with the ULM than with the RSM. For the largest distances downstream, the concentrations of the two model versions start to overlap (Fig. 6c, f) due to approaching a well-mixed state.
Another important result is the slanted profile of CIC and ARCMAX in the UCL for ULM in nearly all stabilities and distances, at least until the whole profile starts to be wellmixed. In contrast, the RSM must make an assumption about the profile below d and this was historically a well-mixed assumption. This is shown in Fig. 6 by extending the profile of RSM uniformly down to the ground, as dotted lines. However, the slanted profiles of ULM in the UCL indicate that the assumption of a uniformly distributed concentration in the UCL does not hold.

Sensitivity Analysis
The chosen UCL parametrizations are based on data with considerable spread. A sensitivity study is therefore conducted in order to assess the impact the newly introduced profile parameters ('UCL parameters' in the following) (see Sect. 4) exhibit on the output (i.e., the modelled concentrations). Due to the non-negligible computational effort per model run, the state-of-the-art approach, i.e., a variance decomposition method, is unfeasible. Instead, we follow Saltelli et al. (2008) and choose 'Morris sampling', after Morris (1991), enhanced by Campolongo et al. (2007). The method is also called the 'elementary effects method' and employs one-factor-at-a-time changes (OAT), but alleviates many of the disadvantages of OAT-based methods (Saltelli et al. 2008). Unfortunately, Morris sampling requires scalar outputs, while in the present case we have concentration distributions. We therefore consider four scalars describing the concentration distribution.
The four scalars are shown in Fig. 7: CIC , the maximum of the crosswind-integrated concentration; width , the width in metres of the CIC peak, defined as the distance between with increasing source distances. The four blue labels describe four scalars characterizing the plume the points where the linearly interpolated CIC curve reaches 75% of CIC ; σ y , the standard deviation of particle spread in the lateral direction at the distance of CIC occurence; and dσ y , the derivative of σ y with respect to the streamwise direction x. We test all UCL parameters, as described in Sect. 4. These are a um , a Re , a u2 , a v2 , a w2 , a Sk , a f , and a ε . Furthermore, we vary λ p from 0 to 1, which governs the reflection at roof level from 0 to 100%. For consistency with the other parameters, it is denoted by a lp in the following. Additionally we also vary some parameters of the original RSL parametrization('RSL parameters' in the following), to give context to the new parameters. These are a d , which varies the zero-plane displacement d from 0.5 to 0.9 of z h ; a zs , which varies z * z h from 1.2 to 5.0; and a a , which is the parameter a in Eq. 1 in Rotach (2001) and determines the shape of the RSL parametrization of the local u * ,l .
For details regarding the implementation of the Morris method, see "Appendix 2". Basically, it returns a summary statistic μ * (Eq. 27 in the "Appendix 2"), which measures the impact a parameter has on the overall model output; lower means less impact. In addition to μ * we used arithmetically averaged μ * to jointly asses the impact of (i) the aforementioned four output scalars ( CIC , width , σ y , dσ y ); (ii) four heights where the concentrations are calculated (1 m, 16 m, 31 m, 46 m), 1 m above ground and roof, equidistant above; (iii) three meteorological situations (see Table 2); (iv) three source heights (2 m, 16 m, 40 m), deep within the canyon, above z h and in the upper part of the RSL (depending on a zs . Figure 8 singles out the three stability cases in the first three panels and combines all cases in the fourth. Generally speaking, the RSL parameters (in lighter colours) have a slightly higher impact on the output than the UCL parameters. This result is fortunate, because the UCL parameters are highly uncertain, due to the large spread of the profiles (Fig. 2). The most important individual parameter is a zs , with the highest impact in the neutral and especially Fig. 8 Aggregated impact, μ * (higher value signify higher influence), of model parameters according to the Morris method for the three stability groups (blue bars) and all combined (orange bars). UCL parameters in dark colours; RSL parameters in lighter colour. Vertical black lines indicate one standard deviation the stable cases and still considerable impact in the convective case. Already Rotach (2001) noted the importance of this parameter and showed that it is less impactful in convective conditions. With the explicit treatment of the UCL, the impact of a d becomes less important than that of a zs , while in the RSM (Rotach 2001) it had been similar. In contrast to a zs , a d has a stronger impact when stability decreases. The RSL parameter a a is most impactful in the neutral case and less so for stable and convective situations. The convective case is especially interesting, because there a a is even less impactful than some of the UCL parameters. This is in spite of a a governing all turbulence profiles in the RSL and-due to continuity at the roof level-therefore also all turbulence profiles in the UCL. However, in convective situations the impact of u * on the RSL parametrizations is diminished in the model, because they additionally depend on w * . The impact of a lp is not overly large, despite a lp 's property of varying the roof-top reflection from 0% to 100%, thereby in the extreme cases even turning the UCL parametrizations on or off. We suspect that this is caused by the plume not dispersing strongly in the stable and neutral cases and therefore not many particles actually reaching the measurement layers far away from the source height. This is corroborated by the high importance in the convective case with much higher vertical dispersion. See also the discussion of Fig. 9.
Of the new parameters, a um and a v2 are the most important. On one hand, a um governs the mean flow thus explaining its importance. On the other hand, the relatively large impact of a v2 is not due to us choosing two scalars describing the lateral dispersion ( σ y and dσ y ), because removing dσ y does not change the relative importance of a v2 much (not shown). Of medium importance among the UCL parameters are a u2 , a w2 , and a f . Since the latter is not based on data (Sect. 4), its medium impact calls for more work with respect to the suitability of describing the velocity p.d.f. in the UCL using the function f (Eq. 1) in combination with the vertical velocity skewness. a Re , a Sk and a ε are of least importance, so that the lack of data for at least a Sk and a ε (Fig. 2) does not add a major source of uncertainty for the UCL simulations. Among the two parameters influencing the Reynolds stress profile, a Re as an UCL parameter is of minor importance, while its RSL counterpart (a a ) is among the most influential. This is because in the RSL parametrization Rotach (2001) all turbulence profiles depend on u * ,l , while in the UCL parametrization of this study, they do not. Consequently, a Re only modifies a profile that rapidly approaches zero with decreasing height and its effect is apparently small. This is further amplified by the relatively small spread of the profiles in Fig. 3b, compared to that for, e.g., the velocity variances in Fig. 3e, f.
When we separate the results of the sensitivity study for the different emission heights (Fig. 9), we see a different pattern. Note that the 'combined' panel is identical to that of Fig. 8. Clearly, some parameters affect model runs with different emission height differently. It is interesting to note that the UCL parameters in the runs with a high source affect the output similarly than in runs with a lower source, despite only affecting the layer below 15 m. The high importance of a zs and the other RSL parameters is due to their impact on the entire profiles over the buildings, which is the majority of the model domain. The parameter a lp shows that the reflection at the roof level is less impactful if the particles already start below this level. Also interesting is a much higher separation between important and less important parameters for lower sources than for higher sources. This is likely related to the fact that in those cases the UCL parametrization directly influences all particles at the start, instead of only when they are transported into the UCL.

Validation
To properly test the UCL parametrization, the model performance is evaluated for a real dispersion experiment. Required input is mean wind speed at an arbitrary height, mean wind direction, u * , w * (if convective), z i , Obukhov length L, mean building height, plan area density λ p , zero-plane displacement d, roughness length z 0 , and the blending height z * . For the validation, the coordinates and strength of the source and the positions and concentration measurements of the samplers are also required. Here we use data from the Basel UrBan Boundary Layer Experiment (BUBBLE, Rotach et al. 2005). We chose BUBBLE, because we wanted a tracer experiment that approximates the steady state of the model and has a source as low as possible, but not below the zero-plane displacement d. Otherwise the comparison with the RSM would be compromised, because it cannot start particles below d. Gryning et al. (2005) describes the BUBBLE campaign in Basel, Switzerland, from 2001 to 2002. This campaign contains both long-term tower measurements as well as tracer experiments (see also Sect. 3). The tracer experiments took place during four convective days with a thermally driven, local wind phenomenon. Release and most measurements were near roof level. Mean building height in the area is about 15 m. Rotach et al. (2004) used the same LPDM to compare with the BUBBLE dataset and investigated three different averaging periods for the concentration measurements. In order to have the largest possible amount of data for the analysis, we choose the original measurement period of 30 min for our simulations, instead of the 3-h average that Rotach et al. (2004) chose. This decreases the statistical agreement between model and measurements, but gives us more data to work with. Meteorological data were collected at a tower with sonic anemometers. The tower stood within the area covered by the samplers for the concentration measurements. All together, 24 30-min convective periods are available for analysis.

Acceptance Criteria After Hanna and Chang (2012)
Hanna and Chang (2012) provide acceptance criteria for urban dispersion model evaluations. They define the fractional bias (FB), the normalized mean square error (NMSE), and the factor-of-2 (F2). See "Appendix 3" for their definitions. Note that these error measures for the acceptance criteria are not calculated for all points, but only for the maximum value along each arc (ARCMAX), which leads to far better statistical agreement than point-topoint comparisons. The measurements were taken from the highest value along each arc during BUBBLE (see Rotach et al. 2004). Furthermore, Hanna and Chang (2012) define the threshold-based normalized absolute difference (TBNAD, called NAD in their work) for point-to-point comparison if both simulated and observed values are above the threshold of three times the limit of quantification of the measuring instrument. The so-defined statistics are summarized in Table 3.
The RSM already shows a strong performance, fulfilling all four criteria when aggregated over all BUBBLE cases. Only the TBNAD-criterion shows non-optimal performance. When looking at the performance measures for each case individually (last column), the FB-, NMSE-, and F2-criteria are accepted in well above the 50% of cases required by Hanna and Chang (2012). Conversely, for both model versions, the TBNAD-criterion is only accepted in 38% of cases, confirming the earlier, worse model performance when looking at TBNAD. Including the UCL does not change the model performance; in general a slight improvement is found.

Visual Comparison with Field Measurements
Figure 10 compares measured and modelled concentrations, both for the RSM (blue crosses) and the ULM (orange squares). Due to the logarithmic nature of the plot and the fact that the model simulates values near zero to floating point precision-far smaller than the experimental detection limit-the scale is limited to the instrument's limit of detection (LOD = 5 ng m −3 ). All values smaller than LOD (simulated or measured) are set to LOD. Generally, the ULM included produces slightly smaller concentrations than the RSM. This will be shown more clearly in Sect. 6.4. However, the distribution of values is highly similar, indicating that the implementation of the UCL does not fundamentally alter the concentration distribution-at least above roof level where the measurements were performed.
To the side and the top of Fig. 10 are marginal histograms of the distributions. Note that the bins of the histogram follow the same logarithmic scale as the scatter plot. These distributions do not correspond to any theoretical expectation, since they are mostly governed by the choice of measurement locations during the BUBBLE campaign. Both models (right panel in Fig. 10) fail to reproduce the local peak in the observations (top panel in Fig. 10) at 3× 10 −2 mg m −3 . At a large number of sites, where this most frequently observed concentration is measured, both models return a value below the LOD (note the roughly doubled number of LODs in the simulations, as compared to the observations). To quantify the difference between the two simulated distributions and the measured distribution, we use the earth mover's distance (EMD, calculated after Werman 2008, 2009). It is a measure for the distance between two distributions. A larger value indicates that the difference between two distributions is larger. When we calculate the EMD for each model's distribution compared to the measured distribution, the ULM has a slightly lower value, indicating a slightly better performance.

Statistical Comparison with Field Measurements
Measured and simulated data points are co-located in space and time, even if co-location in time might be debatable since the LPDM assumes stationarity. We use the Relative Difference (RD), the Fractional Bias (FB), the Normalized Mean Square Error (NMSE), the correlation coefficient (CORR), the Factor of 2 (F2), the Factor of 5 (F5), the Bounded Normalized Mean Square Error (BNMSE), and the Threshold-Based Normalized Absolute Difference (TBNAD), all of which are detailed in the Appendix 3, Eqs. 28-34. We use a bootstrapping procedure to judge significance of these statistical measures, also detailed in Appendix 3. Table 4 shows the error statistics for the simulation of the BUBBLE tracer experiments with the two model versions. The ULM shows better values for all measures and roughly half the improvements are statistically significant (98% level). The RSM overpredicts the concentrations (negative FB). Since the ULM delivers slightly lower concentration than the RSM, this improves the statistical agreement due to reduced overprediction. The NMSE,

Conclusion and Outlook
An existing Lagrangian particle dispersion model (called RSM) with parametrized flow and turbulence characteristics for the urban roughness sublayer (RSL), and model domain starting at the zero-plane displacement, was extended down into the urban canopy layer (UCL). The model with the UCL included is called ULM. For this, spatially averaged profiles of flow and turbulence were collected from the literature, mostly from wind-tunnel and numerical studies, but also including a few full-scale experiments. The profiles show considerable spread, but nevertheless parametrizations for the mean wind speed u; Reynold's stress u w ; dissipation rate of turbulence kinetic energy ε; skewness of the vertical wind component Sk w ; velocity variances u 2 , v 2 , and w 2 ; and a model-specific transition function f are proposed. These parametrizations are continuous to the parametrizations in the RSL aloft and are intrinsic in nature, i.e., only averaged over the fluid volume-excluding the buildings. For this reason a reflection of probability λ p is introduced at the mean building height. The parametrizations are horizontally homogeneous and independent of wind direction. The ULM still fulfills the well-mixed criterion, with the peculiarity in the UCL that the normalized particle concentration assumes the value of 1 − λ p instead of 1, which is due to the particle reflection mentioned before.
A sensitivity analysis using the Morris method shows that the parameters describing the UCL turbulence have a smaller impact on the concentration distribution and plume characteristics than have those parameters which characterize the RSL turbulence profiles. This reduces the uncertainty introduced through the relatively large spread of the parametrized UCL turbulence profiles.
A comparison of the ULM and the RSM shows slight changes in the model output, but generally the same dispersion behaviour. The longitudinal dispersion is somewhat diminished by the inclusion of the UCL parametrizations. An immediate corollary of this finding is the result that the concentration profiles in the UCL appear to be non-constant with height. This means that assuming uniformly height-distributed concentrations below the zero-plane displacement is questionable.
The RSM and the ULM both fulfill the acceptance criteria after Hanna and Chang (2012). When tested against concentration measurements of the BUBBLE dataset, the ULM shows slightly improved concentration predictions compared to the RSM. Some of the improvements in the error measures are statistically significant, but the overall differences are relatively small.
In summary, the analysis has demonstrated that spatially-averaged canopy turbulence can be introduced into an LPDM without violating the well-mixed condition. Due to the fact that the spatial variability of turbulence characteristics is considerable and hence spatial averages bear a huge uncertainty range, application of the model for near-surface sources to determine near-surface (pedestrian level) concentrations will result in correspondingly uncertain point predictions. Comparison to data from an above-canopy urban tracer experiment, conversely, shows that modelled concentrations are not strongly affected (even slightly improved) when allowing for the UCL to be explicitly included in the model domain. This to some degree reflects the results from the Morris sensitivity analysis, which has demonstrated that the uncertainty of the concentration estimates are to a lesser degree affected by the canopy turbulence parameters than those describing the bulk urban fabric. Thus, the ULM is likely more suitable for modeling dispersion from average street level sources than the RSM.
Furthermore, extending the model domain down to the physical surface will make it feasible to use the dispersion model as a core for a footprint model. Footprints are often required for pedestrian or surface level and hence sensitive to near-surface turbulence. This makes a dispersion model with explicit treatment of UCL turbulence particularly suited for this purpose. However, they are not easily validated directly, so the present results add support to the validity of using the ULM for future footprint modeling applications over urban areas.
Source code of the LPDM, scripts used to normalize, scale and fit the UCL parametrizations, as well as scripts used to analyze the model output are available (Stöckl 2021). The datasets used to drive the LPDM model are available (Stöckl 2021). The UCL profile datasets taken from other studies are available from the corresponding author on reasonable request. version of the w -p.d.f. called where A and B stand for up-and down-drafts, respectively (Luhar 1989). They are defined after Baerentsen (1984) through their p.d.f.s where with q = 1 √ 1−ρ 2 and ρ the correlation coefficient between the Gaussian u and w after Rotach et al. (1996), modified from Luhar (1989). Note that these formulations were only ever meant to model positive skewness for convective conditions. It is, however, possible to use the same formulations to achieve a negatively skewed p.d.f. for w if w 3 in Eq. 23 is negative. Compared to a positive value of w 3 , σ B in Eq. 23 increases, while σ A in Eq. 24 decreases when w 3 is negative. This also leads to a smaller B and larger A in Eq. 25. σ A and σ B are both the scale and location parameter of the corresponding p.d.f. P A and P B , respectively (note the reversed sign of the location parameter in P B ). Altogether, this means that there are now smaller areas of stronger downdrafts and larger areas of weaker updrafts in Eq. 20. Consequently, the p.d.f. of w is skewed negatively.

Appendix 2: Morris Method of Global Sensitivity Analysis
This section provides a short overview of the Morris or elementary effects method (EE) used in the sensitivity study, as well as implementation details. The method requires a scalar output and Sect. 5.3 describes how to derive scalar output from the LPDM. Furthermore, the method requires deterministic output, because noise from random output errors is difficult to distinguish from small effects. However, our model-and therefore its output-is stochastic in nature, not deterministic. To alleviate this, we increase the number of particles and thus reduce the random error of the output. Daniel (1973) states that one-at-a-time (OAT) methods, which are similar but less powerful than the Morris method, can only show effect sizes that are larger than four times the standard deviation of the output. To test this, we ran the model 100 times using identical settings for various numbers of particles and calculated the standard deviation of the scalar outputs. The standard deviation decreases with increasing number of particles (not shown), which proves that given enough particles, it is possible to distinguish effect sizes that are larger than the stochastic variation. We chose enough particles to make this matter of no concern, especially since the bulk of our analysis is based on an aggregate measure that does not take the spread into account (see Eq. 27). Morris (1991) explains how to distinguish between three types of effects a model parameter can have on the output: (i) negligible, (ii) linear and additive, and (iii) nonlinear or involved with interactions with other parameters. To find out into which each of the k-many input parameters belongs, an elementary effect d i (X) = ⎧ ⎨ ⎩ y(X 1 ,...,X i−1 ,X +Δ i ,X i+1 ,...,X k )−y(X) Δ 1 σ y when increasing the i-th parameter y(X)−y(X 1 ,...,X i−1 ,X −Δ i ,X i+1 ,...,X k ) Δ 1 σ y when decreasing the i − th parameter (26) is calculated from a pair of scalar model outputs y that depend on the vector of input parameters X, where the i-th parameter was changed by stepping up or down by Δ in its distribution and σ y is the standard deviation of the model output y. Depending on the parameter's direction of change (step-up or step-down), the subtraction is flipped, as shown in Eq. 26 (Saltelli et al. 2008). The distribution of the parameters are assumed as uniform and linearly transformed to the interval [0, 1] with p many discrete positions. We choose p = 6, thus the step size Δ = p 2( p−1) = 3 5 after Morris (1991). Saltelli et al. (2008) give a visual explanation why this Δ leads to uniform coverage of the parameter space. The normalization by σ y is not from Morris (1991), because he only looked at one output scalar y in this work. We, however, want to compare more than one model output and the normalization increases the comparability of the elementary effects of each of the output-types (Sin and Gernaey 2009). Note that unlike Sin and Gernaey (2009), we do not normalize Δ, because we transformed each parameter range to the interval [0, 1] as described by Morris (1991). Since the actual value of Δ does not directly influence the output y, only which discrete point in the parameter distribution is chosen, this is not necessary, as long as Δ remains unit free.
To calculate more than one d i , naively, one could now create a large number of pairs with different input parameter vectors X and compare two at a time. However, the novel concept of Morris (1991) is to randomly generate a 'trajectory' in the k-dimensional parameter space, where each following step changes exactly one parameter and each parameter is only changed once, such that a trajectory has k +1-many points (see Saltelli et al. (2008) for an illustration). This means that only k + 1 simulations lead to k-many d i , one per input parameter. Then the same procedure is repeated r times with different (random) trajectories, until there are r -many d j i and simple statistics over the distribution of d j i can be calculated. For the exact method of generating the trajectories see Morris (1991), except that we also allow a stepup of Δ (instead of only step-down), as implemented in the R-package 'sensitivity' (Iooss et al. 2019). Note that we use an additional improvement over the original method of Morris (1991) and generate 1000 trajectories before selecting the r = 20 with the largest summed Euclidean distance between the selected trajectories for increased coverage of the parameter space (Campolongo et al. 2007). The amount of 20 trajectories is fairly arbitrary and was chosen as a compromise between computational effort and quality of the result. See Khare et al. (2015) for a short review on the number of trajectories used in literature, as well as a study on the suitability of the Campolongo et al. (2007) method to select the trajectories (called OT in Khare et al. 2015).
The arithmetic mean μ i and the sample standard deviation σ i of the d j i distributions (for each i = 1 . . . k) are the classic Morris sensitivity measures and have to be looked at simultaneously. Campolongo et al. (2007) suggest the use of as a summary statistic with some loss of information, but increased ease of use.

Appendix 3: Statistical Treatment
We use the following definitions to measure error statistics: where o is the observed concentration and s the simulated concentration, both observed at a receptor. A F is the number of false-positive and false-negative predictions, where both model and measurements are judged on whether or not they exceed a threshold, defined as three times the measurement device's limit of detection (LOD). A OV is the number of OVerlapping predictions, meaning where both or none of the values exceed the threshold.
To be able to judge the significance of differing statistical measures, we use a bootstrapping procedure (Tukey 1987), where we calculate the statistical measure for 50,000 random samples with replacement of the simulation-observation pair vectors. These 50,000 statistical measures form the bootstrap distribution. The exact procedure is listed in Stöckl (2015), but generally we follow the 'moment bootstrap' of Hanna (1989). Additionally, we do not compare the confidence intervals of the two bootstrap distributions directly, but instead bootstrap the difference of the error measure between the two model variants for the 'same sample' (Irwin 2001). The 'same sample' means that we randomly draw rows of the two simulation and the one observation vectors. If the 98% confidence interval of the difference between the two comparative error measures does not contain the value zero, the difference in the error measure is significant.