Articles | Volume 16, issue 20
Model description paper
17 Oct 2023
Model description paper |  | 17 Oct 2023

QES-Plume v1.0: a Lagrangian dispersion model

Fabien Margairaz, Balwinder Singh, Jeremy A. Gibbs, Loren Atwood, Eric R. Pardyjak, and Rob Stoll

Low-cost simulations providing accurate predictions of transport of airborne material in urban areas, vegetative canopies, and complex terrain are demanding because of the small-scale heterogeneity of the features influencing the mean flow and turbulence fields. Common models used to predict turbulent transport of passive scalars are based on the Lagrangian stochastic dispersion model. The Quick Environmental Simulation (QES) tool is a low-computational-cost framework developed to provide high-resolution wind and concentration fields in a variety of complex atmospheric-boundary-layer environments. Part of the framework, QES-Plume, is a Lagrangian dispersion code that uses a time-implicit integration scheme to solve the generalized Langevin equations which require mean flow and turbulence fields. Here, QES-Plume is driven by QES-Winds, a 3D fast-response model that computes mass-consistent wind fields around buildings, vegetation, and hills using empirical parameterizations, and QES-Turb, a local-mixing-length turbulence model. In this paper, the particle dispersion model is presented and validated against analytical solutions to examine QES-Plume’s performance under idealized conditions. In particular, QES-Plume is evaluated against a classical Gaussian plume model for an elevated continuous point-source release in uniform flow, the Lagrangian scaling of dispersion in isotropic turbulence, and a non-Gaussian plume model for an elevated continuous point-source release in a power-law boundary-layer flow. In these cases, QES-Plume yields a maximum relative error below 6 % when compared with analytical solutions. In addition, the model is tested against wind-tunnel data for a uniform array of cubical buildings. QES-Plume exhibits good agreement with the experiment with 99 % of matched zeros and 59 % of the predicted concentrations falling within a factor of 2 of the experimental concentrations. Furthermore, results also emphasize the importance of using high-quality turbulence models for particle dispersion in complex environments. Finally, QES-Plume demonstrates excellent computational performance.

1 Introduction

Rapid growth of urban populations around the world impacts all sectors of human activity, including industry and transportation. Additionally, growth is also increasing pressure on agricultural systems to boost yields. These trends raise concerns about the deterioration of the environment, a decline in quality of life, or worsening air quality (Britter and Hanna2003). Along with these long-term problems, acute risks including the potential accidental or deliberate release of a chemical or biological agent pose a major threat in densely populated urban areas (Gowardhan et al.2021).

In response to these and other issues, a number of fast-response transport and dispersion models for urban areas and complex terrain have been developed. Fast-response models are characterized by their ability to keep computational costs low while providing realistic representations of the effects of buildings, canopies, and terrain on velocity distributions and the dispersion of scalars (Pardyjak et al.2008; Singh et al.2008; Gowardhan et al.2011). Different classes of models exist for scalar dispersion. First, reduced-order models, like Gaussian dispersion models, use algebraic descriptions of plumes to calculate concentrations (Hanna et al.2003; Philips et al.2013; Prussin et al.2015; Miller et al.2018). Examples of operational reduced-order models include AERMOD (Cimorelli et al.2005) or SIRANE (Soulhac et al.2011). However, the performance of Gaussian models in complex environments is limited as this class of model struggles to capture aspects of the plume such as asymmetry, plume turning, or other critical features especially in urban settings (Hertwig et al.2018; Pirhalla et al.2021). Another class of models follows an Eulerian approach that solves the prognostic equation for concentration, such as the large-eddy-simulation suite PALM (Maronga et al.2020), the hybrid Eulerian–Lagrangian Code Saturne (Archambeau et al.2004; Bahlali et al.2018), and the CMAQ Modeling System (Byun and Schere2006; US EPA Office Of Research And Development2020). However, the sophistication and computational cost result in these models being unsuitable for rapid deployment at urban or local scales. The use of these models is often more appropriate for research or at large, continental, or global scales. A third class, Lagrangian stochastic dispersion models (LSDMs), describe the movement of each particle in turbulent flows using a random-walk approach (Pope1987; Thomson1987; Rodean1996).

Examples of operational LSDMs are QUIC-PLUME (Williams et al.2002), MicroSpray (Tinarelli et al.2013), and GRAMM/GRAL (Oettl2015). LSDMs have the benefit of being easily parallelized thanks to their mathematical formulation (Singh et al.2011). In addition to the dispersion model, accurate estimations of concentrations of airborne substances require mean-wind and turbulence models such as QUIC-URB (Pardyjak and Brown2003), SIRANE (Soulhac et al.2011) or street-network models (Hertwig et al.2018). The importance of these models was emphasized by Carissimo et al. (2021), who evaluated how different numerical modeling approaches resolved the concentration in the wake of buildings.

Another area of study where fast-dispersion models can have a big impact is aerobiology, the study of the aerial dispersion of biological particles in the atmosphere (Legg1983; Aylor2017). In particular, the accurate prediction of temporal and spatial dispersion of pathogens can provide crucial information related to the development of plant disease epidemics or for fungicide resistance management (Thiessen et al.2016). An overview of the use of aerial sampling for the detection of epidemic development can be found in Mahaffee et al. (2023). Similarly, dispersion models can be used to study long-distance transport of pollen and spores (e.g., Aylor2003). In recent years, there has been an increase in interest from growers to adopt techniques from aerobiology to improve yields (Thiessen et al.2016).

Numerical integration of the equations used in LSDMs is challenging due to the stiffness of the model's mathematical formulation (Yee and Wilson2007). An equation is considered stiff when some numerical methods are unstable, unless extremely small time or space steps are used. In particular, for LSDMs a particle may travel large distances over small time steps because of the existence of instabilities in the numerical methods. Simplified versions of the equations are not as unstable because of the reduced number of terms compared to the generalized model. However, simplified models are still numerically unstable. Furthermore, early numerical integration methods were known for violating the well-mixed condition, where uniformly distributed particles become unmixed even under homogeneous conditions because of numerical instabilities (Thomson1984). This phenomenon has been attributed to “rogue” trajectories, where particles accumulate energy and develop arbitrarily large velocity fluctuations (Yee and Wilson2007; Postma et al.2012; Wilson2012; Postma2015).

Several methods have been presented to integrate LSDMs. Yee and Wilson (2007) proposed a fractional step methodology to partially circumvent the stiffness in the generalized model. Ramli and Esler (2016) outlined a rigorous methodology to evaluate numerical schemes for LSDMs where a series of one-dimensional test problems were introduced based on the Fokker–Planck equation. They conclude that if long time steps have to be used – for long-distance transport, for example – it can be beneficial to use a random displacement model approximation rather than classical integration schemes. Bailey (2017) investigated the possibility of using a time-implicit scheme to eliminate rogue trajectories. He showed that numerical instabilities of the temporal integration scheme lead to nonphysical trajectories and that a lagged implicit scheme is unconditionally stable for the generalized model. The conclusion of Bailey (2017) motivated the use of his methodology in the present work.

To address the issues presented above, the Quick Environmental Simulation (QES) framework was developed to provide high-resolution wind and concentration fields in complex urban and agricultural environments. The framework is composed of QES-Winds, a 3D fast-response model that computes mass-consistent wind fields around buildings and vegetation using empirical parameterizations (Bozorgmehr et al.2021; Margairaz et al.2022b); QES-Turb, a turbulence model based on an eddy-viscosity parameterization with a local mixing length (Pope2000); and QES-Plume, a particle dispersion model. The objective of this work is to describe and validate QES-Plume.

In the following, we first describe the mathematical formulation of the LSDM in Sect. 2. Next, Sect. 3 introduces the QES framework, and the results of the model validation are presented in Sect. 4. The importance of the wind flow and turbulence models is discussed in Sect. 5 and, finally, conclusions are presented in Sect. 6.

2 Lagrangian stochastic dispersion model

The motion of passive tracers in turbulent flow can be described by a random-walk model (Thomson1987). The time evolution of a fluid particle's position is given by

(1) d x p , i d t = U i + u i ,

where xp,i is the position of the particle p in a Cartesian coordinate system i{1,2,3}, t is time, Ui is the particle's mean velocity, and ui is its velocity fluctuation. Here, a “particle” is a statistical representation of fluid element containing many molecules (Pisso et al.2019) and not an actual aerosol particle. The fluctuation is modeled by a random-walk process formalized in the Langevin equations (Langevin1908). This stochastic differential equation is given by

(2) d u i = - a u i d t + b i j d W j ,

where Einstein notation is used to imply summation over repeating indices. In Eq. (2), a is a damping coefficient associated with viscous drag (Rodean1996) and bij is a scaling tensor for the three independent random variates dWj representing Brownian (or Wiener) processes (Thomson1987). The random terms dWj follows a Gaussian distribution with zero mean and variance dt (Yee and Wilson2007).

For stationary, homogeneous, and isotropic turbulence, a simplified model is obtained by writing a and bij as a function of the Lagrangian velocity timescale τL. In this case, the simplified Langevin equations (SLEs) are given by

(3) d u i = - 1 τ L u i d t + 2 σ 2 τ L 1 / 2 d W i .

Following Rodean (1996), the timescale τL can be parameterized as

(4) τ L = 2 σ 2 C 0 ε ,

where σ2 is the velocity variance, C0 is a universal constant, and ε is the mean dissipation rate of turbulence kinetic energy (TKE). For C0, Rodean (1991) reported values ranging from 1.6 to 10, with Du (1997) proposing a value of C0=3.0±0.5, and for a constant-stress region, Rodean (1991) derived a semi-analytical value of C0=5.7.

The model in Eq. (3) can be expanded for Gaussian, inhomogeneous, and anisotropic turbulence. Following Thomson (1987), the 3D generalized Langevin equations (GLEs) are given by

(5) d u i = - C 0 ε 2 τ i k - 1 u k d t + τ j - 1 2 d τ i d t u j d t + 1 2 τ i x d t + C 0 ε 1 / 2 d W i ,

where τij is the Reynolds stress tensor and τij-1 is its inverse. One of the consequences of the hypothesis of Gaussian turbulence is that the Reynolds stress tensor must be positive definite, otherwise the probability distribution function of the velocity fluctuations becomes ill-defined (Bailey2017). A tensor is positive definite if and only if all of its eigenvalues are positive. Since the GLEs require the existence of an inverse of τij, the modeled stress tensor has to be positive definite. Therefore, all eigenvalues and principal invariants of the stress tensor have to be positive; i.e., it must be realizable (Vachat1977).

