Joint PDF modelling of turbulent flow and dispersion in an urban street canyon
Abstract
The joint probability density function (PDF) of turbulent velocity and concentration of a passive scalar in an urban street canyon is computed using a newly developed particleincell Monte Carlo method. Compared to moment closures, the PDF methodology provides the full onepoint onetime PDF of the underlying fields containing all higher moments and correlations. The smallscale mixing of the scalar released from a concentrated source at the street level is modelled by the interaction by exchange with the conditional mean (IECM) model, with a micromixing time scale designed for geometrically complex settings. The boundary layer along noslip walls (building sides and tops) is fully resolved using an elliptic relaxation technique, which captures the high anisotropy and inhomogeneity of the Reynolds stress tensor in these regions. A less computationally intensive technique based on wall functions to represent boundary layers and its effect on the solution are also explored. The calculated statistics are compared to experimental data and largeeddy simulation. The present work can be considered as the first example of computation of the full joint PDF of velocity and a transported passive scalar in an urban setting. The methodology proves successful in providing high level statistical information on the turbulence and pollutant concentration fields in complex urban scenarios.
keywords:
Langevin equation; MonteCarlo method; Probability density function method; Scalar dispersion; Urbanscale turbulence1 Introduction
Regulatory bodies, architects and town planners increasingly use computer models to assess ventilation and occurrences of hazardous pollutant concentrations in cities. These models are mostly based on the Reynoldsaveraged NavierStokes (RANS) equations or, more recently, largeeddy simulation (LES) techniques. Both of these approaches require a series of modelling assumptions, including most commonly the eddyviscosity and gradientdiffusion hypotheses. The inherent limitations of these approximations, even in the simplest engineering flows, are well known, and detailed for example by Pope (2000). Although the effects of the assumptions are more pronounced in RANS than in LES models, where they are confined to the smaller (modelled) scales, there is clearly a need to develop higherorder models. In pollutant dispersion modelling it is also desirable to predict extreme events such as peak values or probabilities that concentrations will exceed a certain threshold. In other words, a fuller statistical description of the concentration field is required (Chatwin and Sullivan, 1993; Kristensen, 1994; Wilson, 1995; Pavageau and Schatzmann, 1999). These issues have been explored in grid turbulence and in the unobstructed atmosphere, and models capable of predicting higherorder statistics have also appeared (Yee et al., 1994; Yee and Wilson, 2000; Luhar et al., 2000; Cassiani and Giostra, 2002; Franzese, 2003; Cassiani et al., 2005a, b), but more research is necessary to extend these capabilities to cases of builtup areas.
Probability density function (PDF) methods have been developed mainly within the combustion engineering community as an alternative to moment closure techniques to simulate chemically reactive turbulent flows (Lundgren, 1969; Pope, 1985; Dopazo, 1994). Because manyspecies chemistry is highdimensional and highly nonlinear, the biggest challenge in reactive flows is to adequately model the chemical source term. In PDF methods, the closure problem is raised to a statistically higher level by solving for the full PDF of the turbulent flow variables instead of its moments. This has several benefits: convection, the effect of mean pressure, viscous diffusion and chemical reactions appear in closed form in the PDF transport equation. Therefore these processes are treated mathematically exactly without closure assumptions, eliminating the need for gradienttransfer approximations. The effects of fluctuating pressure, dissipation of turbulent kinetic energy and smallscale mixing of scalars still have to be modelled. The rationale is that, since the most important physical processes are treated exactly, the errors introduced by modelling assumptions for less important processes amount to a smaller departure from reality. Moreover, the higher level description provides more information that can be used in the construction of closure models.
The PDF transport equation is a highdimensional scalar equation. Although techniques of solution based on stochastic Eulerian methods have been developed (Valiño, 1998; Mobus et al., 2001; Soulard and Sabel’nikov, 2006), Lagrangian Monte Carlo methods are a more natural choice because their computational cost increases only linearly with the dimensionality of the problem, favourably comparing to the more traditional finite difference, finite volume or finite element methods.
The numerical development of Lagrangian PDF methods has mainly centred around three distinctive approaches, all of them representing a finite ensemble of fluid particles with Lagrangian particles. A common approach is the standalone Lagrangian method, where the flow is represented by particles whereas the Eulerian statistics are obtained using kernel estimation (Pope, 2000; Fox, 2003). Another technique is the hybrid methodology, which builds on existing Eulerian computational fluid dynamics (CFD) codes based on moment closures (Muradoglu et al., 1999, 2001; Jenny et al., 2001; Rembold and Jenny, 2006). Hybrid methods use particles to solve for certain quantities and provide closures for the Eulerian moment equations using the particle/PDF methodology. A more recent approach is the selfconsistent nonhybrid method (Bakosi et al., 2007, 2008), which also employs particles to represent the flow, and uses the Eulerian grid only to solve for inherently Eulerian quantities (such as the mean pressure) and for efficient particle tracking. Since the latter two approaches extensively employ Eulerian grids, they are particleincell methods (Grigoryev et al., 2002).
The current study presents an application of the nonhybrid method to a simplified urbanscale case where pollution released from a concentrated line source between idealized buildings is simulated and results are compared to data from windtunnel experiments and LES.
PDF methods in atmospheric modelling mostly focus on the simulation of passive pollutants, wherein the velocity field (mean and turbulence) is either assumed or obtained from experiments (Sawford, 2004, 2006; Cassiani et al., 2005a, b, 2007). In contrast, the current model directly computes the joint PDF of the turbulent velocity, characteristic turbulent frequency and scalar concentration, extending the use of PDF methods in atmospheric modelling to represent more physics at a higher statistical level. A computed full joint PDF also has the advantage of providing information on the uncertainty originating from turbulence on a physically sound basis.
In this study the turbulent boundary layers developing along solid walls are treated in two different ways: either fully resolved or via the application of wall functions (i.e. the logarithmic “law of the wall”). The full resolution is obtained using Durbin’s elliptic relaxation technique (Durbin, 1993), which was incorporated into the PDF methodology by Dreeben and Pope (1997a, 1998). This technique allows for an adequate representation of the nearwall lowReynoldsnumber effects, such as the high inhomogeneity and anisotropy of the Reynolds stress tensor and wallblocking. Wall conditions for particles based on the logarithmic “law of the wall” in the PDF framework have also been developed by Dreeben and Pope (1997b). These two types of wall treatments are examined in terms of tradeoff between computational cost and performance, addressing the question of how important it is to adequately resolve the boundary layers along solid walls in order to obtain reasonable scalar statistics.
At the urban scale the simplest settings to study turbulent flow and dispersion patterns are street canyons. Due to increasing concerns for environmental issues and air quality standards in cities, a wide variety of canyon configurations and release scenarios have been studied both experimentally (Hoydysh et al., 1974; Wedding et al., 1977; Rafailids and Schatzmann, 1995; Meroney et al., 1996; Pavageau and Schatzmann, 1999) and numerically (Lee and Park, 1994; Johnson and Hunter, 1995; Baik and Kim, 1999; Huang et al., 2000; Liu and Barth, 2002). Street canyons have a simple flow geometry, they can be studied in two dimensions, and a wealth of experimental and modelling data are available for different streetwidth to buildingheight ratios. This makes them ideal candidates for testing a new urban pollution dispersion model. We validate the computed velocity and scalar statistics with the LES results of Liu and Barth (2002) and the windtunnel measurements of Meroney et al. (1996); Pavageau (1996); Pavageau and Schatzmann (1999). The experiments have been performed in the atmospheric wind tunnel of the University of Hamburg, where the statistics of the pollutant concentration field have been measured in an unusually high number of locations in order to provide fine details inside the street canyon.
The rest of the paper is organized as follows. In Section 2, the exact and modelled governing equations are presented. Several statistics are compared to experimental data and largeeddy simulations in Section 3. Finally, Section 4 draws some conclusions and elaborates on possible future directions.
2 Governing equations
We write the Eulerian governing equations for a passive scalar released in a viscous, Newtonian, incompressible flow as
(1)  
(2)  
(3) 
where , , , , and are the Eulerian velocity, pressure, constant density, kinematic viscosity, scalar concentration and scalar diffusivity, respectively. Based on this system of equations the exact transport equation that governs the onepoint, onetime Eulerian joint PDF of velocity and concentration can be written as (Pope, 1985, 2000),
(4) 
where and denote the sample space variables of the stochastic velocity and concentration fields, respectively, and the pressure is decomposed into its mean and fluctuation part . In Equation (4) the physical processes of advection (second term on the left), viscous diffusion (first term on the right) and transport of in velocity space by the mean pressure gradient (second term on the right) are represented mathematically exactly. The last three terms in the form of conditional expectations have to be modelled: these are respectively the effect of dissipation of turbulent kinetic energy, the effect of fluctuating pressure and the smallscale diffusion of the scalar. After appropriate modelling of the unclosed terms, Equation (4) can, in principle, be solved with a traditional numerical method. However, the highdimensionality makes Monte Carlo methods more appealing. In particular, because Equation (4) is a FokkerPlanck equation, it can be written as an equivalent system of stochastic differential equations (SDEs) (van Kampen, 2004). We use the generalized Langevin model (GLM) of Haworth and Pope (1986) for the velocity increment, and the interaction by exchange with the conditional mean (IECM) model for the scalar. The physics and characteristics of the IECM model are discussed in detail by Fox (1996); Pope (1998) and Sawford (2004). Thus our system of SDEs that solve Equation (4) is written as
(5)  
(6)  
(7) 
Equation (5) governs the Lagrangian particle position , and it consists of advection by the instantaneous particle velocity and molecular diffusion represented by the isotropic Wiener process , which is a given Gaussian process with zero mean and variance . The particle velocity is governed by Equation (6). The second and third terms on the righthand side are a direct consequence of the particular Lagrangian representation of the viscous diffusion in Equation (5) (for a derivation, see Dreeben and Pope, 1997a), thus the first three terms on the righthand side govern the particle velocity in a laminar flow mathematically exactly. The last two terms involving the tensor and jointly model the effect of pressure redistribution and anisotropic dissipation of turbulent kinetic energy. is a secondorder tensor function of the mean velocity gradients , the Reynolds stresses and the rate of dissipation of turbulent kinetic energy , while is a positive constant. Note that the Wiener process appearing in both Equations (5) and (6) is the same process, i.e. the same exact series of random numbers, and is independent of the other process . The concentration of the passive scalar is governed by Equation (7), which represents the physical process of diffusion by relaxation of the particle scalar towards the velocityconditioned scalar mean with a timescale . Note that the Eulerian statistics denoted by angled brackets are to be evaluated at the fixed particle locations , and in particleincell methods, this is usually achieved using an Eulerian grid and computing ensemble averages in grid elements and/or around vertices. The energy dissipation rate is defined as
(8) 
where is a model constant. We adopt the model of van Slooten et al. (1998) for the stochastic characteristic turbulent frequency
(9) 
where is a source/sink term for the mean turbulent frequency
(10) 
where is the production of turbulent kinetic energy, is a scalar valued Wiener process, while and are model constants. To define the micromixing time scale for a scalar released from a concentrated source in a geometrically complex flow domain we follow Bakosi et al. (2007, 2008) and specify the inhomogeneous as
(11) 
where denotes the radius of the source, is a characteristic velocity at that we take as the Euclidean norm of the mean velocity vector, is the location of the source, and and are model constants.
2.1 Elliptic relaxation modelling of nearwall turbulence
The turbulence model GLM (represented by the last two terms of Equation (6) for the velocity increments) is in fact a family of models. A specific definition of and corresponds to a particular closure from a wealth of models, many of which can be made equivalent (at the level of second moments) to popular Reynoldsstress closures (Pope, 1994). To be able to capture the nearwall lowReynoldsnumber effects on the Reynolds stresses in fully resolved boundary layers, we follow Durbin (1993) and Dreeben and Pope (1998) and specify and through the tensor
(12)  
(13) 
where denotes the turbulent kinetic energy and is obtained by solving the elliptic equation
(14) 
where the fourthorder tensor is given by
(15) 
with
(16) 
the Reynolds stress anisotropy
(17) 
and are model constants. The characteristic length scale in Equation (14) is defined by the maximum of the turbulent and Kolmogorov length scales
(18) 
with , where and are model constants, and is the unit wallnormal of the closest wall element pointing outward of the flow domain. For the applied wallboundary conditions in this fullyresolved case the reader is referred to Dreeben and Pope (1998) and Bakosi et al. (2008). More details on the inflow and outflow conditions for the mean pressure and on the wall conditions for the tensor are given in (Bakosi et al., 2008). Equation (14) was developed in conjunction with turbulent channel flow (Durbin, 1993; Dreeben and Pope, 1998). Modifications and different forms of this idea have also been proposed (Whizman et al., 1996; Dreeben and Pope, 1997a, 1998; Wacławczyk et al., 2004), and an application in channel flow with the current nonhybrid model is presented by Bakosi et al. (2007). The model to compute the joint PDF of velocity, scalar, and characteristic turbulent frequency is now complete.
2.2 Wall functions modelling of nearwall turbulence
When the wall region has sufficient resolution, and as defined by Equations (12) and (13) enable the model to adequately capture the nearwall effects on the higherorder turbulence statistics. Full resolution at the wall is strictly required in certain cases such as, for example, computations of heat transfer at walls embedded in a flow or detaching boundary layers with high adverse pressure gradients. In many other realistic simulations, however, full resolution of highReynoldsnumber boundary layers is not always possible and may not be necessary when the flow details close to walls are not important because the analysis focuses on the boundarylayer effects at farther distances. In these cases an alternative option is to model the nearwall turbulence using wall functions instead of the elliptic relaxation technique. Employing wall functions for noslip walls provides a tradeoff between the accuracy of fully resolved boundary layers and computational speed. Wall functions are widely applied in atmospheric simulations, where full wall resolution is usually prohibitively expensive even at the microscale or urban scale (Bacon et al., 2000; Lien et al., 2004). It is worth emphasizing that one of the main assumptions used in the development of wall functions is that the boundary layer remains attached, which is not always the case in simulations of complex flows. However, since wall functions are the only choice for realistic atmospheric simulations, they are still routinely employed with reasonable success.
To investigate the gain in performance and the effect on the results, we implemented the wall treatment for complex flow geometries that has been developed for the PDF method by Dreeben and Pope (1997b). In this case, the tensor is defined by the simplified Langevin model (SLM) (Haworth and Pope, 1986) and is simply a constant:
(19) 
with . In line with the purpose of wall functions, boundary conditions have to be imposed on particles that hit the wall so that their combined effect on the statistics at the first grid point from the wall will be consistent with the universal logarithmic wall function in equilibrium flows, i.e. in boundary layers with no significant adverse pressure gradients. The development of boundary conditions based on wall functions rely on the selfsimilarity of attached boundary layers close to walls. These conditions are applied usually at the first grid point from the wall based on the assumption of constant or linear stress distribution. This results in the wellknown selfsimilar logarithmic profile for the mean velocity. For the sake of completeness the conditions on the particles developed by Dreeben and Pope (1997b) are repeated here. The condition for the wallnormal component of the particle velocity reads
(20) 
where the subscripts and denote reflected and incident particle properties, respectively. The reflected streamwise particle velocity is given by
(21) 
where the coefficient is determined by imposing consistency with the logarithmic law at the distance of the first grid point from the wall, :
(22) 
where is a characteristic velocity scale of the turbulence intensity in the vicinity of , defined as
(23) 
with . , and are, respectively, the mean streamwise velocity, the wallnormal component of the Reynolds stress tensor and the turbulent kinetic energy at . In Equation (22) is the magnitude of the equilibrium value of the mean velocity at and is specified by the logarithmic law
(24) 
where is the Kármán constant and the surface roughness parameter for a smooth wall. The friction velocity is computed from local statistics as
(25) 
with
(26) 
In Equations (2025) the streamwise and wallnormal coordinate directions are defined according to the local tangential and normal coordinate directions of the particular wall element in question. In other words, if the wall is not aligned with the flow coordinate system then the vectors and , and the Reynolds stress tensor , need to be appropriately transformed into the wallelement coordinate system before being employed in the above equations.
The condition on the turbulent frequency is given by
(27) 
with
(28) 
which completes the description of the wall function approach.
In summary, the flow is represented by a large number of Lagrangian particles whose position , velocity , scalar concentration and characteristic turbulent frequency are governed by Equations (5), (6), (7) and (9), respectively. These equations are discretised and advanced in time by the explicit forward EulerMaruyama method (Kloeden and Platen, 1999). The mean pressure, required in Equation (6), is obtained via a pressure projection scheme (Bakosi et al., 2008). Full wall resolution is obtained through Equations (1218), while wall functions are applied through Equations (1928). The pressurePoisson and elliptic relaxation (14) equations are solved using an unstructured Eulerian grid with the finite element method. The grid is also used to track particles throughout the domain and to estimate Eulerian statistics using ensemble averaging. In practical simulations using PDF methods a few hundred particles per element is usually employed. Adequate stability can already be achieved using as little as 50–100 particles, however, 300–500 particles per elements are recommended to exploit the bin structure to compute (Bakosi et al., 2008) and to decrease the statistical error. The numerical algorithm and performance issues are detailed in Bakosi et al. (2008).
3 Modelling the street canyon
#  1  2  3  4  5  6  7  8  9  10  11  12 

