Interactive comment on “ Stochastic fields method for sub-grid scale emission heterogeneity in mesoscale atmospheric dispersion models ”

Chemistry transport models have to deal with simplified descriptions. The neglect of small-scale emission fluctuations (e.g. caused by highways, localized industry) is one of these simplifications. This paper investigates ways to account for this sub-gridscale variability. The specific solution proposed is to do an ensemble simulation in which different concentration fields are created. Each concentration field has a different emission, and as such the ensemble of fields represents the concentration fluctuations that are expected from sub-grid-scale emission variability.


Introduction
As pointed out by Galmarini et al. (2008, emission inventories used in atmospheric modeling may provide information on the emissions that is too detailed to be useful for the atmospheric models themselves.This is partly due to the fact that emission may be described at spatial Correspondence to: M. Cassiani (mc@nilu.no)scales finer than the grid resolution and also due to the incapacity of air quality models to resolve processes below the computational grid size.The wealth of information that an emission inventory can potentially provide is completely lost in the volume averaging process usually performed within transport models and no information on the sub-grid emission heterogeneity is therefore transferred to the atmospheric layers.All this results in an unspecified emission heterogeneity that is lost forever in the advection-diffusion process.In an attempt to overcome this paradox, G2008 proposed a second order closure sub-grid formulation applicable to Reynolds Averaged Navier Stokes (RANS) models working at any scale from the meso to the global one.The model provided a way to translate sub-grid scale emission heterogeneity into concentration fluctuations through turbulent transport theory.Large-eddy simulations (LES) of the dispersion from simple prototype sources were used for evaluation of their closure and to study and obtain the necessary closure constants.The results demonstrated the nature of the problem and that it can significantly influence the sub-grid scale concentration fluctuations.Furthermore the sub-grid model showed that by accounting for the second moment, a substantial amount of information can be recovered.
This work develops in the footsteps of G2008 addressing the same problem but suggesting an alternative solution.As in G2008, emission heterogeneity below the grid size of a RANS model is defined as the deviation in terms of mass from the volume average normally used by RANS, due to the presence of individual and different emitting elements smaller than the grid size.We propose an alternative and in some aspects more adequate way of expressing emission heterogeneity and its impact on the variability of the air concentration based on the stochastic fields method (Valiño, 1998).
Published by Copernicus Publications on behalf of the European Geosciences Union.
M. Cassiani et al.: Stochastic fields method for sub-grid scale emission heterogeneity In Sect.2, the new approach will be described in conceptual terms prior to the formal description provided in Sect.3. Therein the necessary concepts of the stochastic fields method will be introduced and the problem of including sub-grid scale source heterogeneity within this framework will be appropriately formulated.In Sect. 4 we will follow the work of G2008 and evaluate the results obtained by the stochastic fields method using high resolution largeeddy simulation and try to reproduce the case study presented in G2008 that will be treated as benchmark.Finally conclusions and future plans will be outlined in Sect. 5.

The approach
We start by proposing a representation of emissions within a model grid cell by means of an emission probability density function (PDF).The better adequacy of such a representation with respect to the classical one lies in the fact that although RANS models are often referred to as deterministic models they effectively work on the basis of statistical definitions that are limited to the representation of first moments (e.g.mean concentration).The representation of emission inventories in the form of PDFs and the possibility to translate this information into air concentration seems, therefore, an appropriate choice.It fits better with the statistical nature of the turbulent dispersion and provides a complete representation of the emission and concentration range from the statistical viewpoint since the knowledge of PDFs allows recovery of any other moment of the distribution.
The second step of the approach proposed therefore consists of finding a way to translate the emission PDF in atmospheric concentration.A full set of modeling methods, known as PDF methods for turbulent reacting flows (e.g.Pope 1985Pope , 2000;;Dopazo et al., 1997;Valiño, 1998;Fox, 2003;Heinz, 2003;Cassiani et al., 2005aCassiani et al., , 2007a)), exists that provides a way to compute the evolution of the concentration PDF and all its moments within a turbulent flow.In this context it is natural, therefore, to consider the sub-grid emission heterogeneity as represented by a PDF and to use these methods to convert it to an average concentration and all its statistical moments.Some applications of PDF methods to atmospheric dispersion and chemistry have recently appeared in the literature using both the Lagrangian particle approach (Cassiani et al., 2005a(Cassiani et al., , b, c, 2007a, b;, b;Luhar and Sawford, 2005;Sawford, 2006;Dixon and Tomlin, 2007;Bakosi et al., 2009) and the stochastic fields method approach (Garmory et al., 2006(Garmory et al., , 2008)).However, all of these works focused on a single source or few sources treated in great detail to improve our understanding of the turbulent dispersion processes and turbulence chemistry interaction.Here we will follow the stochastic fields method (Valiño, 1998) and we will demonstrate how this can be used as a straightforward extension of commonly used mesoscale chemical transport models (CTM) to include sub-grid emission heterogeneity in dispersion calculations.
Although our simulations and results will be limited to non-reactive substances the inclusion of any chemical reaction in the context of the PDF methods is straightforward and for this property PDF methods are commonly used in engineering fluid mechanics and turbulent combustion (e.g.Pope, 1985Pope, , 2000;;Dopazo et al., 1997;Valiño, 1998;Fox, 2003;Heinz, 2003).The PDF method does not require any closure for the chemistry thus solving the problems connected with the relative speed of chemical reaction and turbulent mixing that would arise in the case of fast chemistry and the use of first or second order closures (e.g.Vila-Guerau de Arellano et al., 2004).

Stochastic fields method
The stochastic fields method (SFM) is a technique for solving the joint scalar PDF transport equation for a turbulent flow.Loosely speaking a stochastic field should be understood as a spatial scalar field member of an ensemble of realizations sharing one-point statistical properties with an ensemble of true realizations of the scalar field dispersing in a turbulent flow.It was originally proposed by Valiño (1998) and first applied to atmospheric dispersion of reactive scalars by Garmory et al. (2006Garmory et al. ( , 2008)).The mathematical and physical foundation of this novel approach has also been investigated and further developed by Sabelnikov and Soulard (2005), who in particular relaxed some of the original assumptions made by Valiño (1998) related to the properties of the stochastic field.Here we will give a brief overview of this method and the first step is to present the transport equation for the one-point joint PDF f φ ≡ f φ (ψ ;x,t), of a set of dynamically passive chemical components with concentration where R is the source term associated with chemical reactions, S any other source term (e.g.emissions), U is the mean velocity vector and is the molecular diffusivity.u i |φ = ψ represents the expected value of the fluctuating velocity conditional on the set of scalar φ taking some specific value ψ. ψ is formally referenced as the sample space variable (e.g.Pope, 2000).This equation and the methods to derive it from the advection diffusion equation of reactive scalars are extensively discussed in the literature (e.g.Pope, 1985Pope, , 2000;;Dopazo et al., 1997;Heinz, 2003;Fox, 2003;Atmos. Chem. Phys., 10, 267-277, 2010 www.atmos-chem-phys.net/10/267/2010/Hauke and Valiño, 2004).Some introductory concepts are reviewed here.First we recall that f φ (ψ ;x,t) dψ is the probability of obtaining simultaneously the scalar concentrations φ at time t and point x between ψ α < φ α (x,t) < ψ α + dψ α .The terms on the left hand side of (1) are closed and do not require any modeling assumption given the knowledge of an averaged velocity field.They are, from left to right, the temporal variation of the PDF, the advection by the mean (averaged) velocity field, the source term associated with chemical reactions, any other source term, and finally the last term on the left represent the diffusion of the PDF in physical space by molecular diffusivity.This is usually neglected for high Reynolds number flow.Although we are not considering reactive scalars it is worth remembering that whatever the complexity (nonlinearity) of the chemical reaction term it is described here in a closed form, meaning that no modeling assumption is necessary since it is a local term not involving multi-point information (e.g.Pope, 1985Pope, , 2000;;Fox, 2003).The first term on the right is the conditional turbulent convection term and represents the convection by the unresolved fluctuating velocity.In Eq. ( 1) this term appears in its unclosed form since the joint scalar PDF does not carry any information about the velocity field.This would be a closed term if the joint velocity-scalar PDF transport equation was to be considered (e.g.Pope 2000;Fox, 2003;Cassiani et al., 2005a).The second term on the right is the conditional scalar dissipation term.It represents a transport in the scalar space, and it is responsible for the dissipation of the scalar fluctuations.Since the one-point joint scalar PDF does not have any information on the spatial derivative of the scalar this is an unclosed term and needs some modeling assumption for its closure.
A possible closed form of the scalar PDF transport equation is (e.g.Pope, 1985Pope, , 2000;;Dopazo et al., 1997;Valiño, 1998;Fox, 2003), where the conditional turbulent convection is closed using a simple gradient diffusion scalar hypothesis with diffusion coefficient K and the conditional scalar dissipation term is closed using the simple interaction by exchange with the mean model (IEM), Villermaux and Devillon (1972), also called linear mean square estimation model (LMSE) by Dopazo and O'Brien (1974).This model is based on a linear relaxation of the local concentration towards the local average value and the relaxation rate is controlled by the time scale T mix .A concise overview of the properties of this mixing model is given in Cassiani et al. (2005a), other more physically consistent and complex models exist (e.g.Dopazo et al., 1997, Fox, 2003;Cassiani et al., 2005a) among which the IECM (interaction by exchange with the conditional mean).However, the implementation of this last model would require the modification of the meteorological model generating the wind field and was not considered here.The most straightforward definition of the time scale is T mix = C φ T , T ≡ e/ε where e is the mean turbulent kinetic energy, ε is the mean turbulent energy dissipation rate and T is therefore the characteristic turbulent time scale.The constant of proportionality should be of the order of one but some variability is observed and the time scale may depend on several factors, see Cassiani et al. (2005a) for a short review.Valiño (1998) firstly demonstrated that the solution of the above Eq.( 2) on a discrete grid is equivalent to the solution of an ensemble of (N ) stochastic partial differential equations (SPDE), where n indicate the members of the ensemble (1 ≤ n ≤ N ) and φ n (x,t) must be interpreted here as a realization of a smooth stochastic field and not as a true scalar field dispersing in a turbulent flow.As remarked by Valiño (1998), the actual and the stochastic scalar fields are statistically equivalent in the sense that they contain the same one-point statistical information, i.e. their one-point PDF is the same and satisfies Eq. ( 2).This also means that no attempt should be made to extract multi-point information from the ensemble of stochastic fields (e.g.two-point scalar correlations).From the ensemble of stochastic fields any one-point moment of the scalar field can be extracted, for example: consideration of different diffusion coefficients depending on anisotropic meshes.However, for our purpose in the following discussion a simplified version of Eq. ( 3) will be introduced; a single dynamically passive and non-reactive scalar will be considered and the assumption of neglecting any effect of turbulent transport in the horizontal directions will be made.This latter approximation is very often made in large or mesoscale CTM.The resulting simplified equation is, In the numerical solution of Eq. ( 6) we will closely follow the suggestion of Garmory et al. (2005Garmory et al. ( , 2008) ) and make use of the fractional step method.More details are given in Sect.4.1, for now it suffices to note that according to this approach a standard solver for the Reynolds averaged advection diffusion equation can be used for the advection diffusion part of Eq. ( 6), while a standard Monte Carlo PDF technique can be used for the stochastic contribution to the equation.In this sense the stochastic fields method can be implemented as a straightforward extension of a standard CTM.

Sub-grid scale source heterogeneity
First it is important to note that in Eq. ( 6) the source term S n can be stochastic in nature, i.e. variable among different stochastic fields.This means that consistently with the PDF transport equation used to obtain Eq. ( 6) we are actually able to naturally include a PDF of the source strength in this equation.Second it should be noted that the stochastic fields method has also been formulated in terms of volume averages (Mustata et al., 2005) obtaining equations formally equivalent to Eqs. ( 2) and (3) where the only difference is that the averages must be interpreted as the resolved (at the grid level) velocity and scalar fields.For our purpose this means that the PDF obtained from the ensemble of stochastic fields can be interpreted as generated by the resolved and sub-grid scale features of the scalar field.
As pointed out by G2008, mesoscale and other large scale air quality models have as available input averaged emission rates for various primary pollutants, which accounts for the volume averaged quantity of mass released per unit time.No other information takes into account the fact that for example a large amount of mass can be emitted by a small portion of the grid surface or by several sources scattered around it.G2008 refer to this as sub-grid emission heterogeneity.The emission heterogeneity can be seen by disaggregating an emission inventory for a specific species over a mesh more refined than the one used in the atmospheric transport model calculation and we will call this a sub-grid emission mesh.Within each element of the sub-grid emission mesh different surfaces will emit potentially different amounts of mass.
Therefore, the problem is formulated as follows: given a dynamically passive pollutant, which is emitted by a series of surfaces within the model grid with different rates, is there any way to transfer the information on this sub-grid emission heterogeneity to the atmospheric concentration of the species?The answer using the stochastic field method approach is that this must be done according to a sub-grid scale surface emission PDF from which the source strength used in any stochastic field calculation, S n in Eq. ( 6), will be extracted.
When considering a rectangular computational mesh element of size x y z the sub-grid emission mesh element can be represented as δ x δ y δ z with δ x = x /N x , δ y = y /N y , δ z = z /N z , where N x,y,z represents the refinement of the emission grid with respect to the computational mesh in each direction.With increasing N x,y,z the emission per unit volume will converge towards the true emission.This general definition allows for emissions not necessarily at the ground surface.However, our present scope is more limited and the assumption of a surface emission vertically diluted in the first computational cell is made, δ z = z .This is not a very restrictive hypothesis since usually the vertical extension of the first ground element of a mesoscale CTM is only of few meters.These concepts were originally introduced in G2008, where it was noticed that defining the emission for a single i − th sub-grid mesh emission surface as where M i is the amount of mass emitted per unit time.An average emission value over the surface of the computational element can be defined as and similarly for higher moments These arguments can be straightforwardly extended and bring to a sub-grid surface emission PDF.We note that the relationship between the surface emission per unit area E and the emission S in Eq. ( 6) is simplyS = E/ z .

Comparison of the stochastic fields method with large eddy simulations
The stochastic fields method (SFM) formulation proposed in the previous section will be tested by means of large-eddy simulation (LES) of controlled emission cases and a three dimensional RANS model is modified to include the Stochastic fields method termed RANS-SFM model.For the sake of simplicity and comparability with previous work we chose the same evaluation strategy used by G2008.Following G2008 the LES model was used as it guarantees a detailed representation of the turbulent flow and dispersion in the atmospheric boundary layer ranging from the peak of maximum spectral energy down to the inertial subrange scales.In such controlled flow conditions it is possible to define specific, detailed and controlled emission scenarios which allow for rigorous evaluation of the approach.The large-eddy simulations are run on a domain size corresponding to a few grid cells of a mesoscale model.Specifically the atmospheric flow contained in a volume of 12 km×12 km×1500 m is simulated by LES using 120×120×60 grid points whereas the RANS-SFM model will use 4×4×60 cells in total, leading to a horizontal grid size of 3 km×3 km (see Fig. 1).Slab averaging operations are performed at all vertical levels to make the LES results comparable to the RANS-SFM.The slabs over which the averages are taken correspond to the dimensions of the individual RANS-SFM model cells.In Fig. 1, six of the sixteen cells of the RANS-SFM model (corresponding to 16 LES sub-domains) have been labeled A through F in order to facilitate the analysis of the results.The 6 cells are representative in terms of source location and flow direction.

Models and simulation details
The LES model is the one used by G2008 and originally developed by Cuijpers and Duynkerke (1993), Siebesma and Cuijpers (1995), Cuijpers and Holtslag (1998) and Vilà-Guerau de Arellano and Cuijpers (2000).The model has evolved over the years and has been successfully used to study many different processes in the atmospheric boundary layer.A full list of references and the related applications can be found in G2008 and the reader is referred to this paper and the reference therein for all the details.The cases analysed here are the same as those of G2008 and details can be found in that paper, but for the sake of clarity we will summarize briefly the case specifications.The LES simulation runs for 3 h with a maximum time step used in the calculations of 0. Following G2008 a RANS model has been selected to host the SFM approach.The three-dimensional RANS model is the mesoscale model described in detail by Martilli (2002).Periodic boundary conditions are used for the dynamics while non-periodic boundary conditions are used for the tracer.Initial conditions for the dynamics are obtained from the LES.The average concentration equation has been modified to allow for the application of the SFM approach for the transport of passive scalars.The numerical solution of Eq. ( 6) is carried out following the suggestion of Garmory et al. (2005Garmory et al. ( , 2008) ) who apply the fractional step method.The first step involves the solution of a standard advection diffusion equation, for which any solver usually available in CTMs can be used.
In the present case the piecewise parabolic method with the correction algorithm explained in Clappier (1998) is used for the advection and an implicit scheme for the diffusion.We remark that this equation must be solved in parallel for all of the N stochastic fields.The N fields resulting at the end of this first fractional step constitute the starting point for the subsequent calculations based on a Monte Carlo or stochastic approach (e.g.Valiño, 1998;Garmory et al., 2005Garmory et al., , 2008) ) and involve therefore the solution of the following stochastic differential equation in any point of the computational grid: As customary in Monte Carlo PDF methods the solution of this equation is further split in different steps (see e.g.Fox, 2003;Garmory et al., 2008).First the source and stochastic term are considered, and afterwards the mixing term is advanced.It is underlined that the mixing step involves the use of the average φ obtained from Eq. (4).We note, as also mentioned by Fox (2003), that if chemical reaction source/sink terms exist they should be advanced separately as the last fractional step.As discussed by Garmory et al. (2006Garmory et al. ( , 2008) ) the numerical simulations of the random Wiener term (dW ) introduces a random Gaussian jump at each time step and it would therefore be possible to simulate negative concentrations/fields.This can in principle be avoided using extremely small time steps, but it comes at the expense of extremely large computation times.We note that in a continuous analytical sense this issue does not exist since the gradient goes to zero when the concentration goes to zero (e.g.Valiño, 1998).Here, to avoid this issue, the method suggested and discussed in Garmory et al. (2006Garmory et al. ( , 2008) ) is used; basically the random increment is limited to values not exceeding the magnitude of the local concentration.
The mixing time scale is defined as T mix = T ≡ e/ε, thus implicitly assuming the value of unity for the proportionality constant without attempting any kind of adjustment.The turbulent time scale is provided by the mesoscale model where it is defined as, T = l ε /(e 1/2 C ε ) and it is therefore a function of the elevation above ground.Here l ε is a length scale from Bougeault and Lacarrère (1989) and C ε (= 0.125) is a proportionality constant.The time step used to discretise Eq. ( 11) is chosen to be min(T mix )×10 −2 , i.e. a small fraction of the smaller mixing time scale in the domain.This is smaller than the global time scale used to update Eq. ( 10) and the two equations are synchronized at every global time step.
The source emission strength S n is defined according to the sub-grid scale (SGS) emission PDF that must be known a-priori.Following G2008 the sub-grid scale emission heterogeneity is generated considering a single area source emitting within one of the RANS-SFM cells that varies in size from a total coverage of the cell to, possibly, the minimum resolvable size for the LES.Four different emission scenarios were actually simulated.For the first the passive tracer is released over the entire RANS-SFM grid element.The remaining three cover the grid elements corresponding to 64, 44 and 28% of the grid surface thus leading to an increasing sub-grid scale variability of the emission and corresponding emission variance (see Fig. 1).The uniform emitting surface releases continuously with an emission rate of 0.1 ppb m s −1 .In the LES simulations the overall flux is maintained constant in spite of the change in the size of the emitting surface.Therefore the source strength is defined as S = 0.1/(A z)ppbs −1 where A(= 1, 0.64, 0.44, 0.28) represents the surface fraction of the RANS-SFM grid cell covered by the source.In the RANS-SFM context this must be interpreted as follows: i) for the uniform source we have only one possible value therefore for any of the N stochastic fields the source is deterministic, S n = 0.1/( z)ppbs −1 ii) for the non-uniform sources at the sub-grid level two values are possible, S = 0 with a probability corresponding to the surface fraction of the non emitting surface, and S = 0.1/(A z)ppbs −1 with a probability corresponding to the surface fraction of the emitting surface.This means that a fraction of the N stochastic fields equal to N×(1 − A), have no emission while the remaining fraction, equal to N×A, emit S = 0.1/(A z)ppbs −1 .This can be rewritten more compactly as, Formally this source function could be rewritten as a discrete two value density function S(ψ) = δ(ψ)(1 − A) + δ(ψ − ψ 0 )A with ψ 0 = 0.1/(A z) where δis the Dirac's delta function.