The GLEs are considered stiff because of the presence of the wide range of timescales. In particular, a particle may travel a significant distance over a small time step because of the existence of instabilities in the numerical solution. The SLEs drastically reduce the number of terms in the GLEs and are considered less unstable, however, numerical instabilities still exist for the SLEs (Yee and Wilson2007).

On the ground and at building surfaces, reflecting boundary conditions are used. At the domain top, an outlet condition is used, indicative of particles traveling past the top of the domain. Yet, the top boundary conditions can be more complex if the top of the atmospheric boundary layer is considered (see Thomson et al.1997).

3 Method

To calculate the trajectory of fluid particles, the QES framework is composed of three main modules: (1) QES-Winds, a mass-consistent wind model; (2) QES-Turb, a turbulence model; and (3) QES-Plume, an LSDM. Figure 1 shows the workflow of the QES framework. First, the steady-state wind field is computed. The resulting velocity vector (Ui) is passed to the turbulence model, where the Reynolds stress tensor (τij) and the mean dissipation rate of TKE (ε) are calculated. Finally, all variables are passed to QES-Plume. QES's modules are coded in C++ and NVIDIA's CUDA application programming interface (Margairaz et al.2022a). The input and output take advantage of the netCDF (Network Common Data Form) library developed by UCAR/Unidata (Rew et al.1989). In addition, QES-Plume can also run as a stand-alone model with pre-computed velocity and turbulence fields. In this mode, QES-Plume is initialized directly from netCDF files. This mode is mainly used for verification and validation purposes.

Figure 1Simulation workflow of the QES framework.


3.1 QES-Winds: mass-consistent wind model

The wind field is obtained using QES-Winds. This model is based on the framework introduced by Sherman (1978) and Röckle (1990) to simulate a divergence-free steady-state 3D wind field around buildings. The framework uses a combination of a mass-consistent diagnostic wind model and various empirical parameterizations to describe the dynamics of the flow around buildings or through vegetated canopies. The mass-consistent solver employs a variational-analysis technique introduced by Sasaki (1970) where the Euler–Lagrange equation is used to minimize the error between the guess and final flow field under the divergence-free constraint. The minimization process is achieved through the solution of a Poisson equation with a successive over-relaxation (SOR) solver parallelized on GPUs using a red–black successive over-relaxation method (Bozorgmehr et al.2021). The impact of momentum transport on flow around buildings is parameterized using models inspired by the QUIC-URB model (Pardyjak and Brown2003), namely, the building rooftop recirculation or velocity attenuation, the upstream recirculation zone, and the downwind recirculation zone and velocity deficit wake (Röckle1990; Singh et al.2008; Gowardhan et al.2010; Brown et al.2013). In addition, a street canyon parameterization has been implemented based on Singh et al. (2008). Flow through vegetated canopies is parameterized using bulk canopy attenuation (Pardyjak et al.2008) and isolated-tree (Margairaz et al.2022b) and row-organized canopy models (Ulmer et al.2023).

3.2 QES-Turb: turbulence model

To accurately represent particle motion in turbulent flows, the LSDM needs the stress tensor and the dissipation rate of TKE. The QES framework calculates turbulence variables using a local-mixing model based on Prandtl’s mixing-length and the Boussinesq eddy-viscosity hypotheses (Pope2000). The mixing-length model was selected following a fast-response philosophy, favoring low-cost algebraic models rather than models based on transport equations. From the turbulent-viscosity hypothesis, the Reynolds stresses (τij=uiuj, where represents the time-averaging operator) are given by

(6) τ i j = 2 3 k δ i j - 2 ν T S i j ,

where k=12(uiui) is the TKE, δij is the Kronecker delta tensor, νT is the eddy viscosity, and Sij is the mean strain-rate tensor given by

(7) S i j = 1 2 U i x j + U j x i .

The eddy viscosity is modeled as a product of a characteristic length scale and velocity scale. For boundary-layer flows, the classical approach is to calculate the eddy viscosity as

(8) ν T = m 2 U z ,

with the mixing length specified as m=κz, where κ≈0.4 is the von Kármán constant and z is the height above the ground. Smagorinsky (1963) proposed a generalization for the eddy viscosity calculated using the mean strain-rate tensor. This formulation is given by

(9) ν T = m 2 2 S i j S i j 1 / 2 .

For complex flows, such as flows around buildings, the mixing length is computed as m=κdW, where dW is the distance to the closest wall.

From the eddy-viscosity model, the TKE can be defined as

(10) k = ν T C T m 2 ,

where CT is a model constant. In QES-Turb, this constant is set to CT=0.55, to obtain the correct behavior for the log-law region (Pope2000, Chap. 10). Following the same philosophy, the mean dissipation rate ε scales as a velocity scale cubed and then divided by a length scale. Hence, ε is modeled as

(11) ε = C D 3 k 3 / 2 m ,

where CD is a model constant. To be consistent with the definition of the TKE, this constant is related to CT such that CD=CT (i.e., CD=0.55; Pope2000, Chap. 10).

The strain-rate tensor from Eq. (7) is calculated from the local gradients of the steady-state velocity fields simulated by QES-Winds. The velocity gradients are computed using central finite differences away from the walls. At the walls, the derivatives are computed with first-order forward finite differences.

For boundary-layer flows, there is ample evidence suggesting that the normal stresses are anisotropic and that for flow aligned with the x1 direction τ11>τ22>τ33. However, the stress formulations in Eq. (6) do not guarantee that the modeled stress emulates these observations. Hence, to ensure that the stress reflects the anisotropic condition, the following formulation is used:


where the constants C{u,v,w} are defined as the velocity component standard deviations normalized by the friction velocity (i.e., Cu=σu/u*, Cv=σv/u*, and Cw=σw/u*). Table 1 summarizes the value of these constants for different types of flows found in the literature for neutral stability. The current version of the turbulence model does not account for atmospheric stability, as this study is focused on the urban canopy sublayer, where the atmosphere is mostly neutral (Ramamurthy et al.2007). Future versions should include corrections to these coefficients in the surface layer following Monin–Obukhov similarity theory (Stiperski and Calaf2023).

Roth (2000)Roth (2000)Britter and Hanna (2003)Brunet (2020)Finnigan (2000)Brunet (2020)Finnigan (2000)Brunet (2020)Castelli et al. (2001)

Table 1Summary of C{u,v,w} values from the literature for different types of flows (* value not reported).

Download Print Version | Download XLSX

Because the turbulence model is a local-mixing model, it relies heavily on the magnitude of the local velocity gradients to estimate the stress tensor. This is problematic in regions where the velocity gradients are small and leads to the model predicting negligible turbulence, for example at the core of a street canyon. While multiple methods exist to address these issues (e.g., Williams et al.2002), a non-local turbulence mixing coefficient Cnlm is added to the diagonal elements of the stress tensor everywhere in the domain to enhance the mixing and account for processes such as advection of turbulence into regions with small gradients. The result of this addition is an increase in the magnitude of the turbulence while maintaining the influence of local velocity gradients in the stress tensor. The non-local-mixing coefficient, Cnlm, may not be valid for general applicability as this may yield unreasonable results unless validated extensively with more experimental data. For the array of cubical buildings presented in this paper, this coefficient is set to Cnlm≈0.3 m2 s−2 based on background mixing from the wind-tunnel data (see Appendix A2).

3.3 QES-Plume: Lagrangian stochastic dispersion model

3.3.1 Numerical implementation

A common approach is the use of an explicit scheme to integrate the GLEs. The scheme is conditionally stable with a condition on the time step related to the velocity variance and mean dissipation rate. However, due to the stiffness of the equation, extremely small time steps are required to maintain stability and numerical errors can inject more energy than the viscosity can dissipate. As a consequence, the particle velocity can become arbitrarily large, leading to rogue trajectories (Yee and Wilson2007). Bailey (2017) proposed using an implicit numerical scheme to eliminate the numerical instability. However, a fully implicit scheme for the GLE is challenging because of forcing terms. In order to simplify the problem and avoid coupling between the position and Langevin equations, the forcing terms are evaluated at the current particle position. This approach, called “lagged” forcing, is used to avoid a costly iterative scheme. The implicit scheme for the velocity fluctuations from t(n) to t(n+1)=t(n)+Δt is given by

(15) u i ( n + 1 ) = u i ( n ) + [ - C 0 ε 2 τ i k - 1 ( n ) u k ( n + 1 ) + τ j - 1 2 Δ τ i Δ t ( n ) u j ( n + 1 ) + 1 2 τ i x ( n ) ] Δ t + C 0 ε ( n ) 1 / 2 Δ W i ,

where ui(n) is the fluctuation at time t(n) and ui(n+1) is the fluctuation at time t(n+1). Note that all the forcing terms are lagged (i.e., evaluated at time t(n)). In particular, the time derivative of the stress tensor is approximated by

(16) d τ i j d t Δ τ i j Δ t ( n ) = τ i j ( n ) - τ i ( n - 1 ) Δ t .

This scheme is unconditionally stable but the system of equations in Eq. (15) requires the solution of a 3×3 matrix. Similar to the observation made by Bailey (2017), the authors have never observed the matrix to become singular, provided that the forcing terms are well defined.

Finally, a forward Euler scheme is used to update the particle position. The position xi(n+1) is given by

(17) x i ( n + 1 ) = x i ( n ) + U i ( n ) + u i ( n + 1 ) Δ t ,

where Ui(n) is the mean velocity at the position xi(n).

Figure 2Workflow of the QES-Plume model using the implicit scheme to solve the 3D GLE. The particle position is advanced from its position xi(n) at t(n) to its new position xi(n+1) at t(n+1)=t(n)+Δt. The algorithm repeats the advection loop until the time reaches the final time Tfinal.


The workflow of the QES-Plume model is presented in Fig. 2. At each time step, all sources emit particles based on each source's parameters. Typically, the particle's initial fluctuations are obtained by interpolating the velocity variances at the source location. Then, the motion of all particles is computed in the advection routine. First, the value of the velocity, stress tensor, stress tensor divergence, and mean dissipation rate are interpolated at the particle position using a trilinear interpolation method. To ensure that the model is well defined, the stress tensor must be made realizable if it is not, which can be either because of the turbulence model or interpolation. The stress tensor is realizable if all three of its principal invariants are strictly larger than a tolerance of 10−5 (Bailey2017). If this condition is not satisfied, the diagonal elements of the stress tensor are incremented by 5 % of the mean TKE. This procedure is repeated until all three principal invariants are larger than the selected tolerance. Enforcing this condition guarantees that the stress tensor can be inverted. Next, the velocity fluctuations are computed by solving Eq. (15) following the procedure in Fig. 2. Here, both the stress tensor τij and the matrix Aij are explicitly inverted using the adjugate method (Horn and Johnson2012). Finally, the position of the particle is updated using Eq. (17) and the boundary conditions are applied, such as a reflecting boundary condition at the ground or on building faces (see Sect. 3.3.3) or an outlet condition at the sides and top of the domain (where particles are deleted once they travel outside of the domain). Once all particle positions have been updated, the relevant statistics are computed; in particular, the airborne concentration, and the time is incremented.