0.5  1  1.5  2  2  2  2  2  2  2.5  2  2  
2  2  2  2  1.93  1.5  1.33  1  0.67  0.5  0.33  0.17  
#  13  14  15  16  17  18  19  20  21  22  23  24 
4  4  4  4  4  4  4  4  4  4.5  5  5.5  
0.17  0.33  0.5  0.67  1  1.33  1.5  1.93  2  2  2  2 
1.85  0.63  5  0.25  6  0.134  72  1.4  0.1  0.5  0.73  0.02  0.7 
Street canyons are often used to study flow and pollutant dispersion patterns in urban areas. A wealth of experimental data for this simplified urbanscale setting is available from windtunnel and LES data (Meroney et al., 1996; Pavageau and Schatzmann, 1999; Liu and Barth, 2002), making it a natural choice to validate the current, newly developed method. We will simulate the “urban roughness” case of Meroney et al. (1996), which is a model for a series of street canyons in the streamwise direction. The simulations are performed for statistically twodimensional flow geometry, with periodic inflow and outflow boundary conditions in the free stream above the buildings (i.e. the particles crossing the outflow boundary are reinjected at the inflow boundary). The Reynolds number based on the maximum free stream velocity and the building height was . This corresponds to based on the friction velocity and the free stream height, , if the free stream above the buildings is considered as the lower part of an approximate fully developed turbulent channel flow. The velocityconditioned scalar mean required in Equation (7) has been computed using the general method described by Bakosi et al. (2007) using a bin structure of . After the flow has reached a statistically stationary state, time averaging is used to collect velocity statistics and a continuous scalar is released from a street level line source at the centre of the canyon (corresponding to a point source in two dimensions). Particles travelling through the source of diameter are assigned a unit source strength. The scalar field is also time averaged after it has reached a stationary state.
The simulations with the full resolution model have been made with the constants given in Table 2, using 300 particles per element. The Eulerian mesh used for this simulation is displayed in Figure 1, which shows the considerable refinement along the building walls and tops necessary to resolve the boundary layers. In this case, the high anisotropy and inhomogeneity of the Reynolds stress tensor in the vicinity of walls are captured by the elliptic relaxation technique, using Equation (14).
The simulations using wall functions were performed on the Eulerian mesh displayed in Figure 2, also using 300 particles per element. The particleboundary conditions described in Sec. 2.2 were implemented for arbitrary geometry. Note that the first grid point where the boundary conditions based on wall functions are to be applied should not be closer to the wall than , where is the nondimensional distance from the wall in wall units, but should still be sufficiently close to the wall to lie in the inertial sublayer (Dreeben and Pope, 1997b). Accordingly, the grid in Figure 2 only contains the domain stripped from the wall region at .
Turbulence and scalar statistics are obtained entirely from the particles that represent both the flow itself and the scalar concentration field. The Eulerian meshes displayed in Figure 1 for the full resolution and in Figure 2 for the wall function cases are used to extract the statistics, to track the particles throughout the domain and to solve the Eulerian equations, namely Equation (14) and the meanpressurePoisson equation in the fully resolved case and only the latter in the wall function case.
In Figure 3, the mean velocity vector field and the isocontours of the turbulent kinetic energy are displayed for both fully resolved and wall function simulations. It is apparent that the full resolution captures even the smaller counterrotating eddies at the internal corners of the canyon, while the coarse grid resolution with wall functions only captures the overall flow pattern characteristic of the flow, such as the large steadily rotating eddy inside the canyon. The turbulent kinetic energy field is captured in a similar manner. Both methods reproduce the highest turbulence activity at the building height above the canyon, with a maximum at the windward building corner. The full resolution simulation shows a more detailed spatial distribution of energy, whereas the coarse resolution of the wallfunction simulation still allows tha capture of the overall pattern.
In Figure 4 the normalized turbulent intensities and are displayed for both simulation cases and compared with the LES results of Liu and Barth (2002). In the large eddy simulations the filtered momentum equations are solved by the Galerkin finite element method using brick threedimensional elements, while the residual stresses are modelled using the Smagorinsky closure.
The full resolution simulation shows very good agreement with LES. The contour plots of correctly display two local maxima, at the windward external and at the leeward internal corners. The contour plots of show distributed high values at the building level above the canyon, along the windward internal corner and wall, and at the street level downstream of the source. By contrast, the wallfunction contour plots are in general less detailed, failing to reproduce the internal maximum of , and showing a more uniform representation of .
Several windtunnel measurements have been carried out for a scalar released from a street level continuous line source at the centre of the canyon, providing concentration statistics above the buildings, at the walls, and inside the canyon (Meroney et al., 1996; Pavageau, 1996; Pavageau and Schatzmann, 1999). To examine the concentration values along the building walls and tops, we sampled the computed mean concentration field at the locations depicted in Figure 1 and listed in Table 1.
The excellent agreement of the results using both full resolution and wall functions with a number of experiments is shown in Figure 5. The concentration peak is precisely captured at the internal leeward corner and the model accurately reproduces the pattern of concentration along both walls including the higher values along the leeward wall.
In Figure 6, the first two statistical moments of the concentration inside the canyon are compared with experimental data and LES. The agreement with observations indicates that both the fluid dynamics and the micromixing components of the model provide a good representation of the real field.
Because the onepoint onetime joint PDF contains all higherorder statistics and correlations of the velocity and scalar fields resulting from a close, lowlevel interaction between the two fields, a great wealth of statistical information is available for atmospheric transport and dispersion calculations. As an example, in Figure 7 the cumulative distribution functions (CDF) of scalar concentration fluctuations
(29) 
are depicted after time averaging at selected locations of the domain for the full resolution case. No experimental data are available to assess the distributions in Figure 7, although the irregular shape of the CDF at , , which corresponds to a multimodal PDF, is likely to be an artifact of the micromixing model.
The performance gain obtained by applying wall functions as opposed to full resolution was about two orders of magnitude already at our moderate Reynolds number. The gain for higher Reynolds numbers is expected to increase faster than linearly.
4 Discussion
We have used an Eulerian unstructured grid, consisting of triangular element types, to estimate Eulerian statistics, to track particles throughout the domain, and to solve for inherently Eulerian quantities in conjunction with a PDF method. The boundary layers developing close to solid walls are fully captured with an elliptic relaxation technique, but can also be represented by wall functions, which use a coarser grid resolution and require significantly less particles, resulting in substantial savings in computational cost. We found that the onepoint statistics of the joint PDF of velocity and scalar are wellcaptured by the wall function approximation. In view of its affordable computational load and reasonable accuracy, this approximation appears to hold a realistic potential for application of the PDF method in atmospheric simulations, where the natural extension of the work is the implementation of the model in three spatial dimensions.
In hybrid PDF models developed for complex chemically reacting flows, numerical treatments for boundary conditions have been included for symmetric, inflow, outflow and freeslip walls employing the ghostcell approach common in finite volume methods (Rembold and Jenny, 2006). The representation of noslip boundaries adds a significant challenge to the above cases. This is partly due to the increased computational expense entailed by the higher Eulerian grid resolution that is required to fully resolve the boundary layers. In addition, there is an increased complexity in specifying the noslip particle conditions for both fully resolved and wall function representations. We presented an implementation of both approaches to treat noslip boundaries with unstructured grids in conjunction with the finite element method. This also obviates further complications with ghost cells.
In the case of full wall resolution we employed the Lagrangian equivalent of a modified isotropisation of production (IP) model as originally suggested by Dreeben and Pope (1998). The elliptic relaxation technique, however, allows for the application of any turbulence model developed for highReynoldsnumber turbulence (Durbin, 1993; Whizman et al., 1996). The standard test case for developing nearwall models is the fully developed turbulent channel flow. In this case, we explored the simpler Rotta (1951) model, which is the Eulerian equivalent of the simplified Langevin model (SLM) in the Lagrangian framework (Pope, 1994). This is simply achieved by eliminating the term involving the fourthorder tensor from the righthand side of Equation (14). While the SLM makes no attempt to represent the effect of rapid pressure (Pope, 2000) (in fact it is strictly correct only in decaying homogeneous turbulence), it is widely applied due to its is simplicity and robustness. Our experience showed a slight degradation of the computed velocity statistics (as compared to direct numerical simulation) using SLM for the case of channel flow. Since we experienced no significant increase in computational expense or decrease in numerical stability, we retained the original IP model.
Similarly, in the case of wall functions, several choices are available regarding the employed turbulence model. The methodology developed by Dreeben and Pope (1997b) uses the SLM, but it is general enough to include other more complex closures, such as the Haworth and Pope (1986, 1987) models (HP1 and HP2), the different variants of the IP models (IPMa, IPMb, LIPM) (Pope, 1994) or the Lagrangian version of the SSG model of Speziale et al. (1991). All these closures can be collected under the umbrella of the generalized Langevin model, by specifying its constants as described by Pope (1994). These models have been all developed for highReynoldsnumber turbulence and need to be modified in the vicinity of noslip walls. Including them in the wallfunction formulation is possible by specifying the reflected particle frequency at the wall as
(30) 
instead of Equation (28). This involves the additional computation of the statistics and at , which does not increase the computational cost significantly, but may result in a numerically less stable condition since the (originally constant) parameter that appears using the SLM has been changed to a variable that fluctuates during simulation. We implemented and tested all the above turbulence models using the wall function technique. Without any modification of the model constants we found the IPMa and SLM to be the most stable, providing very similar results. Thus we retained the original (and simplest) SLM along with Equation (28).
The most widely employed closure to model the small scale mixing of the passive scalar in the Lagrangian framework is the interaction by exchange with the mean (IEM) model of Villermaux and Devillon (1972) and Dopazo and O’Brien (1974). This simple and efficient model, however, fails to comply with several physical constraints and desirable properties of an ideal mixing model (Fox, 2003). The interaction by exchange with the conditional mean (IECM) model overcomes some of the difficulties inherent in the IEM model. In this study we justify the use of the IECM model by its being more physical and more accurate, but we acknowledge that it markedly increases the computational cost (Bakosi et al., 2008).
5 Acknowledgment
It is a pleasure to ackowledge the fruitful discussions during the course of the present work with Dr. Guido Cervone at George Mason University.
References
 Bacon et al. (2000) Bacon, D. P., Ahmad, N. N., Boybeyi, Z., Dunn, T. J., Hall, M. S., Lee, P. C. S., Sarma, R. A., Turner, M. D., III., K. T. W., Young, S. H., Zack, J. W., 2000. A dynamically adapting weather and dispersion model: The Operational Multiscale Environment Model with Grid Adaptivity (OMEGA). Mon. Weather Rev. 128, 2044–2076.
 Baik and Kim (1999) Baik, J.J., Kim, J.J., 1999. A numerical study of flow and pollutation dispersion characteristic in urban street canyons. J. Appl. Meteorol. 38, 1576–1589.