Comparison of mean concentrations
Figure 2 shows the mean concentration predicted by the LES and the RANS-SFM with 3 km resolution.The RANS-SFM results have no dependence on the SGS emissions PDF and the data reported in Fig. 2 are representative for all the fraction of surface source coverage.The LES results are weakly dependent on the source surface coverage and some noticeable differences arise between the different source coverage's in sub-domain A and B, while no significant difference arise in sub-domain C and D, where the value reported is the one corresponding to 100% source surface coverage and it is representative for all the source configurations.Sub-domains A and B are only weakly influenced by mean advection and present very low values for the mean concentration.The LES simulations with 100% source surface coverage generate the highest value of mean concentration since the boundary of the source (sub-domain C) is in direct proximity of the boundary of A and B. For the smaller area sources the mean concentration decreases due to the increasing distance between sub-domains A and B and the source boundary (see Fig. 1), consistent with the decrease in the surface coverage.In sub-domain B the mean concentration value for the 64%, 44% and 28% sources are reported (thin dashed lines) together with the 100% coverage (thick dashed line).The 28% source has very low mean concentration in this sub-domain.In sub-domain A only the 100% source coverage has an average concentration significantly different from zero, while the smaller source can be barely observed on a linear scale.
In general, and as expected following the results obtained by G2008 using second order RANS simulations, the mean concentration predicted by the RANS-SFM model and the LES are always in very good agreement in sub-domains C and D, corresponding to the source and downwind cell respectively, and show some departure in sub-domain A and B when the surface of the source shrink according to the resolution limits of the RANS-SFM model.This confirms that the application of the stochastic terms to the RANS does not alter the level of agreement obtained by the original model and that the numerical treatment of the stochastic term as proposed by Garmory et al. (2008) is satisfactory.The source of the small discrepancies between the RANS-SFM model and the LES model with 100% source coverage can be related to small differences in the wind velocity field reproduced by the LES and the RANS model (G2008).We remark again that the RANS-SFM model, as any RANS model, neglects the spatial localization at the sub-grid scale and cannot capture differences in the mean field arising from SGS source heterogeneity.As previously shown by G2008 no significant concentration values are predicted within the sub-domains E and F by both the LES and the RANS model and they will not be referred to in the following discussions.This is due to the progressive rotation of the mean wind vector towards the left of the assumed higher elevation westerly wind while progressively moving closer to the surface.