The new implementation of the LSDM has been tested following the procedure proposed by Bailey (2017). These tests checked that the well-mixed condition is respected for a wide range of forcing conditions and time steps. To verify that unmixing did not occur, particles were uniformly distributed within the domain, the LSDM was run with different flow conditions and time steps, and the final distribution of particle positions was checked. Three forcing conditions have been considered: (i) synthetic data with zero flow and sinusoidal vertical stress, (ii) channel flow data from direct numerical simulation (Kim et al.1987), and (iii) data from large-eddy simulation of the atmospheric boundary layer from Stoll and Porté-Agel (2006) (see Bailey2017, for more details about the tests).

3.3.2 Dynamic time step

The method presented eliminates the possibility of the calculation of a rogue trajectory during particle advection. However, some physical processes, such as reflection or deposition, require the particle to travel only one Eulerian grid cell (i.e., QES-Winds velocity grid) at a time. To control the total distance traveled, a Courant-number-based algorithm is used. The time step is reduced as the particle moves close to a wall. The new time step is calculated following the procedure presented in Algorithm 1. The user-specified time step Δt has to be progressively reduced as the particle moves closer to the wall using a user-specified Courant number CN. The progressive reduction in the time step is required to avoid carrying large fluctuations into a smaller time increment.

Algorithm 1Courant-Number-based time step reduction. Δt is the user-specified time step. CN is the user-specified Courant Number (CN<1).

d distance from the particle to the closest wall;
if d>6UΔt then
if d>3UΔt then
end if
end if

This procedure is used to divide the user-defined time step into smaller time increments to ensure that the particle travels small distances. Hence, other algorithms requiring that the particle only travels one grid cell can be executed correctly.

3.3.3 Reflecting boundary condition

Perfect rebound is used as wall boundary conditions for the particles (Brzozowska2013; Bahlali et al.2020). Figure 3 represents the trajectory of a particle with a double reflection in a corner. The wall-reflection method is presented in Algorithm 2. This method is called every time a particle crosses the fluid–solid interface.

Algorithm 2Wall reflection.

while Particle in solid do
Find closest face along trajectory
Find intersection with face p
d is the portion of the trajectory beyond the intersection with the face
end while

Figure 3Example of double reflection in a corner from particle position x(n) to x(n+1). Dotted lines represent the part of the trajectory in the solid, represented by d. Dashed arrow represents the reflection r=2(dn^)n^, where n^ is the wall normal. The particle position after first reflection is x=p+r+d, where p is the intersection between the trajectory and the solid face.


This specular rebound approach is the most common treatment of the wall boundary condition for LSDM. However, Bahlali et al. (2020) remarked that this approach does not account for momentum exchange in the wall-tangential directions. Following Dreeben and Pope (1997) and Bahlali et al. (2020), the proper treatment of the boundary condition requires a precise evaluation of the near-wall Reynolds stresses, which is not provided by the turbulence model. The reflection method can also be generalized for an arbitrary surface (Brzozowska2013).

3.3.4 Concentration

In this study, fluid particles are considered massless, and, therefore, the Eulerian concentrations are obtained by counting the number of particles N(t) present in a sampling box of volume VxsΔysΔzs at each time step during an averaging time T. The average concentration is given by

(18) C ( x , y , z ) = 1 V T N ( t ) Δ t .

To compare concentration results between QES-Plume, analytical solutions, and experimental data, all concentrations are normalized using the following formula

(19) C * = C U H 2 Q ,

where C* is the normalized concentration; H is a reference height, often chosen as the height of the point of release; U is a reference velocity, often chosen as the velocity at the point of release; and Q is the source strength with units corresponding to the dimensional concentration.

4 Model evaluation

The performance of QES-Plume has been evaluated against two idealized test cases and a wind-tunnel test case for a 7×11 cubical array of buildings (Brown et al.2001). The primary goal is to examine the acceptability of QES-Plume with the newly implemented GLE solver when compared to available analytical solutions and wind-tunnel data.

4.1 Continuous release in uniform flow

The normalized concentration profiles from the 3D-GLE model computations have been compared to a classical Gaussian solution for an elevated continuous point-source release in a steady-state, horizontally homogeneous, neutrally stable atmosphere with constant wind speed and constant eddy diffusivity (Seinfeld and Pandis2016). The analytical concentration field is given by

(20) C ( x , y , z ) = Q 2 π U σ y σ z exp - y 2 2 σ y 2 exp - ( z - z s ) 2 2 σ z 2 ,

where Q is the source strength and U is the wind speed at source height zs. The plume standard deviations, σy and σz, in the crosswind and vertical direction, respectively, are given by σy=σvFt and σz=σwFt with F=[1+(t/Ti)1/2]-1 and Ti=(2.5u*/zi)-1, where t=x/U is the time of flight and zi is the boundary-layer height.

For this test case, the flow was prescribed by a horizontally and vertically uniform wind speed and friction velocity of U=2 m s−1 and u*=0.174 m s−1, respectively. The boundary-layer height was set to zi=1000 m. In addition, the turbulence model had to be simplified for the horizontally homogeneous, constant eddy diffusivity and neutral stability conditions following Rodean (1996). The stress tensor was computed using the following algebraic expressions:


The mean dissipation rate of TKE was given by

(26) ϵ ¯ = u * 3 κ z ( 1 - 8.5 z / z i ) 3 / 2 .

The particles were continuously released from a point source at height H=70 m at a rate of 200 particles per second with a time step of 1 s (Δt=1 s) for a duration of 2100 s. To obtain statistically stationary concentration estimates, the concentration was averaged over 1800 s with a starting time of 300 s after the beginning of the release. The physical domain was broken up into 20, 49, and 69 sampling boxes in the x , y, and z directions, respectively, over a domain size of 100 m × 140 m × 140 m. The simulation parameters are summarized in Table 2.

Table 2Simulation parameters for continuous release in a horizontally and vertically uniform flow test case.

Download Print Version | Download XLSX

Figure 4 shows the lateral and vertical normalized concentration profiles at three streamwise locations (x/H=0.464, x/H=0.821, and x/H=1.036) for QES-Plume and the Gaussian analytical solution (Eq. 20). The QES-Plume concentrations are in good agreement with the analytical solution with a coefficient of determination of r2=0.996, a root-mean-square error (RMSE) for the normalized concentration of 0.165, and a maximum relative error of 5.91 % over all profiles with x/H0.2. This threshold was chosen because the Gaussian plume model is known to be inaccurate close to the source (Stockie2011; Seinfeld and Pandis2016).

Figure 4Profiles of the normalized concentration for the classical Gaussian plume model (line) and the QES-Plume model (diamonds) at three different x/H locations. Panels (a)(c) are lateral profiles at z/H=1, and (d)(f) are vertical profiles at y/H=0. Concentrations are normalized following Eq. (19).


In addition, this test case corresponds to the dispersion from a point source in statistically stationary isotropic turbulence presented in Pope (2000, Chap. 12), following the Lagrangian approach introduced by Taylor (1921). To further validate QES-Plume, particle-trajectory statistics can be compared to this canonical example of turbulent dispersion. This approach shows that (i) trajectories close to the source consist essentially of straight-line motions, (ii) farther downstream the plume spread adheres to the square root of travel time, and (iii) the scaling variable is the Lagrangian timescale, τL. To produce well-converged Lagrangian statistics, 100 000 particles were released from the source, the domain was extended to Lx=12 km and Ly=8 km, the source was placed at yS=4 km, and the total simulation time was increased to 7800 s.

The timescale, τL, can be calculated based on the ratio of TKE to dissipation rate for the SLE (Eq. 4), with σ2=2/3k. However, this method, which yields τL≈136 s, is not valid for QES-Plume, which solves the 3D GLE. The timescale is computed using the autocorrelation function of the velocity fluctuation ρ. In this case, τL corresponds to the decorrelation timescale (or relaxation time) such that ρ(s)exp(-s/τL), where s represents the lags for the autocorrelation. Due to the limited length of the trajectories and large variability in the model, the autocorrelation function is averaged over all particles and the timescale is computed by fitting an exponential decay on lags smaller that 400 s. This method yields τL≈121 s, which is very close to the results from the SLE. According to Pope (2000, Chap. 12), the standard deviation of the positions follows the approximated form given by

(27) σ Y ( t ) u t , for t τ L , 2 u 2 τ L t , for t τ L ,

where σY is calculated on the spanwise position and u is the root-mean-square value of the velocity fluctuation from Eq. (5). Thus, the trajectories exhibit two distinct regimes with a region of linear spread and a region where spread follows a square root. Figure 5 presents the trajectories obtained in QES-Plume, for homogeneous isotropic turbulence from the uniform flow test case. The standard deviation of the spanwise position σY (Fig. 5a) matches both linear and square-root scaling regions. This behavior is visible in the sample of 200 trajectories which demonstrates the linear spread for t<τL (Fig. 5b) and the square-root spread for t>τL (Fig. 5c). Note that for t/τL>50, particles are exiting the domain, affecting the calculation of σY, hence the degradation of the scaling for large times. In conclusion, QES-Plume reproduced the canonical dispersion behavior from the Langevin model for homogeneous isotropic turbulence (Pope2000).

Figure 5Scaling of trajectories in the spanwise direction for homogeneous isotropic turbulence showing the linear spread for t<τL and the square-root spread for t>τL. Panel (a) shows the scaling of the standard deviation σY(t) of the horizontal spread. Panels (b) and (c) show a sample of 200 trajectories in the region of linear spread and square-root spread, respectively, where the black lines represent ±σY from Eq. (27).


4.2 Continuous release in a power-law atmospheric-boundary-layer flow

The next test case examines the performance of the QES-Plume model against an existing analytical solution for a continuous point-source release in a boundary-layer flow. The source was relatively close to the ground (H=4 m) to allow reflection of the emitted particles off the ground. The QES-Plume normalized concentration profiles are compared to the classical non-Gaussian solution (Huang1979; Brown et al.1993, 1997) for a steady-state, horizontally homogeneous, neutral atmospheric stability, power-law wind profile u(z)=azp and power-law eddy diffusivity Kz(z)=bzn. The analytical solution for the concentration is given by