Bakosi et al. (2007)
Bakosi, J., Franzese, P., Boybeyi, Z., 2007. Probability density function
modeling of scalar mixing from concentrated sources in turbulent channel
flow. Phys. Fluids 19 (11), 115106.
URL http://link.aip.org/link/?PHF/19/115106/1  Bakosi et al. (2008) Bakosi, J., Franzese, P., Boybeyi, Z., 2008. A nonhybrid method for the PDF equations of turbulent flows on unstructured grids. J. Comput. Phys. 227 (11), 5896–5935.
 Cassiani et al. (2005a) Cassiani, M., Franzese, P., Giostra, U., 2005a. A PDF micromixing model of dispersion for atmospheric flow. Part I: development of the model, application to homogeneous turbulence and to neutral boundary layer. Atmos. Environ. 39 (8), 1457–1469.
 Cassiani et al. (2005b) Cassiani, M., Franzese, P., Giostra, U., 2005b. A PDF micromixing model of dispersion for atmospheric flow. Part II: application to convective boundary layer. Atmos. Environ. 39 (8), 1471–1479.
 Cassiani and Giostra (2002) Cassiani, M., Giostra, U., 2002. A simple and fast model to compute concentration moments in a convective boundary layer. Atmos. Environ. 36 (30), 4717–4724.
 Cassiani et al. (2007) Cassiani, M., Radicchi, A., Albertson, J. D., 2007. Modelling of concentration fluctuations in canopy turbulence. BoundaryLayer Meteorol. 122 (3), 655–681.
 Chatwin and Sullivan (1993) Chatwin, P. C., Sullivan, P. J., 1993. The structure and magnitude of concentration fluctuations. BoundaryLayer Meteorol. 62, 269–280.
 Dopazo (1994) Dopazo, C., 1994. Recent developments in pdf methods. In: Libby, P. A. (Ed.), Turbulent reactive flows. Academic, New York, pp. 375–474.
 Dopazo and O’Brien (1974) Dopazo, C., O’Brien, E. E., 1974. An approach to the autoignition of a turbulent mixture. Acta Astronaut. 1, 1239–1266.