Comparison of concentration standard deviation
Figure 3 shows the comparison for the standard deviation of concentration for the three sources with partial surface coverage ranging from 63% (top panels), 44% (middle panels) to 28% (bottom panels).Focusing on the sub-domain C, which contains the emitting source, an overall good agree- ment is observed over the whole boundary layer.However some noticeable differences emerge between the LES and the RANS-SFM which shows larger fluctuations than the LES above z/ h ≈ 0.2 and lower fluctuations below this elevation.The magnitudes of the discrepancies at the first grid level are below 20% with higher values for the source with smaller surface cover.The discrepancies can be higher above z/ h ≈ 0.2 and reach a maximum at z/ h ≈ 0.8.Looking jointly to Fig. 2 and Fig. 3 the higher standard deviation value above z/ h ≈ 0.2 generated by the RANS-SFM can be related to higher concentrations of the tracer, while the lower value below z/ h ≈ 0.2 must be attributed to a faster dissipation of concentration fluctuations with respect to the LES.This last discrepancy could be corrected using a different constant of proportionality between T mix and T , nonetheless the level of agreement is already fully satisfactory.The situation in the downwind sub-domain D, is different.Here the effect of the advecting wind field is more significant and, although the agreement between the standard deviation profiles in the bulk of the boundary layer is acceptable, the values close to the surface are under predicted.This behavior is similar to the one observed in G2008 for their second order closure RANS.
In Fig. 4a comparison between the results of the second order model of G2008 and the present RANS-SFM model is shown for the intermediate surface cover (44%) which is indicative for all the emitting surfaces.Being subject to identical wind field and numerical solver this last comparison further validates the finding of the second order closure model of G2008.Going back to Fig. 3 and the comparison between LES and RANS-SFM simulations in the sub-domain A and B the value of the standard deviations are very low reflecting the low value of the mean concentration, and therefore the amount of tracer advected in these cells.However, differences between the stochastic model and the LES emerge when the effect of the source coverage on the standard deviation is considered.As explained previously, smaller surface coverage in the LES means that its spatial distance from the boundaries of the sub-domain A and B increase.This is reflected by lower values of mean concentration in these cells, as discussed above and shown in Fig. 2. Correspondingly the value of standard deviation also decreases when the surface coverage is reduced from 63% to 44% and 28%.This behavior cannot be captured by the RANS-SFM where the spatial resolution is 3 km irrespective of the source configuration and the sub-grid scale emission variability has no spatial localization but only a probabilistic representation.Therefore, the decrease in source coverage affects only the concentration fluctuations with the standard deviation increasing in both A and consistent with the increase observed in the sub-domain C. The second order RANS model of G2008 has a similar behavior consistent with the lack of sub-grid scale spatial localization.
Comparing Fig. 3 and Fig. 2 (sub-domains A and B) the overall picture shows that the value of concentration intensity (σ φ / φ ) for the LES are actually higher than those for the RANS-SFM model consistent with the fact that the sub-domains B and A are at larger crosswind distance from the emitting source (see e.g.Sawford, 2004) when the source surface shrinks.However, any discrepancies in subdomains A and B pertain to values of mean concentration and standard deviation that are more than one order of magnitude smaller than that observed in the emitting cell.