(28) C ( x , y , z ) = Q 2 π σ y exp - y 2 2 σ y 2 ( z z s ) ( 1 - n ) / 2 b α x × exp - a ( z α + z s α ) b α 2 x I - ν 2 a ( z z s ) α / 2 b α 2 x ,

where α=2+p-n, ν=(1-n)/α and Iν is the modified Bessel function of first kind of order ν. The lateral standard deviation σy is given by σy=0.32x0.78 (Seinfeld and Pandis2016). In this test case, the turbulence stress gradients were present due to the velocity gradients in the vertical direction. The friction velocity is given by u*=κpazp. The stress tensor becomes τxx=(2.5u*)2, τyy=(2.3u*)2, τzz=(1.3u*)2, τxz=-(u*)2, and τxy=τyz=0. The mean dissipation rate of TKE is ε=5.7u*3/κz.

To obtain near-statistically stationary concentration estimates, 420 000 particles were continuously released from a point source at an emission rate of 200 particles per second with a time step of 1 s (Δt=1 s) for a duration of 2100 s. The power-law exponent for the velocity profile was 0.15 with a reference velocity U=5.90 m s−1 at a reference height of H=4 m. The concentration was averaged over 1800 s with a starting time of 300 s after the beginning of the release. The number of sampling boxes in the x, y, and z directions were 36, 59, and 20, respectively, over a domain size of 200 m × 100 m × 20 m. The source was specified to be at xs=20 m, ys=50 m, and zs=H=4 m. The simulation parameters are summarized in Table 3.

Table 3Simulation parameters for continuous release in a power-law boundary-layer flow test case.

Download Print Version | Download XLSX

Figure 6 shows the lateral and vertical normalized concentration profiles at four streamwise locations (x/H=4.03, x/H=10.97, x/H=19.31, and x/H=37.36) for the QES-Plume model and the non-Gaussian analytical solution. The model concentrations are in good agreement with the analytical solution. At all locations, the horizontal spread of the concentration along the centerline (z/H=1) is captured well by QES-Plume. At x/H=4.03, the concentration from QES-Plume matches the analytical solution except at the peak of the plume where there is a 19.69 % error in the maximum. At the second profile downstream of the release (x/H=10.97), the QES-Plume concentration peak has a significantly improved match with the theoretical model with only a 2.6 % overprediction. Further downstream, the prediction below the centerline (z/H=1) exhibits significant deviations from Eq. (28), and at x/H=19.31 and beyond the model fails to reproduce both the height and the value of the peak. The QES-Plume predictions exceed the theoretical concentration at the ground by a factor 3 at x/H=19.31. Additionally, the profile calculated from the model shifts to a monotonous profile around x/H=23.47, whereas the analytical solution does not exhibit the same behavior (Fig. 6h). The non-Gaussian solution is known to underpredict the concentration close to the ground (Brown et al.1993, 1997) likely accounting for a significant number of the observed deviations in this region.

Figure 6Profiles of the normalized concentration for the non-Gaussian plume analytical solution (line) and the QES-Plume model (diamonds) at four different x/H locations. The top panels are lateral profiles at z/H=0.875, and the bottom panels are vertical profiles at y/H=0. Concentrations are normalized following Eq. (19).


4.3 Array of cubical buildings

The idealized test cases are useful to evaluate QES-Plume's GLEs solution methodology because they control for external factors that impact the quality of the dispersion model including turbulence parameterization, mean velocity specification, and boundary conditions. The limitation is that they do not fully engage the GLEs because they do not activate all components of the stress and velocity gradient tensors. To more fully examine the performance of QES-Plume’s GLE implementation, it is compared to dispersion data from a wind-tunnel experiment for a 7×11 array of cubical buildings. The wind-tunnel experiment was conducted in a United State Environmental Protection Agency (EPA) meteorological wind tunnel (Brown et al.2000, 2001). For concentration estimates, high-purity ethane (C2H6, chemically pure grade, minimum purity 99.5 mole percent) was used as a tracer, which is slightly heavier than air (molecular weight 30). This tracer may be regarded as neutrally buoyant owing to the high turbulence level and the release rate of the tracer. A perforated plastic sphere with a diameter of 10 mm was used for continuously releasing the tracer at ground level behind the first centerline building of the 7×11 array. Figure 7a presents the layout of the array within the domain. Figure 7b shows the location of the source (diamond) and the location of vertical concentration profile measurements (circle) in the wind tunnel. In addition, spanwise transverse concentration measurements were made at a height of z/H=0.1 at the same x/H locations, except at x/H=2.5 because this location interferes with the buildings. The minimum quantifiable non-dimensional concentration for the wind-tunnel measurements is assumed to be 10−3 (following the recommendation from Chang and Hanna2004).

Figure 7Simulation setup for the 7×11 array of cubical buildings. Panels are (a) the full simulation domain of 26.6H×26.6H and (b) a zoomed-in section of the domain showing the location of the source behind the first building at z/H=0.1 (diamond) and the locations of wind-tunnel measurements at y/H=0 and y/H=1 for multiple downstream locations (circles).

For this test case, QES-Plume was driven using the flow field computed by QES-Winds. Following Singh et al. (2008), the street-canyon parameterization needed to be modified to address some of the shortcomings of the original street-canyon parameterization from Röckle (1990) used by QUIC-URB. The modified model adds a blending region at the edge of the canyons to resolve issues related to erroneous gradients, a feature of the flow which is critical to obtain the correct particle dispersion in and out of the street canyons. The turbulence fields were computed by QES-Turb, with the addition of a non-local background mixing (Sect. 3.2).

A total of 1 950 000 particles were released continuously for 3900 s (i.e., Q=500 particles per second) from a point source to obtain near-statistically stationary concentration estimates with QES-Plume. The concentration was averaged over 3600 s with a starting time of 300 s after beginning the release. The size of the sampling boxes in x, y, and z directions was set to 1.5, 1.5, and 1.5 m, respectively, over a domain size of 400 m ×400 m ×60 m. The source was placed behind the first centerline building of the array at x/H=1.067, y/H=0, and z/H=0.067 to be in agreement with the wind-tunnel experiments. During the simulation, no rogue trajectories were detected, confirming the stability of the integration scheme in a complex environment.

Table 4Simulation parameters for the array of cubical buildings test case.

Download Print Version | Download XLSX

Figure 8Non-dimensional concentration field from QES-Plume: (a) vertical slice at y/H=0 and (b) horizontal slice at z/H=0.3. The concentration has been normalized using Eq. (19), the color axis is in logarithmic scale, the contour lines represent each decade of the concentration scale, and the gray contour line (C*=10-3) corresponds to the minimum concentration threshold measurement in the wind-tunnel experiment.


Cross sections of the concentration results of dispersion through the building array are shown in Fig. 8. Many processes characteristic of urban dispersion at the neighborhood scale are illustrated (Belcher2005). To simplify the discussion, the spaces in between the buildings along the x and y directions are designated as street channels and street canyons, respectively. Importantly, the minimal non-dimensional concentration computed by the model is C*=1.5×10-4. For comparison with the wind-tunnel data, values smaller than the experimental minimum threshold of 10−3 – marked by the gray contour line in the figure – are regarded as zeros. Figure 8a shows a vertical slice along the y/H=0 centerline and Fig. 8b a horizontal slice at z/H=0.3. As particles advect downstream, the plume spreads laterally into the adjacent street canyons. Between the first and the third street canyons, the plume width grows linearly in the first street channel (y/H=±1.5). This linear growth stage is characteristic of the near-field region (Belcher2005). In addition, the recirculating flow present in the street canyons traps particles within the regions between buildings. High concentrations in street canyons behave similar to a source, leading to a phenomenon referred to as a secondary source (Belcher2005). Above the array, the vertical spread of the plume is enhanced by the vertical transport out of each street canyon. Farther downstream, the plume width is nearly constant once the particles have reached the second street channel (y/H=±2.5). Belcher (2005) calls this region the far field; it is characterized by a slower rate of spread past the first few rows of buildings. Street intersections also play an important role in the observed dispersion pattern. For example, T-junctions lead to very different behavior than four-way intersections. In the former, topological dispersion occurs due to dividing streamlines and is not observed in the latter type of intersections (Belcher2005). In the array considered in the current simulation, there is no topological dispersion at the street intersections. Additionally, the near-field and far-field regions seem akin to the different scaling regions presented in Fig. 5, suggesting that the region between 1<x/H<6 and -1.5<y/H<1.5 corresponds to the region of linear growth with τL<1. However, further investigation is needed and would be beyond the scope of this paper.

To begin a quantitative comparison, QES-Plume results are plotted with respect to the concentrations from the wind-tunnel dataset in Fig. 9 for vertical measurements along the centerline (y/H=0, panel a) and in between the row of buildings (y/H=1, panel b) as well as transverse measurements close to the ground in the street canyons (z/H=0.1, panel c) for multiple downstream locations.

Figure 9Comparison of QES-Plume concentration (lines) with wind-tunnel data (squares) for the 7×11 cubical array study. The panels in (a) show vertical profiles at y/H=0, the panels in (b) show vertical profiles at y/H=1, and the panels in (c) show horizontal profiles at z/H=0.1. The downstream locations are reported in the panels. The concentrations have been normalized using Eq. (19). Missing panels correspond to locations where the measurements would have coincided with buildings.


In the first street canyon (first graph of Fig. 9a), the concentration is measured in the middle of the street canyon, close to the source. The experimental profile increases from the ground to its maximum near the building height (z/H=0.85) and then decreases with height. In contrast, the QES-Plume profile contains two maxima: one in the middle of the street canyon (z/H=0.40) and one above the top of the building (z/H=1.1). The location of the street-canyon vortex and the shear-layer growth are critical to the dynamics of dispersion in the first street canyon because of the source location. Small deviations can have a great consequence for the concentration profile. A stronger back flow in the lower part of the street canyon is observed in the wind tunnel, leading to the higher location of the maximum concentration compared to QES-Plume. The vertical spread above the building height is comparable between the QES-Plume and the wind-tunnel data (see Appendix A1). Overall, the model prediction is acceptable considering the strong sensitivity to source location with an RMSE of the normalized concentration of 2.305 (relative RMSE of 0.272). The RMSE and relative RMSE of all the different profiles are compiled in Table 5.