Dreeben and Pope (1997a)
Dreeben, T. D., Pope, S. B., 1997a. Probability density function
and Reynoldsstress modeling of nearwall turbulent flows. Phys. Fluids
9 (1), 154–163.
URL http://link.aip.org/link/?PHF/9/154/1 
Dreeben and Pope (1997b)
Dreeben, T. D., Pope, S. B., 1997b. Wallfunction treatment in pdf
methods for turbulent flows. Phys. Fluids 9 (9), 2692–2703.
URL http://link.aip.org/link/?PHF/9/2692/1  Dreeben and Pope (1998) Dreeben, T. D., Pope, S. B., 1998. Probability density function/Monte Carlo simulation of nearwall turbulent flows. J. Fluid Mech. 357, 141–166.
 Durbin (1993) Durbin, P. A., 1993. A Reynolds stress model for nearwall turbulence. J. Fluid Mech. 249, 465–498.

Fox (1996)
Fox, R. O., 1996. On velocityconditioned scalar mixing in homogeneous
turbulence. Phys. Fluids 8 (10), 2678–2691.
URL http://link.aip.org/link/?PHF/8/2678/1  Fox (2003) Fox, R. O., 2003. Computational models for turbulent reacting flows. Cambridge University Press.
 Franzese (2003) Franzese, P., 2003. Lagrangian stochastic modeling of a fluctuating plume in the convective boundary layer. Atmos. Environ. 37, 1691–1701.
 Geuzaine and Remacle (2009) Geuzaine, C., Remacle, J.F., 2009. Gmsh: a threedimensional finite element mesh generator with builtin pre and postprocessing facilities (accepted for publication). Int. J. Numer. Meth. Eng.
 Grigoryev et al. (2002) Grigoryev, Y. N., Vshivkov, V. A., Fedoruk, M. P., 2002. Numerical “particleincell” methods: theory and applications. Utrecht, Boston.