The effect of sub-grid scale source probability distribution function (SGS source PDF) choice
Any calculations up do this point have been performed using the SGS emission PDF reported in Eq. ( 12), i.e. a two-value distribution.However, it is of interest to investigate the influence of the particular functional form for the PDF on second order and higher order statistics.Most of the time in practical applications the exact form of the SGS source PDF will be unknown and an approximate form will be used to extract the realizations.Here we test the influence of the emitting PDF comparing the results obtained using Eq. ( 12) with a clipped Gaussian PDF (Lewellen and Sykes, 1986) where the value randomly extracted from a Gaussian PDF are assigned a value of zero when negative values occurs.The mean and variance of the Gaussian PDF are appropriately chosen such that the overall PDF, including the occurrence of zero values, has the mean and variance corresponding to the values from the PDF in Eq. ( 12).This type of PDF is only one of many possible PDFs and other forms may also be applied, see e.g. the short review by Yee et al. (1993).Figure 5, reports the values obtained for the standard deviation, the cubic root of the third centered moment, (φ − φ ) 3 1/3 , and the skewness Sk= (φ − φ ) 3 σ −3 for the two different SGS source PDFs and the LES for the intermediate case with 44% source surface coverage.The standard deviation is almost unaffected by the change in the SGS source PDF with only minor changes in sub-domain C and hardly any difference in the other sub-domains.The third centered moment is more affected by the PDF choice and particularly in sub-domain C some noticeable difference between the two PDFs arise.The two-value PDF agrees slightly better with the LES results.Note that for the third centered moment the agreement between the LES and RANS-SFM results in sub-domain A and B indicates a much higher level of fluctuation in the LES, given that the mean concentration generated by the LES for the 44% source was shown to be much lower in these sub-domains (Fig. 2).As expected the skewness amplifies the differences between RANS-SFM and LES and it is much more sensitive to the change of the SGS PDF.However, the differences in the skewness profiles between the two PDF sources are not dramatic with values within 50% of each other in all cases.