In the second street canyon (x/H=3.5), the experimental profile is mostly constant between the ground and z/H=0.5, then increases to reach its maximum at z/H=1.1, and finally decreases to the top of the plume at z/H=2.0. The profile calculated by QES-Plume shows a significantly different shape with two large spikes at z/H=0.8 and z/H=1.1. The QES concentration close to the ground is 50 % higher than the experimental concentration and is not constant in the lower part of the street canyon. It decreases until z/H=0.7 and then spikes at z/H=0.8. Interestingly, there is only a 2 % difference between the mean concentrations within the street canyon. Hence, even if QES does not capture the shape of the profile, the overall average concentration is estimated correctly. The different spikes observed in the profile are linked to the parameterization used by QES-Winds, an overestimation of the velocity variances, and the presence of sharp gradients at the top of the street canyon (0.7<z/H<1.1). A similar shape is observed for the concentration profile at the center of the third street canyon (x/H=5.5) with an RMSE equal to 0.098. In addition, the average concentrations in the street canyon are also within a 2 % margin. The shape of the profile in the fifth street canyon (x/H=9.5) is very similar to the profile discussed previously but with the experimental values consistently higher than values calculated by QES-Plume. Although the values match close to the surface, the mean experimental concentration in the canyon is 20 % higher. Likewise, QES-Plume significantly underestimates the concentration above the buildings from z/H=1 to z/H=2. The rest of the profiles are in good agreement. At this location, the RMSE is equal to 0.068. Past the last building (x/H=13.5), the shape of the calculated profile does not exhibit large spikes and is closer to the experimental data. However, the model consistently underestimates the concentration at this location, leading to an RMSE of 0.075. In the far wake (x/H=16.8), the model underestimates the concentration compared to the experimental data with an RMSE of 0.048. The oscillations observed at the last two locations correspond to the transition between the different wake regions used by QES-Winds. For example, the spike at z/H0.6 for x/H=16.8 corresponds to the edge of the far-wake parameterization.

The middle panel (Fig. 9b) shows the vertical profiles in the first street channel (y/H=1). At x/H=1.5, the model underestimates the concentration, predicting almost zero concentration at all heights. This is an indication of the lack of lateral dispersion into the channel. At x/H=2.5, the modeled concentration is no longer 0 but still underestimates the data, with concentrations about half of the wind-tunnel measurements for z/H<0.5 and almost 0 for z/H>0.5, resulting in an RMSE=0.197. Farther downstream, at x/H=3.5, QES-Plume shows significant lateral dispersion and is in good agreement with the experimental data close to the ground (z/H<0.2). However, the simulated concentrations decrease much faster with height compared to the wind-tunnel data for z/H>0.2 and significantly underestimate concentrations for z/H>0.5 to z/H<1.5, resulting in an RMSE value of 0.130. At x/H=5.5, the measured concentration profiles are below z/H<0.5 and decrease to reach 0 at z/H1.5. The simulated concentrations exhibit a similar behavior, even if the lower plateau is smaller (z/H<0.3). They also reach 0 at roughly the same height, resulting in an RMSE of 0.070. At x/H=9.5, QES-Plume predictions are in relatively good agreement with wind-tunnel measurement, with RMSE=0.034. However, the measured concentration profile is not monotonous, unlike the QES-Plume profile, with a maximum at z/H0.75. The model failed to capture this feature of the profile. Past the last building, at x/H=13.5 and x/H=16.8, the simulated concentrations show good agreement with the wind-tunnel measurements for all heights with RMSE=0.018 and RMSE=0.028, respectively.

Lateral measurements of the concentration are presented in Fig. 9c. In the first street canyon, the concentration from QES-Plume contains two spikes at y/H=±0.5 which are not observed in the experimental data. These spikes are located at the edges of the upstream building. The QES-Winds canyon model does not extend past the edge of the building (Singh et al.2008), which leads to sharp gradients at these locations, preventing some lateral dispersion. The concentration at the center of the street canyon and the width of the plume match between the model and the experiment. The RMSE over the whole profile is 0.335. Similar behavior is observed in the second and third street canyons with RMSE=0.152 and RMSE=0.077, respectively. Farther downstream, the fifth street canyon is past the near field and the plume has spread beyond the second street channel and into the third. At this point, the modeled and measured concentrations are in good agreement with an RMSE=0.034. The oscillations observed in the simulated data are linked to the different parameterization used by QES-Winds. Past the last building, the plume has spread from y/H<-3 to y/H>3. Although, QES-Plume still has some underestimation of the concentration, the total width of the plume matches well.

Some of the discrepancies found in Fig. 9b are linked to the sharp concentration drop-off observed in the horizontal profiles which, for locations in the first half of the array where the plume width is growing linearly (near field), can be found at the middle of the first channel in between the buildings (y/H=1). The location of this drop-off is very sensitive to velocity and turbulence fields. In the second half of the array, a sharp drop-off is not observed in the horizontal concentration profiles. Overall, QES-Plume predicts the width of the plume and the value of the concentration adequately, both within the array and downstream of the last building. The mean relative RMSE over all the transects and profiles is 15.6 %.

Table 5Summary of QES-Plumes's RMSEs and relative RMSEs at various locations for the 7×11 cubical array study.

Download Print Version | Download XLSX

Figure 10Signed relative error (RE) between the concentration field from QES-Plume and the wind-tunnel measurements where positive values (red markers) represent an overestimation by the model. The panels in the figure are (a) the vertical slice at y/H=0 and (b) the horizontal slice at z/H=0.1. The contour lines represent each decade of the concentration form 10 to 10−3, and the gray contour line (C*=10-3) corresponds to the minimum concentration threshold measurement in the wind-tunnel experiment.


Figure 10 shows the signed relative error between the QES-Plume and the wind-tunnel data, where positive values represent overestimation by the model. Large relative errors can be observed at the edge of the plume, where the absolute concentrations are close to the lower limit of the measurements. The deviations observed in the first street canyon (1x/H2 and -0.5y/H0.5) are linked to the sensitivity of the mean flow model related to the location of the street canyon vortex, emphasizing the discussion related to the first graph of Fig. 9a. Both vertical and horizontal slices show that large relative errors are found in the shear zones behind buildings, in particular at y/H=1.5 between x/H=3 and x/H10. These observations are related to the parameterization used by the wind model, where different zones are defined based on geometrical considerations (Singh et al.2008). A comprehensive comparison between the mean velocity field and velocity variances from QES to the wind-tunnel data is presented in Appendix A. Moreover, Figs. 9 and 10 illustrate the inherent asymmetry present in the data and model due to the variability in the system, where small deviations between the left (y/H>0) and right (y/H<0) sides can yield significant relative errors between the data and model.

Finally, the local-mixing turbulence model relies heavily on the magnitude of the local velocity gradients (see Eqs. 10 and 6). This is problematic in regions where velocity gradients are small and the model predicts negligible turbulence. On the other hand, regions with sharp velocity gradients lead to unrealistically large stresses and TKE, due to the use of a diagnostic wind model and a local-mixing model, as illustrated in Appendix A. A future improvement in the turbulence model would be to limit the magnitude of the velocity gradients, especially the vertical gradient of streamwise velocity. In addition, the Appendix presents the importance of adding non-local mixing to the turbulence model to enhance lateral and vertical spread. Although this observation is related to wind-tunnel data, similar considerations hold for atmospheric flows, where background turbulence levels need to be evaluated from sources such as meteorological measurements or weather prediction models.

Figure 11Paired scatterplot of the wind-tunnel concentration data and QES-Plume concentrations for the 7×11 cubical array study. The one-to-one line is represented by the black line, a factor of 2 by the dashed lines, a factor of 5 by the dot–dashed lines, and a factor of 10 by the dotted lines.


In summary, QES-Plume is capable of reproducing concentration levels in this complex mock-urban setting despite weak performance in the first street canyon. In particular, the results show the expected linear behavior near the source and constant region farther away as reported by Belcher (2005). Similarly, critical metrics such as plume width and height are reproduced by the model. Specifically, out of the 894 measurements in the experimental data, 156 are below the concentration threshold and QES-Plume matches 154 of those points, yielding 99 % matched zeros. Figure 11 shows a paired scatterplot of the vertical and horizontal concentration profiles. The predicted concentrations by QES-Plume are represented on the abscissa and the wind-tunnel concentration data are represented on the ordinate of the scatterplot. This plot shows that the model underestimates the concentration with 79 % of QES concentrations being lower than the experimental data. However, 59 % of the predicted concentrations fall within a factor of 2 (with an RMSE of 0.300), 76 % within a factor 5 (with an RMSE of 0.529), and 80 % within a factor 10 (with an RMSE of 0.519). Of the 20 % of data points outside of the factor 10 threshold, QES-Plume failed to compute any concentration at 80 % of these locations, corresponding to a 15 % false negative rate. The statistics are summarized in Table 6. In addition, the RMSEs in the factor 5 and factor 10 are dominated by five points in the profile at y/H=0 and x/H=1.5. If the corresponding profile is removed from the calculation of the RMSEs, the values become 2 to 5 times smaller. All the data points outside of the factor 10 class correspond to normalized concentrations smaller than 5×10-2 and coincide with the edge of the plume where both measurement and simulated values have large uncertainties.

Table 6Summary statistics comparing QES-Plume results with the wind-tunnel data for the 7×11 cubical array study.

Download Print Version | Download XLSX

Finally, QES-Plume exhibits excellent computational performance even without parallelization. The 7×11 array simulation runs for 3900 s of simulation time and contains about 88 000 active particles at each time step. The average total execution time on an AMD EPYC 7543 32-core processor is around 1100 s. This is equal to an execution time of less than 4×10-6 s per particle per time step.

5 Discussion

The method implemented to partially solve the stiffness problem from the GLEs (Sect. 3.3) guarantees numerical stability of the velocity fluctuations, which means that no energy is artificially added to the system. However, large gradients in the velocity field or regions of unrealistically large turbulence quantities can lead to nonphysical fluctuations. Equation (15) contains both temporal and spatial derivatives of the stress tensor, and, therefore, the model remains sensitive to large increments in these two quantities. The same observation can be made for the mean dissipation rate of TKE. This quantity needs to stay within realistic bounds, especially since it multiplies the random term in the model. Finally, the particle positions are updated using Euler's forward scheme, which does not guarantee numerical stability. The authors have not encountered stability issues related to the last step of the model if the fluctuations obtained from the GLE take realistic values. The numerical stability issues associated with Eq. (17) do not lead to a degradation of the solution as truncation errors are not propagated in time.

The test cases presented by Bailey (2017) showed that the model produced good results even for large time steps. However, some aspects of the model are sensitive to time step sizes. Concentration calculations are impacted if the size of the collection box is too small for the particle motion. Similarly, the reflective boundary condition can lead to errors. For example, if the particle displacements exceed the size of some obstacles in the domain, the particles might appear to have warped through the obstacle.