Haworth and Pope (1986)
Haworth, D. C., Pope, S. B., 1986. A generalized Langevin model for turbulent
flows. Phys. Fluids 29 (2), 387–405.
URL http://link.aip.org/link/?PFL/29/387/1 
Haworth and Pope (1987)
Haworth, D. C., Pope, S. B., 1987. A pdf modeling study of selfsimilar
turbulent free shear flows. Phys. Fluids 30 (4), 1026–1044.
URL http://link.aip.org/link/?PFL/30/1026/1  Hoydysh et al. (1974) Hoydysh, W. G., Griffiths, R. A., Ogawa, Y., 1974. A scale model study of the dispersion of pollution in street canyons. In: 67th Annual Meeting of the Air Pollution Control Association, Denver, Colorado, APCA Paper No. 74157.
 Huang et al. (2000) Huang, H., Akutsu, Y., Arai, M., Tamura, M., 2000. A twodimensional air quality model in an urban street canyon: Elevation and sensitivity analysis. Atmos. Environ. 34, 689–698.
 Jenny et al. (2001) Jenny, P., Pope, S. B., Muradoglu, M., Caughey, D. A., 2001. A hybrid algorithm for the joint PDF equation of turbulent reactive flows. J. Comput. Phys. 166, 218–252.
 Johnson and Hunter (1995) Johnson, G. T., Hunter, L. J., 1995. A numerical study of dispersion of passive scalar in city canyons. BoundaryLayer Meteorol. 75, 235–262.
 Kloeden and Platen (1999) Kloeden, P. E., Platen, E., 1999. Numerical Solution of Stochastic Differential Equations. Springer, Berlin.
 Kristensen (1994) Kristensen, L., 1994. Recurrence of extreme concnetrations. In: S.V. Grining, M. M. Millán (Eds.), Air Pollution Modelling and Its Application. X. Plenum Press, New York.
 Lee and Park (1994) Lee, I. Y., Park, H. M., 1994. Parameterization of the pollutant transport and dispersion in urban street canyons. Atmos. Environ. 28, 2343–2349.
 Lien et al. (2004) Lien, F. S., Yee, E., Cheng, Y., 2004. Simulation of mean flow and turbulence over a 2d building array using highresolution CFD and a distributed drag force approach. J. Wind Eng. Ind. Aerod. 92 (2), 117–158.
 Liu and Barth (2002) Liu, C.H., Barth, M. C., 2002. Largeeddy simulation of flow and scalar transport in a modeled street canyon. J. Appl. Meteorol. 41 (6), 660–673.
 Luhar et al. (2000) Luhar, A. K., Hibberd, M. F., Borgas, M. S., 2000. A skewed meanderingplume model for concentration statistics in the convective boundary layer. Atmos. Environ. 34, 3599–3616.