Skewness predictions RANS-SFM versus LES
The comparison of the skewness predicted by the RANS-SFM and the skewness predicted by the LES needs further comment.The skewness is a measure of the asymmetry of the concentration PDF and, being a normalized moment, it is extremely sensitive to the spatial position within a dispersing plume with extremely high values and spatial variation rates near the plume edges, where intermittency in the concentration is large (e.g.Sawford, 2004;Luhar and Sawford, 2005).
The RANS-SFM model shows only a partial agreement with the LES results and this is due to the fact that the LES dispersing field and the RANS-SFM dispersing field are actually at slightly different horizontal and vertical positions, i.e. there is a small misalignment of the two dispersing fields due to the small differences in the velocity field and the lack of spatial localizations at the SGS in the RANS-SFM model, as already outlined in G2008.These small differences are relatively unimportant close to the location of higher mean concentration, (e.g. in any downwind position close to the ground), but become more significant in the more peripheral dispersing locations.However, since the concentration in these peripheral dispersing locations are very low these differences are not obvious until the moments are normalized (e.g. by the mean or the standard deviation).
The disagreement in sub-domain A and B is related to the larger distance of subdomain A and B from the source in the LES with respect to the RANS-SFM and this is a consequence of the lack of spatial localization of the SGS.The discrepancies in sub-domain C are limited to higher elevation (z/h>0.25)and due to small differences in the vertical rotation rate of the mean velocity in the RANS-SFM and LES.The comparison in sub-domain D is less affected by any of these issues since the tracer is spread over the entire vertical extension of the atmospheric boundary layer and it is approximately located downwind from the source, i.e. there is no peripheral dispersion zone included in sub-domain D. The LES and RANS-SFM models predict in this sub-domain quite similar values of skewness ranging between 2 and 4. Luhar and Sawford (2005), using a Lagrangian PDF model based on the joint velocity concentration PDF, obtained a similar value of 3 throughout the vertical extension of the Boundary layer, at a comparable downwind distance from an extended line source.
It must be remarked that in the position of higher mean concentration (i.e.sub-domain C and D close to the ground) where the population is more exposed to pollutants the skewness predicted by the RANS-SFM and the LES are in fair agreement despite the small differences in the mean field and the lack of spatial localization of the SGS.This is a remarkable result and means that the RANS-SFM is able to predict with a satisfactory accuracy the tails of the concentration distribution and therefore the range of possible concentration to which the population may be exposed.