The local-mixing turbulence model relies heavily on the magnitude of the local velocity gradients (see Eqs. 10 and 6). This is problematic in regions where velocity gradients are small and the model predicts negligible turbulence. On the other hand, regions with sharp velocity gradients lead to unrealistically large stresses and TKE. Appendix A presents the mean velocity field and velocity variances for the array of cubical buildings. Comparison with the wind-tunnel data is also included in the Appendix. The results illustrate the consequences of using a local-mixing model with a diagnostic wind solver. A future option would be to limit the magnitude of the velocity gradients, especially the vertical gradient of streamwise velocity. In addition, the Appendix presents the importance of adding non-local mixing to the turbulence model to enhance lateral and vertical spread. Although this observation is related to wind-tunnel data, similar considerations hold for atmospheric flows, where background turbulence levels need to be evaluated from sources such as meteorological measurements or weather prediction models.

6 Conclusions

Due to the presence of a large number of terms in the GLEs, SLEs are employed in most mainstream dispersion models. The presence of numerical instabilities due to the stiffness of the GLEs has also been a problem for the explicit integration of the GLEs into Lagrangian dispersion models. Although commonly used numerical methods for solving the SLEs are still numerically unstable, the SLEs are considered slightly more stable compared to the GLEs due to the drastic reduction in the number of terms. This paper discussed the implementation of the GLEs with the implicit time-integration method from Bailey (2017) to alleviate the stiffness problem from the GLEs and eliminate rogue trajectories (Yee and Wilson2007). This implicit scheme was implemented in the QES framework and is dynamically coupled with a mass-consistent wind model and a local-mixing turbulence model.

QES-Plume has been validated against analytical solutions in idealized conditions where the model yielded good results. In particular, the overall maximum relative error was under 6 %, and the model captured horizontal and vertical plume width accurately. Discrepancies between the model and the analytical solutions matched known shortcomings of the latter. During the comparison with experimental data, QES-Plume performed well, highlighted by the concentration contours showing good levels of lateral and vertical dispersion, as well as a well-mixed plume in the street channels. The results showed 99 % matched zeros, a factor 2 concentration prediction of 59 %, and only a 15 % false negative rate. However, the 7×11 test case emphasized the sensitivity to the mean-wind and turbulence models, as emphasized by Bahlali et al. (2019). In particular, the street-canyon model proposed by Singh et al. (2008) is used in this study and yielded better mean flow, turbulence, and dispersion results compared to the base model proposed by Röckle (1990). Moreover, the local-mixing-length turbulence model necessitated the addition of a constant to the diagonal elements of the stress tensor to enhance the turbulence mixing within the street canyons and in the free stream. The constant of 0.3 m2 s−2 was added to match the turbulence level of the experimental data within the free stream. However, background turbulence levels should be adjusted to reflect realistic conditions and can be evaluated from meteorological measurements or weather prediction models. In general, the results are dependent on the wind and turbulence models. Further improvement in QES-Turb are planned to address some of the shortcomings observed in the turbulence fields, such as sharp gradients in the shear zone. In addition, QES-Plume does not require QES-Winds and QES-Turb; the code can be used with inputs from other models.

Finally, QES-Plume was implemented from the ground up as an object-oriented C++ code and has demonstrated excellent computational performance. Future versions of QES-Plume are likely to use a GPU-based implementation to enable much faster than real-time simulations. The potential use cases of a model like QES-Plume are numerous. For example, the model can be used to run simulations for decision makers for tabletop exercises or to study particulate dispersion in complex environments, such as spore or smoke transport in agricultural fields, as well as urban pollution and air quality.

Appendix A: Modeled mean-wind and turbulence fields for the 7×11 array of cubical buildings

In this section, the performance of QES-Winds and QES-Turb is briefly discussed in the context of the 7×11 array of cubical buildings. As dispersion results rely heavily on mean flow and turbulence stress fields, it is important to understand the quality of the fields used to drive the dispersion model. This section is not intended as a strict validation of the mean flow model (see Singh et al.2008, for a validation of the street canyon model) or the local-mixing turbulence model.

A1 Mean flow field

Figure A1 shows three profiles of the streamwise and vertical velocity upwind of the first building. The profiles at x/H=-1 illustrate that the inlet velocity field used in QES matches the experimental data. An upstream recirculation zone starts at x/H-0.8 and grows to reach z/H0.7 on the face of the building. The QES profiles in this upwind region show the deceleration of the streamwise velocity and the displacement of the streamlines, in good agreement with the experimental data, on average. The u-velocity profile contains oscillations at the edge of the parameterization region, which is typical of the kind of model used. The vertical velocity at x/H=-0.3 illustrates the sensitivity to parameterization choices of the location of topological flow features such as the upwind vortex (Hayati et al.2019). The experimental data indicate that this vortex is located further away from the building compared to the location predicted by QES. Hence, the w velocity is positive in the lower part of the profile (0<z/H<0.5).

Figure A2 compares the profiles at the top of the first building along the centerline of the array. Profiles at x/H=0.1 illustrate that the flow field predicted by QES at this location is in good agreement with the experimental data. At x/H=0.5, the streamwise velocity overshoots the experimental data at z/H1.25 as the rotation strength predicted by the rooftop vortex parameterization is higher than observed in the wind-tunnel data. Similarly, the parameterization overpredicts the w velocity close to the back edge of the building x/H=0.9.

Profiles at three locations in the first street canyon are presented in Fig. A3. The u-velocity profiles from QES contain sharp vertical gradients and a very shallow shear layer at the top of the street canyon (z/H=1). In contrast, the experimental data have less steep gradients and a deeper shear layer. The street canyon model proposed by Singh et al. (2008) accounts for shear-layer growth from the edge of the building. However, in its current formulation, the model only parameterizes the part of the shear layer below z/H=1 and does not extend vertically past the building top. Still, the street canyon model does a relatively good job of predicting the value of the streamwise velocity within the canyon. The vertical-velocity profiles illustrate that QES underestimates w in the second half of the street canyon. These observations indicate that the intensity of the street canyon vortex might be underestimated and that its location might not be predicted correctly.

Figure A4 shows profiles downstream of the last building. Similar to the street canyon, the u-velocity profiles (x/H=13.1 and x/H=13.5) illustrate sharp vertical gradients close to the top of the building with a very shallow shear layer. The gradients are less steep in the experimental data and the shear layer is deeper, extending from z/H0.8 to z/H1.5. These observations are consequences of the formulation of the wake model as it does not extend vertically past the top of the building and does not model a shear layer that grows above the top of the building. At x/H=16.5, the downstream recovery calculated by QES is faster than observed in the experimental data. Finally, the w-velocity profiles show good agreement for all profiles.

Overall, the mean flow along the centerline of the array is in good agreement with the experimental data and illustrates that the model proposed by Singh et al. (2008) has been correctly implemented in QES.

Figure A1Comparison of QES-Winds (lines) with wind-tunnel data (squares) for the streamwise velocity u (a) and the vertical velocity w (b) upwind of the first building. All profiles are taken along the centerline (y/H=0).


Figure A2Comparison of QES-Winds (lines) with wind-tunnel data (squares) for the streamwise velocity u (a) and the vertical velocity w (b) at the roof-top level of the first building. All profiles are taken along the centerline (y/H=0).


Figure A3Comparison of QES-Winds (lines) with wind-tunnel data (squares) for the streamwise velocity u (a) and the vertical velocity w (b) within the first street canyon. All profiles are taken along the centerline (y/H=0).


Figure A4Comparison of QES-Winds (lines) with wind-tunnel data (squares) for the streamwise velocity u (a) and the vertical velocity w (b) downwind of the last building. All profiles are taken along the centerline (y/H=0).


A2 Turbulence fields

In this section, the velocity variances from the QES-Turb model and the wind-tunnel data are compared at selected locations to qualitatively evaluate the performance of the model.

Figure A5 indicates that QES-Turb overestimates both the streamwise-velocity and the vertical-velocity variances upwind of the first and second buildings. Within the upwind recirculation zone (-0.8<x/H0 and 0z/H<0.7), the model performs better. In addition, the signature of the oscillation observed in the velocity field can be seen in the variances as expected from Eq. (6).

Figure A6 compares the QES-Turb model and the wind-tunnel data for the velocity variances at the centerline over the rooftop of the first and second buildings. The profiles illustrate that QES overestimates the streamwise-velocity variance above the first building (1<z/H<2). The vertical-velocity variance is overestimated near the leading edge of the first building. The model yields better results past the middle of the building. Over the second building, both variances are better reproduced by the model. The same observation can be made at all available downstream roof-top locations (not shown). All profiles show that the modeled free-stream turbulence levels match the experimental data well.

Figures A7 and A8 compare the QES-Turb model and the wind-tunnel measurements for the velocity variances at the centerline of the first street canyon and downwind of the last building. The calculated variances are in good agreement with the experimental data for all profiles far above the buildings. The same observation can be made in the lower half of the street canyon (z/H<0.5). However, the region in between shows that the turbulence model overestimates the variances. In particular at z/H=1, where the transition between the parameterizations occurs, the variances are overestimated due to the steep gradients presented in Sect. A1. Similarly, in the wake and recovery zones downstream of the last building, both variances contain a large peak at this transition.

In summary, the turbulence model used to drive QES-Plume shows good agreement with the wind-tunnel data. Importantly, the local-mixing model needed the addition of a constant Cnlm to enhance the turbulence mixing within the street canyons and in the free stream (see Sect. 3.2). Without the non-local-mixing constant, particles had the tendency to follow closely the mean wind near the source until they are ejected out of the first street canyon without much mixing within the street canyon and channel. The constant was selected to match the turbulence level of the experimental data within the street canyons and in the free stream. Although the turbulence levels at the edges of the street canyons are increased with the non-local-mixing constant added, the turbulence levels within the street canyons matched the test data well. The result of this modification was to augment the magnitude of the turbulence within the street canyons while maintaining the influence of the velocity gradients in the stress tensor.

Figure A5Comparison of QES-Turb (lines) with wind-tunnel data (squares) for the streamwise-velocity variance σu (a) and the vertical-velocity variance σw (b) upwind of the first building. All profiles are taken along the centerline (y/H=0).


Figure A6Comparison of QES-Turb (lines) with wind-tunnel data (squares) for the streamwise-velocity variance σu (a) and the vertical-velocity variance σw (b) at the roof top of the first and second buildings. All profiles are taken along the centerline (y/H=0).


Figure A7Comparison of QES-Turb (lines) with wind-tunnel data (squares) for the streamwise-velocity variance σu (a) and the vertical-velocity variance σw (b) within the first street canyon. All profiles are taken along the centerline (y/H=0).