Lundgren (1969)
Lundgren, T. S., 1969. Model equation for nonhomogeneous turbulence. Phys.
Fluids 12 (3), 485–497.
URL http://link.aip.org/link/?PFL/12/485/1  Meroney et al. (1996) Meroney, R. N., Pavageau, M., Rafailidis, S., Schatzmann, M., 1996. Study of line source characteristics for 2d physical modelling of pollutant dispersion in street canyons. J. Wind Eng. Ind. Aerod. 62 (1), 37–56.
 Mobus et al. (2001) Mobus, H., Gerlinger, P., Bruggemann, D., 2001. Comparison of eulerian and lagrangian monte carlo pdf methods for turbulent diffusion flames. Combust. Flame 124, 519–534.
 Muradoglu et al. (1999) Muradoglu, M., Jenny, P., Pope, S. B., Caughey, D. A., 1999. A consistent hybrid finitevolume/particle method for the PDF equations of turbulent reactive flows. J. Comput. Phys. 154, 342–371.
 Muradoglu et al. (2001) Muradoglu, M., Pope, S. B., Caughey, D. A., 2001. The hybrid method for the PDF equations of turbulent reactive flows: consistency conditions and correction algorithms. J. Comput. Phys. 172, 841–878.
 Pavageau (1996) Pavageau, M., 1996. Concentration fluctuations in urban street canyons – groundwork for future studies. Tech. rep., Meteorological Institute of the University of Hamburg.
 Pavageau and Schatzmann (1999) Pavageau, M., Schatzmann, M., 1999. Wind tunnel measurements of concentration fluctuations in an urban street canyon. Atmos. Environ. 33, 3961–3971.
 Pope (1985) Pope, S. B., 1985. PDF methods for turbulent reactive flows. Prog. Energ. Combust. 11, 119–192.