Effect of increased RANS resolution
Figure 6 reports the results of the RANS-SFM model with an increase in the spatial resolution from 3 km to 1 km while conserving the same source size and two-value SGS emitting source PDF.The increase of the resolution has some influence on the predicted mean concentration field (top panels) that now seems to agree even better with the results for the LES at 100% source surface coverage in sub-domains A, B and C. In sub-domain B the thin dashed line refers to the 44% source while the 100% source surface coverage is the same as the RANS-SFM predictions.We note that with respect to the mean concentration field the results of the LES with 100% source surface coverage are the most accurate results that can be achieved with the RANS-SFM since no spatial localization is considered at the SGS, see also the discussion relating to Fig. 2 in Sect.4.2.
The standard deviation for the 44% SGS source is shown in the middle panels and demonstrates that with the increased resolution a better match is possible in both the subdomain C, and the downwind sub-domain D. No great improvement is observed in the sub-domain A and B, although now the ratio of the standard deviation and the mean in B seems to be closer to the value predicted by the LES.The skewness is reported in the bottom panels for 44% SGS source, and some overall improved agreement is observed for the sub-domain B, with values now comparable between the LES and the RANS-SFM.An improvement is also visible in sub-domain C where the values are now similar up to an elevation of about z/ h ≈ 0.7.The agreement remains good for the sub-domain D.