Figure A8Comparison of QES-Turb (lines) with wind-tunnel data (squares) for the streamwise-velocity variance σu (a) and the vertical-velocity variance σw (b) downwind of the last building. All profiles are taken along the centerline (y/H=0).


Code availability

The Quick Environmental Simulation (QES) software has been developed as part of a collaboration between the University of Utah, the University of Minnesota Duluth, and Pukyong National University (, Margairaz et al.2022a). The software is written in C++ and NVIDIA's CUDA. The latest version of QES is publicly accessible on GitHub (, last access: 6 December 2022) under the GNU General Public License (version 3). More information is available on GitHub about the installation of QES with or without CUDA. QES version v2.0.0 was used to produce the results used in this paper and is archived on Zenodo (Margairaz et al.2022a).

Data availability

The data from the EPA for the 7×11 array wind-tunnel experiment from Brown et al. (2000) are publicly available. The dataset can be obtained upon request to the EPA ( The dataset can also be downloaded from the QES-Public GitHub repository (, last access: 6 December 2022).

Author contributions

Conceptualization was done by RS and ERP. The code was written by JAG, LA, and FM. Test cases were developed by BS and FM. The analysis and interpretation of the data were carried out by FM. The figures were produced by FM. The original draft of the paper was written by FM, with edits, suggestions, and revisions provided by JAG, BS, RS, and ERP. Grants supporting the project leading to this publication were awarded to RS and ERP.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to acknowledge the funding provided by the National Institute of Environment Research (NIER), funded by the Ministry of Environment of the Republic of Korea; the United States Department of Agriculture National Institute for Food and Agriculture (USDA-NIFA) Specialty Crop Research Initiative (SCRI); and the United States Department of Agriculture Agricultural Research Service (USDA-ARS) through a research support agreement.

Financial support

This research has been supported by the National Institute of Environmental Research (grant no. NIER-SP2019-312), the United States Department of Agriculture National Institute for Food and Agriculture (grant nos. USDA-NIFA-SCRI 2018-03375 and USDA-NIFA-SCRI 2020-02656), and the United States Department of Agriculture Agricultural Research Service (research support agreement no. 58-2072-0-036).

Review statement

This paper was edited by Leena Järvi and reviewed by Bertrand Carissimo and Jérémy Bernard.


Archambeau, F., Méchitoua, N., and Sakiz, M.: Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows – Industrial Applications, International Journal on Finite Volumes, 1, (last access: 26 September 2023), 2004. a

Aylor, D.: Aerial Dispersal of Pollen and Spores, The American Phytopathological Society, St. Paul, Minnesota, USA,, 2017. a

Aylor, D. E.: Spread of plant disease on a continental scale: role of aerial dispersal of pathogens, Ecology, 84, 1989–1997,, 2003. a

Bahlali, M. L., Dupont, E., and Carissimo, B.: A hybrid CFD RANS/Lagrangian approach to model atmospheric dispersion of pollutants in complex urban geometries, Int. J. Environ. Pollut., 64, 74–89,, 2018. a

Bahlali, M. L., Dupont, E., and Carissimo, B.: Atmospheric dispersion using a Lagrangian stochastic approach: Application to an idealized urban area under neutral and stable meteorological conditions, J. Wind Eng. Ind. Aerod., 193, 103976,, 2019. a

Bahlali, M. L., Henry, C., and Carissimo, B.: On the Well-Mixed Condition and Consistency Issues in Hybrid Eulerian/Lagrangian Stochastic Models of Dispersion, Bound.-Lay. Meteorol., 174, 275–296,, 2020. a, b, c

Bailey, B. N.: Numerical Considerations for Lagrangian Stochastic Dispersion Models: Eliminating Rogue Trajectories, and the Importance of Numerical Accuracy, Bound.-Lay. Meteorol., 162, 43–70,, 2017. a, b, c, d, e, f, g, h, i, j

Belcher, S. E.: Mixing and transport in urban areas, Philos. T. Roy. Soc. A, 363, 2947–2968,, 2005. a, b, c, d, e, f

Bozorgmehr, B., Willemsen, P., Gibbs, J. A., Stoll, R., Kim, J.-J., and Pardyjak, E. R.: Utilizing dynamic parallelism in CUDA to accelerate a 3D red-black successive over relaxation wind-field solver, Environ. Modell. Softw., 137, 104958,, 2021. a, b

Britter, R. E. and Hanna, S. R.: Flow and Dispersion in Urban Areas, Annu. Rev. Fluid Mech., 35, 469–496,, 2003. a, b

Brown, M., Lawson, R., Decroix, D., and Lee, R.: Mean Flow and Turbulence Measurement around a 2-D Array of Buildings in a Wind Tunnel, in: 11th Joint Conference on the Applications of Air Pollution Meteorology with the AWMA, Long Beach, CA, 9–14 January 2000, (last access: 26 September 2023), 2000. a, b

Brown, M., Lawson, R., DeCroix, D., and Lee, R.: Comparison of Centerline Velocity Measurements Obtained Around 2D and 3D Building Arrays in a Wind Tunnel, (last access: 26 September 2023), 2001. a, b

Brown, M. J., Arya, S. P., and Snyder, W. H.: Vertical Dispersion from Surface and Elevated Releases: An Investigation of a Non-Gaussian Plume Model, J. Appl. Meteorol., 32, 490–505,<0490:vdfsae>;2, 1993. a, b

Brown, M. J., Arya, S. P., and Snyder, W. H.: Plume: Descriptors derived from a non-Gaussian concentration model, Atmos. Environ., 31, 183–189,, 1997. a, b

Brown, M. J., Gowardhan, A. A., Nelson, M. A., Williams, M. D., and Pardyjak, E. R.: QUIC transport and dispersion modelling of two releases from the Joint Urban 2003 field experiment, Int. J. Environ. Pollut., 52, 263–287,, 2013. a

Brunet, Y.: Turbulent Flow in Plant Canopies: Historical Perspective and Overview, Bound.-Lay. Meteorol., 177, 315–364,, 2020. a, b, c

Brzozowska, L.: Validation of a Lagrangian particle model, Atmos. Environ., 70, 218–226,, 2013. a, b

Byun, D. and Schere, K. L.: Review of the Governing Equations, Computational Algorithms, and Other Components of the Models-3 Community Multiscale Air Quality (CMAQ) Modeling System, Appl. Mech. Rev., 59, 51–77,, 2006. a

Carissimo, B., Castelli, S. T., and Tinarelli, G.: JRII special sonic anemometer study: A first comparison of building wakes measurements with different levels of numerical modelling approaches, Atmos. Environ., 244, 117798,, 2021. a

Castelli, S. T., Ferrero, E., and Anfossi, D.: Turbulence Closures In Neutral Boundary Layer Over Complex Terrain, Bound.-Lay. Meteorol., 100, 405–419,, 2001. a

Chang, J. C. and Hanna, S. R.: Air quality model performance evaluation, Meteorol. Atmos. Phys., 87, 167–196,, 2004. a

Cimorelli, A. J., Perry, S. G., Venkatram, A., Weil, J. C., Paine, R., Wilson, R. B., Lee, R. F., Peters, W. D., and Brode, R. W.: AERMOD: A Dispersion Model for Industrial Source Applications. Part I: General Model Formulation and Boundary Layer Characterization, J. Appl. Meteorol., 44, 682–693,, 2005. a

Dreeben, T. D. and Pope, S. B.: Probability density function and Reynolds‐stress modeling of near‐wall turbulent flows, Phys. Fluids, 9, 154–163,, 1997. a

Du, S.: Universality of the Lagrangian Velocity Structure Function Constant (C0) Across Different Kinds of Turbulence, Bound.-Lay. Meteorol., 83, 207–219,, 1997. a

Finnigan, J.: Turbulence in Plant Canopies, Annu. Rev. Fluid Mech., 32, 519–571,, 2000. a, b

Gowardhan, A. A., Brown, M. J., and Pardyjak, E. R.: Evaluation of a fast response pressure solver for flow around an isolated cube, Environ. Fluid Mech., 10, 311–328,, 2010. a

Gowardhan, A. A., Pardyjak, E. R., Senocak, I., and Brown, M. J.: A CFD-based wind solver for an urban fast response transport and dispersion model, Environ. Fluid Mech., 11, 439–464,, 2011. a

Gowardhan, A. A., McGuffin, D. L., Lucas, D. D., Neuscamman, S. J., Alvarez, O., and Glascoe, L. G.: Large Eddy Simulations of Turbulent and Buoyant Flows in Urban and Complex Terrain Areas Using the Aeolus Model, Atmosphere, 12, 1107,, 2021. a

Hanna, S. R., Britter, R., and Franzese, P.: A baseline urban dispersion model evaluated with Salt Lake City and Los Angeles tracer data, Atmos. Environ., 37, 5069–5082,, 2003. a

Hayati, A. N., Stoll, R., Pardyjak, E. R., Harman, T., and Kim, J.: Comparative metrics for computational approaches in non-uniform street-canyon flows, Build. Environ., 158, 16–27,, 2019. a

Hertwig, D., Soulhac, L., Fuka, V., Auerswald, T., Carpentieri, M., Hayden, P., Robins, A., Xie, Z.-T., and Coceal, O.: Evaluation of fast atmospheric dispersion models in a regular street network, Environ. Fluid Mech., 18, 1007–1044,, 2018. a, b

Horn, R. A. and Johnson, C. R.: Matrix Analysis, Cambridge University Press, 2 edn.,, 2012. a

Huang, C.: A theory of dispersion in turbulent shear flow, Atmos. Environ., 13, 453–463,, 1979. a

Kim, J., Moin, P., and Moser, R: Turbulence statistics in fully developed channel flow at low Reynolds number, J. of Fluid Mech., 177, 133–166,, 1987. a

Langevin, M. P.: Sur la théorie du mouvement brownien, Note de M. P. Langevin, présentée par M. Mascart, Comptes rendus hebdomadaires des séances de l’Académie des sciences, Série physique, Séance du 9 mars 1908, tome 146, 530–533, (last access: 26 September 2023), 1908. a

Legg, B. J.: Movement of plant pathogens in the crop canopy, Philos. T. Roy. Soc. Lon. B, 302, 559–574,, 1983. a

Mahaffee, W. F., Margairaz, F., Ulmer, L., Bailey, B. N., and Stoll, R.: Catching Spores: Linking Epidemiology, Pathogen Biology, and Physics to Ground-Based Airborne Inoculum Monitoring, Plant Disease, 107, 13–33,, 2023. a

Margairaz, F., Bozorgmehr, B., Gibbs, J., Singh, B., Willemsen, P., Pardyjak, E., and Stoll, R.: UtahEFD/QES-Public: v2.0.0, Zenodo [code],, 2022a. a, b, c

Margairaz, F., Eshagh, H., Hayati, A. N., Pardyjak, E. R., and Stoll, R.: Development and evaluation of an isolated-tree flow model for neutral-stability conditions, Urban Climate, 42, 101083,, 2022b. a, b

Maronga, B., Banzhaf, S., Burmeister, C., Esch, T., Forkel, R., Fröhlich, D., Fuka, V., Gehrke, K. F., Geletič, J., Giersch, S., Gronemeier, T., Groß, G., Heldens, W., Hellsten, A., Hoffmann, F., Inagaki, A., Kadasch, E., Kanani-Sühring, F., Ketelsen, K., Khan, B. A., Knigge, C., Knoop, H., Krč, P., Kurppa, M., Maamari, H., Matzarakis, A., Mauder, M., Pallasch, M., Pavlik, D., Pfafferott, J., Resler, J., Rissmann, S., Russo, E., Salim, M., Schrempf, M., Schwenkel, J., Seckmeyer, G., Schubert, S., Sühring, M., von Tils, R., Vollmer, L., Ward, S., Witha, B., Wurps, H., Zeidler, J., and Raasch, S.: Overview of the PALM model system 6.0, Geosci. Model Dev., 13, 1335–1372,, 2020. a

Miller, N. E., Stoll, R., Mahaffee, W. F., and Neill, T. M.: Heavy particle transport in a trellised agricultural canopy during non-row-aligned winds, Agr. Forest Meteorol., 256, 125–136,, 2018. a

Oettl, D.: Evaluation of the Revised Lagrangian Particle Model GRAL Against Wind-Tunnel and Field Observations in the Presence of Obstacles, Bound.-Lay. Meteorol., 155, 271–287,, 2015. a

Pardyjak, E., Speckart, S., Yin, F., and Veranth, J.: Near source deposition of vehicle generated fugitive dust on vegetation and buildings: Model development and theory, Atmos. Environ., 42, 6442–6452,, 2008. a, b

Pardyjak, E. R. and Brown, M.: QUIC-URB v1.1 Theory and User’s Guide, Los Alamos National Laboratory, (last access: 26 September 2023), 2003. a, b

Philips, D. A., Rossi, R., and Iaccarino, G.: Large-eddy simulation of passive scalar dispersion in an urban-like canopy, J. Fluid Mech., 723, 404–428,, 2013. a

Pirhalla, M., Heist, D., Perry, S., Tang, W., and Brouwer, L.: Simulations of dispersion through an irregular urban building array, Atmos. Environ., 258, 118500,, 2021. a

Pisso, I., Sollum, E., Grythe, H., Kristiansen, N. I., Cassiani, M., Eckhardt, S., Arnold, D., Morton, D., Thompson, R. L., Groot Zwaaftink, C. D., Evangeliou, N., Sodemann, H., Haimberger, L., Henne, S., Brunner, D., Burkhart, J. F., Fouilloux, A., Brioude, J., Philipp, A., Seibert, P., and Stohl, A.: The Lagrangian particle dispersion model FLEXPART version 10.4, Geosci. Model Dev., 12, 4955–4997,, 2019. a

Pope, S. B.: Consistency conditions for random-walk models of turbulent dispersion, Phys. Fluids, 30, 2374–2379,, 1987. a

Pope, S. B.: Turbulent Flows, Cambridge University Press, Cambridge,, 2000. a, b, c, d, e, f, g

Postma, J. V.: Timestep Buffering to Preserve the Well-Mixed Condition in Lagrangian Stochastic Simulations, Bound.-Lay. Meteorol., 156, 15–36,, 2015. a

Postma, J. V., Yee, E., and Wilson, J. D.: First-Order Inconsistencies Caused by Rogue Trajectories, Bound.-Lay. Meteorol., 144, 431–439,, 2012. a

Prussin, A. J., Marr, L. C., Schmale, D. G., Stoll, R., and Ross, S. D.: Experimental validation of a long-distance transport model for plant pathogens: Application to Fusarium graminearum, Agr. Forest Meteorol., 203, 118–130,, 2015. a

Ramamurthy, P., Pardyjak, E. R., and Klewicki, J. C.: Observations of the Effects of Atmospheric Stability on Turbulence Statistics Deep within an Urban Street Canyon, J. Appl. Meteorol. and Climatology, 46, 2074–2085,, 2007. a

Ramli, Huda Mohd. and Esler, J. G.: Quantitative evaluation of numerical integration schemes for Lagrangian particle dispersion models, Geosci. Model Dev., 9, 2441–2457,, 2016. a

Rew, R., Davis, G., Emmerson, S., Cormack, C., Caron, J., Pincus, R., Hartnett, E., Heimbigner, D., Appel, L., and Fisher, W.: Unidata NetCDF, UCAR [code and software],, 1989. a

Rodean, H. C.: The universal constant for the Lagrangian structure function, Phys. Fluids A-Fluid, 3, 1479–1480,, 1991. a, b

Rodean, H. C.: Stochastic Lagrangian Models of Turbulent Diffusion, American Meteorological Society, Boston, MA,, 1996. a, b, c, d

Roth, M.: Review of atmospheric turbulence over cities, Q. J. Roy. Meteor. Soc., 126, 941–990,, 2000. a, b

Röckle, R.: Bestimmung Der Strömungsverhältnisse Im Bereich Komplexer Bebauungsstrukturen, PhD thesis, Technischen Hochschule Darmstadt, (last access: 6 October 2023), 1990. a, b, c, d

Sasaki, Y.: Some Basic Formalisms in Numerical Variational Analysis, Mon. Weather Rev., 98, 875–883,<0875:sbfinv>;2, 1970. a

Seinfeld, J. H. and Pandis, S. N.: Atmospheric chemistry and physics: from air pollution to climate change, 3rd edn., Wiley, New York, ISBN 978-1-119-22117-3, 2016. a, b, c

Sherman, C. A.: A Mass-Consistent Model for Wind Fields over Complex Terrain, J. Appl. Meteorol., 17, 312–319,<0312:amcmfw>;2, 1978. a

Singh, B., Hansen, B. S., Brown, M. J., and Pardyjak, E. R.: Evaluation of the QUIC-URB fast response urban wind model for a cubical building array and wide building street canyon, Environ. Fluid Mech., 8, 281–312,, 2008. a, b, c, d, e, f, g, h, i, j

Singh, B., Pardyjak, E., Norgren, A., and Willemsen, P.: Accelerating urban fast response Lagrangian dispersion simulations using inexpensive graphics processor parallelism, Environ. Modell. Softw., 26, 739–750,, 2011. a

Smagorinsky, J.: General circulation experiments with primirive equations, Mon. Weather Rev., 91, 99–164,<0099:gcewtp>;2, 1963. a

Soulhac, L., Salizzoni, P., Cierco, F.-X., and Perkins, R.: The model SIRANE for atmospheric urban pollutant dispersion; part I, presentation of the model, Atmos. Environ., 45, 7379–7395,, 2011. a, b

Stiperski, I. and Calaf, M.: Generalizing Monin-Obukhov Similarity Theory (1954) for Complex Atmospheric Turbulence, Phys. Rev. Lett., 130, 124001,, 2023. a

Stockie, J. M.: The Mathematics of Atmospheric Dispersion Modeling, SIAM Rev., 53, 349–372,, 2011. a

Stoll, R. and Porté-Agel, F.: Dynamic subgrid-scale models for momentum and scalar fluxes in large-eddy simulations of neutrally stratified atmospheric boundary layers over heterogeneous terrain, Water Resour. Res., 42, 2719–2728,, 2006. a

Taylor, G. I.: Diffusion by Continuous Movements, P. Lond. Math. Soc., s2-20, 196–212,, 1921. a

Thiessen, L. D., Keune, J. A., Neill, T. M., Turechek, W. W., Grove, G. G., and Mahaffee, W. F.: Development of a grower‐conducted inoculum detection assay for management of grape powdery mildew, Plant Pathol., 65, 238–249,, 2016. a, b

Thomson, D. J.: Random walk modelling of diffusion in inhomogeneous turbulence, Q. J. Roy. Meteor. Soc., 110, 1107–1120,, 1984. a

Thomson, D. J.: Criteria for the selection of stochastic models of particle trajectories in turbulent flows, J. Fluid Mech., 180, 529–556,, 1987. a, b, c, d

Thomson, D. J., Physick, W. L., and Maryon, R. H.: Treatment of Interfaces in Random Walk Dispersion Models, J. Appl. Meteorol., 36, 1284–1295,<1284:toiirw>;2, 1997. a

Tinarelli, G., Mortarini, L., Castelli, S. T., Carlino, G., Moussafir, J., Olry, C., Armand, P., and Anfossi, D.: Review and Validation of MicroSpray, a Lagrangian Particle Model of Turbulent Dispersion. In Lagrangian Modeling of the Atmospherem, edited by: Lin, J., Brunner, D., Gerbig, C., Stohl, A., Luhar, A., and Webley, P.,, 2012.  a

Ulmer, L., Margairaz, F., Bailey, B. N., Mahaffee, W. F., Pardyjak, E. R., and Stoll, R.: A fast-response, wind angle-sensitive model for predicting mean winds in row-organized canopies, Agr. Forest Meteorol., 329, 109273,, 2023. a

US EPA Office Of Research And Development: CMAQ, Zenodo [software],, 2020. a

Vachat, R. D.: Realizability inequalities in turbulent flows, Phys. Fluids, 20, 551–556,, 1977. a

Williams, M. D., Brown, M. J., and Pardyjak., E. R.: Development of a dispersion model for flow around buildings, in: Fourth Symposium on the Urban Environment, Norfolk, VA, 19–24 May 2002, (last access: 26 September 2023), 2002. a, b

Wilson, J. D.: “Rogue Velocities” in a Lagrangian Stochastic Model for Idealized Inhomogeneous Turbulence. In Lagrangian Modeling of the Atmosphere, edited by: Lin, J., Brunner, D., Gerbig, C., Stohl, A., Luhar, A., and Webley, P.,, 2012. a

Yee, E. and Wilson, J. D.: Instability in Lagrangian stochastic trajectory models, and a method for its cure, Bound.-Lay. Meteorol., 122, 243–261,, 2007. a, b, c, d, e, f, g

Short summary
The Quick Environmental Simulation (QES) tool is a low-computational-cost fast-response framework. It provides high-resolution wind and concentration information to study complex problems, such as spore or smoke transport, urban pollution, and air quality. This paper presents the particle dispersion model and its validation against analytical solutions and wind-tunnel data for a mock-urban setting. In all cases, the model provides accurate results with competitive computational performance.