Pope (1994)
Pope, S. B., 1994. On the relationship between stochastic Lagrangian models
of turbulence and secondmoment closures. Phys. Fluids 6 (2), 973–985.
URL http://link.aip.org/link/?PHF/6/973/1  Pope (1998) Pope, S. B., 1998. The vanishing effect of molecular diffusivity on turbulent dispersion: implications for turbulent mixing and the scalar flux. J. Fluid Mech. 359, 299–312.
 Pope (2000) Pope, S. B., 2000. Turbulent flows. Cambridge University Press, Cambridge.
 Rafailids and Schatzmann (1995) Rafailids, S., Schatzmann, M., 1995. Physical modelling of car exhaust disperision in urban street canyons. In: Proc. 21st Int. Meeting on Air Pollution Modelling and Its Applications, Baltimore, Nov. 610.
 Rembold and Jenny (2006) Rembold, B., Jenny, P., 2006. A multiblock joint pdf finitevolume hybrid algorithm for the computation of turbulent flows in complex geometries. J. Comput. Phys. 220, 59–87.
 Rotta (1951) Rotta, J. C., 1951. Statistiche theorie nichthomogener turbulenz. Z. Phys. 129, 547.
 Sawford (2004) Sawford, B. L., 2004. Micromixing modelling of scalar fluctuations for plumes in homogeneous turbulence. Flow Turbul. Combust. 72, 133–160.