Conclusions
The stochastic fields method has been applied to the problem of sub-grid emission heterogeneity in mesoscale and largescale dispersion models.Under this general formalism a transport equation for the concentration PDF is solved thus allowing a straightforward and consistent inclusion of subgrid emission heterogeneity in a probabilistic way.
The approach is rather general in nature and accounts not only for the sub-grid emission heterogeneity but also for any sub-grid concentration variability generated by the turbulent mixing processes (e.g.Cassiani et al., 2005;Garmory et al., 2006).Moreover, a property of this approach, well known in the turbulent combustion community where this method has already been applied, is that any chemical reaction is included in a formally closed way, meaning that no additional closure issues arise due to non-linear chemical reactions of any complexity.This of course makes this approach particularly suitable for chemical transport models where second order chemical reactions are often of great importance and pose a significant closure challenge (e.g.Vilà-Guerau de Arellano et al., 2004).Future development of the work will involve the study of chemical reactive substances and real world simulations to evaluate the influence of sub-grid emission heterogeneity and the related concentration fluctuations on atmospheric chemistry.
The method is straightforward to implement in any chemical transport model.In the present case an ensemble of N = 100 realizations was used but a smaller ensemble of N = 50 realizations gave almost indistinguishable results, further decreasing to N = 25 produced noticeably different results.Finally we note that when non-reacting scalars are involved then the mean concentration needed in Eq. ( 11) could be pre-calculated using a standard RANS model, thus allowing the use of the stochastic fields method with serial simulations instead of parallel ones.This simplification is not possible for reactive scalars which require parallel realizations to include chemical reactions in closed form (e.g.Fox, 2003;Cassiani et al., 2005;Garmory et al., 2008) and this was the method followed here.The RANS-SFM is a stochastic method and therefore the inclusion of more chemical reactive tracers increases the computational burden linearly (e.g.Dopazo et al., 1997;Valiño, 1998;Pope, 2000;Fox, 2003;Sabel'nikov and Soulard, 2005).It is not expected that emission variability in adjacent grid cells will increase significantly the computational requirement.However, this is a newly proposed method and extensive test will be carried on in the future to correctly evaluate the model performance under different emission scenarios.
The proposed approach makes available to chemical transport model users a full set of high order statistics that are not currently considered.These statics could be used, for example, to define sub-grid ranges of concentration values to which gridded populations are most likely exposed instead of using a single average value.The consequence of this for future exposure studies is yet to be addressed.
5 s.The surface sensible heat flux is set at +0.05 K m s −1 .A constant westerly wind of 5 ms −1 has been imposed at higher elevation and Coriolis terms are included in the simulation.The initial potential temperature profile has a constant value of 288 K below 662.5 m and increases by 6 K above 712.5 m.The surface roughness length z 0 is set to 0.01 m.These conditions generated what can be considered a convective boundary layer.Periodic lateral boundary conditions are assumed.At the end of the first hour temperature and wind profiles are provided to the RANS-SFM model as initial condition for its simulation.LES data are averaged over the last simulation hour before they are compared with the RANS-SFM results.Non-periodic boundary conditions are used for the tracer.The simulations run for 3 h for the dynamics and 2 h for the tracer dispersion.All the results relate to the last hour of the simulation.

Figure 1 .Fig. 1 .
Figure 1.Simulation domain and source configuration for the RANS-SFM and the LES models.The source is located in sub-domain C. The dominant wind direction is from left to right.

Figure 2 .Fig. 2 .
Figure 2. Mean concentration field.Thick dashed line: LES with 100% source surface coverage.Thick continuous line: RANS-SFM model irrespective of the source surface coverage.In A and B the thin dashed lines are LES results for 64%, 44% and 28% source surface coverage respectively.Note the scale difference for the concentration axis between sub-domains A and B, and C and D.

Figure 4 .Fig. 4 .
Figure 4. Standard deviation of concentration for the RANS-SFM (continuous lines) and the second order RANS of Galmarini et al.(2008) (long dashed lines) .44% source surface coverage of sub-domain C.

Figure 5 . 30 Fig. 5 .
Figure 5.Standard deviation (top panels), cubic root of third centered moment (middle panel) and Skewness (bottom panels) for two different SGS source PDFs and the LES.Two-value distribution defined in Equation (12) (continuous lines) and Clipped Gaussian (dot dashed lines).Source surface coverage is 44%.

Figure 6 . 31 Fig. 6 .
Figure 6.Top panels: mean concentration field from the LES (thick dashed line) and the RANS-SFM model with 1 km resolution (thick continuous line).In B two dashed lines are reported, corresponding to 100% surface coverage (thick) and 44% surface coverage (thin).Middle panels: standard deviations for the 44% source surface coverage for the RANS-SFM model with 1 km resolution and the LES.Bottom panels: skewness for the 44% source surface coverage for the RANS-SFM model with 1 km resolution and the LES.