Sawford (2006)
Sawford, B. L., 2006. Lagrangian modeling of scalar statistics in a double
scalar mixing layer. Phys. Fluids 18 (8), 085108.
URL http://link.aip.org/link/?PHF/18/085108/1  Soulard and Sabel’nikov (2006) Soulard, O., Sabel’nikov, V. A., 2006. Eulerian Monte Carlo method for the joint velocity and massfraction probability function in turbulent reactive gas flows. Combust. Explo. Shock. 42, 753–762.
 Speziale et al. (1991) Speziale, C. G., Sarkar, S., Gatski, T. B., 1991. Modelling the pressurestrain correlation of turbulence: an invariant dynamical systems approach. J. Fluid Mech. 227, 245–272.
 Valiño (1998) Valiño, L., 1998. A field Monte Carlo formulation for calculating the probability density function of a single scalar in a turbulent flow. Flow Turbul. Combust. 60, 157–172.
 van Kampen (2004) van Kampen, N. G., 2004. Stochastic processes in physics and chemistry, 2nd Edition. North Holland, Elsevier B.V., Amsterdam, The Netherlands.

van Slooten et al. (1998)
van Slooten, P. R., Jayesh, Pope, S. B., 1998. Advances in PDF modeling for
inhomogeneous turbulent flows. Phys. Fluids 10 (1), 246–265.
URL http://link.aip.org/link/?PHF/10/246/1  Villermaux and Devillon (1972) Villermaux, J., Devillon, J. C., 1972. Représentation de la coalescence et de la redispersion des domaines de ségrégation dans un fluide par un modèle d’interaction phénoménologique. In: Proceedings of the 2nd Intl. Symp. on Chemical Reaction Engineering. Elsevier, New York, pp. 1–13.

Wacławczyk et al. (2004)
Wacławczyk, M., Pozorski, J., Minier, J.P., 2004. Probability density
function computation of turbulent flows with a new nearwall model. Phys.
Fluids 16 (5), 1410–1422.
URL http://link.aip.org/link/?PHF/16/1410/1  Wedding et al. (1977) Wedding, J. B., Lambert, D. J., Cermak, J. E., 1977. A wind tunnel study of gaseous pollutants in city street canyons. Air Pollut. Control Assoc. J. 27, 557–566.
 Whizman et al. (1996) Whizman, V., Laurence, D., Kanniche, M., Durbin, P. A., Demuren, A., 1996. Modeling nearwall effects in secondmoment closures by elliptic relaxation. Int. J. Heat Fluid Fl. 17, 255–266.
 Wilson (1995) Wilson, D. J., 1995. Concentration fluctuations and averaging time in vapor clouds. Center for Chemical Process Safety. American Institute of Chemical Engineers, New York.
 Yee et al. (1994) Yee, E., Chan, R., Kosteniuk, P. R., Chandler, G. M., Biltoft, C. A., Bowers, J. F., 1994. Incorporation of internal fluctuations in a meandering plume model of concentration fluctuations. BoundaryLayer Meteorol. 67 (1–2), 11–39.
 Yee and Wilson (2000) Yee, E., Wilson, D. J., 2000. A comparison of the detailed structure in dispersin tracer plumes measured in gridgenerated turbulence with a meandering plume model incorporating internal fluctuations. BoundaryLayer Meteorol. 94 (2), 253–296.