On the Convergence and Capability of the Large-Eddy Simulation of Concentration Fluctuations in Passive Plumes for a Neutral Boundary Layer at Infinite Reynolds Number

Large-eddy simulation (LES) experiments have been performed using the Parallelized LES Model (PALM). A methodology for validating and understanding LES results for plume dispersion and concentration fluctuations in an atmospheric-like flow is presented. A wide range of grid resolutions is shown to be necessary for investigating the convergence of statistical characteristics of velocity and scalar fields. For the scalar, the statistical moments up to the fourth order and the shape of the concentration probability density function (p.d.f.) are examined. The mean concentration is influenced by grid resolution, with the highest resolution simulation showing a lower mean concentration, linked to larger turbulent structures. However, a clear tendency to convergence of the concentration variance is observed at the two higher resolutions. This behaviour is explained by showing that the mechanisms driving the mean and the variance are differently influenced by the grid resolution. The analysis of skewness and kurtosis allows also the obtaining of general results on plume concentration fluctuations. Irrespective of grid resolution, a family of Gamma p.d.f.s well represents the shape of the concentration p.d.f. but only beyond the peak of the concentration fluctuation intensity. In the early plume dispersion phases, the moments of the p.d.f. are in good agreement with those generated by a fluctuating plume model. To the best of our knowledge, our study demonstrates for the first time that, if resolution and averaging time are adequate, atmospheric LES provides a trustworthy representation of the high order moments of the concentration field, up to the fourth order, for a dispersing plume.


Introduction
The dispersion of substances from a small punctiform source in an atmospheric turbulent flow is a physical phenomenon of capital importance in many ecological, environmental, and industrial applications. The concentration field of the dispersing substance displays a random behaviour that, for a specific point in space and time, can be fully characterized by its single-point probability density function (p.d.f.) (e.g., Pope 2000). The mean state of many processes with only a linear dependence on the concentration can be well characterized by the knowledge of the first moment of this p.d.f., the mean concentration. However in other cases, e.g., olfactory processes (Balkovsky and Shraiman 2002;Schauberger et al. 2011), or the release of toxic and flammable substances (Hilderman et al. 1999), non-linearity is observed. In such cases, knowledge of the higher moments of the concentration p.d.f. is needed. Certain modelling approaches such as the two-particle Lagrangian model (Durbin 1980;Thomson 1990;Franzese and Borgas 2002), the Eulerian model resolving the concentration-variance balance equation (Milliez and Carissimo 2008;Yee 2009), and heuristic Lagrangian methods based on single particle models (Ferrero et al. 2016;Cassiani 2013;Manor 2014) are able to estimate only the first two moments of the concentration. The meandering plume approach (Gifford 1959;Cassiani and Giostra 2002) can potentially estimate all concentration moments and the p.d.f.. However, the assumption behind the model formulation becomes less justified further away from the source where empirical parametrization and assumptions about the p.d.f. shape have to be made (Yee and Wilson 2000;Luhar et al. 2000;Cassiani and Giostra 2002;Franzese 2003;Mortarini et al. 2009;Marro et al. 2015). To the best of our knowledge, there are two modelling methods that potentially allow the explicit simulation of the concentration p.d.f. and its higher moments in real world conditions and at very high Reynolds numbers (Re = U δ/ν with U a mean flow velocity scale, δ the boundary-layer height and ν the viscosity) typical of atmospheric flows. These models are the Lagrangian or Eulerian p.d.f. (micromixing) method (Pope 2000;Luhar and Sawford 2005;Cassiani et al. 2005a, b, c;Garmory et al. 2006;Cassiani et al. 2007Cassiani et al. , 2010Postma et al. 2011;Amicarelli et al. 2012;Leuzzi et al. 2012) and the large-eddy simulation (LES) method Sykes and Henn 1992;Weil et al. 2004;Xie et al. 2004a;Vinkovic et al. 2006;Xie et al. 2007;Philips et al. 2013).
Among these modelling approaches, LES is nowadays often viewed as the reference method to simulate the atmospheric flow and dispersion. Large-eddy simulation provides access to the three-dimensional turbulent flow field and it is sometimes used as a replacement for experimental measurements at high Reynolds number. Large scales of the turbulent flow are explicitly simulated by LES but the subfilter scales need to be parametrized using a subgrid scale (SGS) model (e.g., Deardorff 1973;Moeng 1984;Meneveau et al. 1996;Pope 2000;Celik et al. 2009). Implicit or explicit approaches to LES filtering are possible (e.g., Sagaut 2000;Celik et al. 2009) based on the decomposition of the unknown filtered correlations involving or not explicit filtering operations. In the context of implicit filtering the filter width enters into the modelling of the SGS stress tensor and depends on the grid size but it can still be made larger than the grid size to ensure results are independent of the grid and numerics (e.g., Mason and Callen 1986;Geurts and Frohlich 2002;Pope 2004;Geurts 2006). However, in most applications the filter is implicitly defined to be exactly or very close to the grid size, which makes the LES results dependent on the grid definition and on the numerical methods used (e.g., Pope 2004;Geurts 2006;Kemenov et al. 2012).
In engineering applications, various methods have been proposed to estimate the quality of the turbulent flow simulated by an LES model (e.g., Celik et al. 2005Celik et al. , 2009). However, we note that these methods are mostly suited for relatively low Reynolds number flows as they involve the actual viscosity of the fluid and convergence to direct numerical simulation (DNS) results. Large-eddy simulation of atmospheric flow is different as it involves the further assumption of infinite Reynolds number, totally neglecting the molecular viscous term in the equations, and the use of a wall-stress model (e.g., Moeng 1984;Porté-Agel et al. 2000;Chow et al. 2005;Sullivan and Patton 2011;Maronga et al. 2015). This implies that the viscous layer is totally unresolved or that the surface is rough but with unresolved roughness scale (Brasseur and Wei 2010).
The dispersion of scalars introduces further approximations associated with the numerical methods and SGS model used for the scalars. Subgrid-scale scalar fluctuations in LES can be obtained from algebraic and transport equation methods (e.g., Colucci et al. 1998;Pierce and Moin 1998;Jimenez et al. 2001;Balarac et al. 2008;Kaul et al. 2009). In atmospheric dispersion applications, the SGS scalar fluctuations are usually not explicitly modelled by a transport equation (e.g., Mironov et al. 2000;Heinze et al. 2015).
Despite the uncertainties related to the filter and numerical implementation, in recent years LES has become a tool used in many applied atmospheric dispersion studies, including the urban environment, and for critical applications such as the release of toxic gas substances (e.g., Fossum et al. 2012;Nakayama et al. 2013;Lateb et al. 2016). However, very few thorough evaluations of LES for plume dispersion and concentration fluctuations from small (point-like) sources in very high Reynolds number boundary layers are available. Henn and Sykes (1992) and Sykes and Henn (1992) studied dispersion and concentration fluctuations in the convective and neutral boundary layers, respectively. The maximum achievable resolution was limited by the available computing power, but source sizes smaller than grid resolution could be simulated by using a puff model for scalar dispersion with a semiempirical parametrized SGS puff expansion (Sykes et al. 1984). This was empirically set by Sykes and Henn (1992) to match the Fackrell and Robins (1982) wind-tunnel experiment, thus allowing a satisfactory comparison with this dataset. Xie et al. (2004aXie et al. ( , 2007 investigated plume dispersion and concentration fluctuations (including extremes) in a neutral boundary layer with the highest achievable LES grid resolution capable to resolve the scalar source size with one grid cell. More recently, Boppana et al. (2012) studied the relatively simpler case of dispersion from a line source in channel flow at a finite Reynolds number of 10800 and still underlines the difficulty implied in investigating scalar dispersion from small sources by means of LES, and the need of thorough validations.
Validating LES involving scalar fluctuations for a dispersing plume, at the very high Reynolds number typical of atmospheric flows, is a complicated task. For instance, the comparison with DNS is not feasible. Atmospheric dispersion experiments (e.g., Mylne and Mason 1991;Mylne 1992Mylne , 1993Jørgensen and Mikkelsen 1993;Yee et al. 1993a, b;Mole and Jones 1994) providing concentration-fluctuation statistics have some characteristics that make them unsuitable for the purpose of validating LES results. These are, (1) the extremely small scalar sources used which are unreachable even for nowadays LES, (2) sensors located quite far downwind from the source so that the effects related to the source size have been mostly forgotten, and (3) emissions and measurements are located very close to the ground, where the effect of the wall-model parametrization in the LES are critical because most of the energy may not be explicitly resolved. A viable alternative, also used in previous LES studies that are mentioned above, is to use the few available wind-tunnel dispersion studies resembling atmospheric boundary-layer conditions (Fackrell and Robins 1982;Nironi et al. 2015). In these experiments the emitting sources are small compared to the boundary-layer thickness but within LES possibility (comparable to a few metres in the atmosphere) and are also placed at relatively high elevation where the LES explicitly resolves most of the turbulent kinetic energy (TKE). One general advantage of the wind-tunnel experiments compared to atmospheric measurements is that they are non-affected by unsteadiness and can use very long averaging time, which for higher moments of concentration is necessary. One potential disadvantage is that wind-tunnel data may be affected by an Re dependence if Re is not very high.
As recognized by many authors, evaluation of LES results should always involve a wide range of grid resolutions (e.g., Pope 2004;Celik et al. 2005;Klein 2005;Klein et al. 2008;Sullivan and Patton 2011;Kemenov et al. 2012) but these studies are rare in the literature. To our knowledge, there exists no study that systematically examines the grid dependence of the LES results for concentration fluctuations from a continuous small finite (point-like) source in an infinite-Re neutral boundary-layer configuration.
The tendency of the error to decrease towards zero by increasing grid resolution can potentially be investigated for low Reynolds number LES with explicit wall simulation that converge to a DNS solution when their resolution increases (e.g., Kemenov et al. 2012). However, as mentioned above, typical atmospheric LES are not part of this category because the Reynolds number is far above the reach of DNS feasibility. Nonetheless, exploring the dependence of the results on the grid resolution and comparison to experimental measurements are fundamental for determining whether the LES results for the scalar converge, in some sense, and to evaluating the expected range of variability in the results.
We take a heuristic but rational and comprehensive approach and explore many statistics, and characteristic scales of the velocity and scalar fields to find evidences of convergence in the turbulent fields of our LES. The analysis of the statistical characteristics of the fluctuating concentration field includes moments up to the fourth order and the concentration p.d.f. with an unprecedented level of detail for a LES. A clear indication on the range of validity of the Gamma p.d.f. model (e.g., Duplat and Villermaux 2008) for the concentration p.d.f. is also obtained. We use the freely available and widely used open source Parallelized LES Model (PALM) (Maronga et al. 2015) to make our analysis directly useful towards practical applications.
The paper is organized as follows. Section 2 provides a description of the numerical method and the simulated cases, and subsequently, turbulent flow statistics are discussed in detail in Sect. 3.1. The spectra of the TKE and two-point statistics are addressed in Sect. 3.2, and the analysis continues with the investigation of the scalar field. The mean concentration, concentration variance and its budget are discussed, respectively in Sects. 4.1 and 4.2, while the mechanisms generating concentration fluctuations, and how this are altered by not adequate resolution, are examined in Sect. 4.3. In Sect. 4.4, we examine the shape of the concentration p.d.f. using the scaled moments of scalar concentration, including the ratio of standard deviation and mean, the skewness and the kurtosis. Finally, a summary and discussion are presented in Sect. 5.

Methods
PALM is an open source model code developed mainly at University of Hannover (Maronga et al. 2015). Here the code is used to solve the non-hydrostatic, filtered incompressible Navier-Stokes equations in Boussinesq-approximated form at formally infinite Reynolds number (Re = u ∞ δ/ν, with u ∞ the freestream velocity) due to the neglect of molecular viscous stress (e.g., Geurts and Frohlich 2002;Stevens et al. 2014), with an advection-diffusion equation for the transport of a passive scalar. Details of the governing equations are reported in Appendix 1 for completeness.
For the advective term in the incompressible LES model equations, the Piacsek and Williams (1970) second-order, formally energy conserving scheme, was chosen over the fifth-order dissipative scheme proposed by Wicker and Skamarock (2002). Heinze et al. (2015) compared these two schemes in PALM and found that the Wicker and Skamarock (2002) scheme is much more dissipative, artificially reducing the scale-interaction term and therefore the transfer of kinetic energy from the resolved scale part of the energy spectrum to the subgrid scale part. Indeed the use of high-order finite-difference dissipative schemes in LES generates excessive damping of the high frequencies and often shows numerical dissipation larger than subgrid scale dissipation (e.g., Beudan and Moin 1994;Mittal and Moin 1997;Sagaut 2000;Park et al. 2004). This can be avoided if the mixing length is made suitably larger than the grid size (e.g., Mason and Callen 1986) but this is seldom done, especially in practical applications, as grid resolution is usually at the limit of available computational resources. On the other hand the Piacsek and Williams (1970) second-order scheme has been successfully used in LES (e.g., Sykes and Henn 1988;Cuijpers and Duynkerke 1993;Mason and Brown 1999;Brown et al. 2000;Dosio et al. 2003). For the dispersing scalar plume, PALM allows the use of two possible schemes. The monotone, locally modified version of Bott's advection scheme proposed by Chlond (1994) is used, ensuring positive definite scalar values and mass conservation. The non-monotone scheme of Wicker and Skamarock (2002) is also available, but was found to be inadequate for simulating plume dispersion A three-dimensional domain of 4.8 m× 0.8 m× 0.8 m, respectively in the along-wind (x), crosswind (y) and vertical (z) directions is simulated to reproduce the wind-tunnel experiment of Nironi et al. (2015) where a boundary-layer thickness of δ = 0.8 m was measured. Some aspects of the turbulent field are also comparable to the wind-tunnel experiment of Fackrell and Robins (1982) (hereafter F&R) if appropriately scaled, as extensively discussed in Nironi et al. (2015). Therefore, the measurements of F&R will also be used when possible.
The neutral boundary layer is simulated as an incompressible half-channel flow at infinite Reynolds number with a strictly symmetric, stress-free condition at the top of the domain (∂u/∂z = 0, ∂v/∂z = 0). The flow is driven by a constant mean pressure gradient. Although this flow configuration is not realizable in nature, it is often used in order to simulate windtunnel-generated boundary layers and the atmospheric boundary layer, neglecting Coriolis effects (e.g., Schumman 1989;Porté-Agel et al. 2000;Xie et al. 2004b;Cassiani et al. 2008;Brasseur and Wei 2010;Margairaz et al. 2018). The constant mean pressure gradient is defined as ∂ p/∂ x = −u 2 * /δ, where u * = 0.185 m s −1 is the friction velocity measured in the wind tunnel. On the bottom wall the roughness length z 0 = 1.1 × 10 −4 m is also set equal to that estimated in the wind tunnel of Nironi et al. (2015). The roughness length enters in the wall model as the wall is not explicitly resolved but a constant-flux layer is used as is commonly done in atmospheric simulations (e.g., Moeng 1984;Maronga et al. 2015). For the velocity, periodic boundary conditions are used on the lateral boundaries, while non-periodic boundary conditions are set for the passive scalar (see Appendix 2 for more information about the boundary conditions).
The number of nodes of the computational grid (N x × N y × N z , with N x , N y , and N z being the number of grid points in along-wind, crosswind, and vertical directions, respectively) is ranging from 256 × 64 × 64 ≈ 10 6 nodes to 2048 × 512 × 512 ≈ 5 × 10 8 nodes. The size of the source is d s = 12 mm = 0.015δ in the vertical and crosswind directions. This initial source size is in the middle of the range of the source sizes investigated in F&R, where d s /δ = 0.0025, 0.007, 0.0125, 0.0208, 0.0291 were considered, but it is larger than the range investigated in Nironi et al. (2015), d s /δ = 0.00375 and 0.0075. This initial source size was chosen so that the grid convergence of the statistics is fully investigated among four different grid resolutions used to resolve the flow. The emission source is simulated with exactly one grid cell in N x = 256, 2 3 , 4 3 , and 8 3 grid cells in N x = 512, N x = 1024, and N x = 2048 simulations, respectively. Hereafter, the four simulations are referred to according to their N x value. The source is located at y s , corresponding to the centre of computational domain in crosswind direction, and at the elevation of z s = 0.19δ which corresponds to the elevation used in both Nironi et al. (2015) and F&R. Table 1 lists some important quantities of the wind-tunnel and numerical experiments.
For a horizontally homogeneous boundary-layer flow discussed here, both time averages and plane averages can be used to calculate one-point flow statistics. However, the scalar field from a small (point-like) source is fully non-homogeneous with only a statistical symmetry, with respect to source location along the y direction, potentially usable to increase the samples. In the following, we will use only time averages to calculate the scalar statistics unless otherwise explicitly stated and explained. Large averaging times are necessary to obtain convergence in higher-order statistical moments if the plane average is not used. The averaging time used for the statistics is 150 s, i.e. between 600 and 800 times the estimated Lagrangian time scale as further discussed below. A spin-up time of 120 s was used to ensure that flow statistics are in steady state before starting the time average. These spin-up and averaging times imply that about 8 × 10 5 core hours were used typically for one 2048 simulation.
In the following we adopt a standard notation with the overbar() denoting a resolved scale (filtered) variable, the single prime () a subfilter scale fluctuation, the angle brackets () a space and/or time average and the double prime () a fluctuation from this average. Any flow variable φ can be decomposed as: φ = φ +φ + φ . Meteorological or index notation are used as convenient so u 1 = u, u 2 = v, u 3 = w represent the velocity components in the along-wind x 1 = x, crosswind x 2 = y, and vertical x 3 = z, directions respectively. Vectors are represented in bold characters, e.g. x = (x 1 , x 2 , x 3 ).

Turbulent Velocity Field
To correctly interpret the scalar dispersion it is necessary to understand the statistical characteristics of the velocity field. Here the first-and second-order statistical moments of velocity and turbulent structures are examined, including their change with grid resolution.

Turbulent Flow Statistics
The LES flow is compared with the wind-tunnel measurements of Nironi et al. (2015) and F&R. Averaging is done over the horizontal plane and in time. The mean wind speed at the top of the boundary layer (u ∞ ) is equal to 4.8 m s −1 in the LES; it is 5 m s −1 in Nironi et al. (2015), and 4 m s −1 in F&R, as listed in Table 1. We remark that the Nironi et al. (2015) experiment has the same values for roughness length, boundary-layer thickness and a very similar friction velocity as in our simulations and should also have very similar mean wind profiles having very similar u ∞ . Figure 1a shows the mean wind profile ū scaled by friction velocity. The LES results are weakly sensitive to resolution with small but visible changes. By increasing the grid resolution, the mean wind speed increases at the source elevation, although it tends to be similar at the top of the boundary layer (see Table 1). The LES overestimates mean wind speed profiles compared to the wind-tunnel measurements in the lower part of the domain. This Table 1 Boundary-layer characteristics: freestream velocity u ∞ , friction velocity u * , boundary-layer thickness δ, roughness length z 0 , mean velocity at source elevation u s , Reynolds number of the flow Re = u ∞ δ/ν, Reynolds number based on SGS viscosity (ν SG S ) in the middle of the boundary layer, Re overestimation of mean wind speed in LES using wall models, compared to the logarithmic law profile, has been investigated by many authors (Xie et al. 2004b;Brasseur and Wei 2010;Hultmark et al. 2013;Ercolani et al. 2017) and some correction methods have been proposed (e.g., Xie et al. 2004b;Hultmark et al. 2013). The factors influencing this aspect are the SGS model, the grid aspect ratio and the wall model. Based on the SGS model closure, grid resolution and grid aspect ratio, all the simulations here but the coarsest fall in the high accuracy zone estimated by Brasseur and Wei (2010), their Fig. 7, in relation to the LES capability of predicting the logarithmic law of the wall. More recent literature (Ercolani et al. 2017) found that the SGS model used here gives overestimation of mean wind speed when used with an almost isotropic grid aspect ratio between one and two, and this is confirmed by our results. Ercolani et al. (2017) suggested that the overestimation is linked to the over dissipative nature of the SGS model and that this behaviour is somewhat masked by the use of an anisotropic grid with an aspect ratio of (Δx = Δy)/Δz = 4. This was not attempted here as the recommended aspect ratio is not ideal for the plume dispersion modelling that ideally requires the crosswind and vertical grid size to be the same. The wall model used in PALM (see Appendix 2) assumes that the log law is valid for the instantaneous resolved wind velocity magnitude. Wall models belonging to this class are commonly used in atmospheric LES (e.g., Moeng 1984;Maronga et al. 2015) and are known to overestimate the mean wind speed compared to the log law as discussed in Hultmark et al. (2013). However, again from the plume dispersion perspective, which is our main focus, the difference in the mean wind speed at source elevations is small and this is the main advection velocity for plumes. The difference of mean wind speed among different resolutions is quite small here, but we mention that an aspect that influences the convergence (among grid resolutions) of mean wind profiles in the wall-model LES is the height used to inject information in the wall model (Kawai and Larsson 2012). The PALM code uses the standard LES approach and takes advantage of the available grid resolution feeding information to the wall model at Δz/2. Kawai and Larsson (2012) found that convergence is improved if the information is transferred to the wall model at a constant height irrespective of grid resolution. Figure 1b shows the SGS and total TKE. The SGS TKE (e = 1 2 u i u i ) is calculated in PALM by solving a model equation (e.g., Moeng 1984;Maronga et al. 2015) as it is common in many atmospheric LES models. The total TKE is E +e where E is the resolved scale TKE, E = 1 2 ū iū i . The total TKE is in very good agreement with both experimental datasets at all grid resolutions. Pope (2000) suggested that LES is properly resolved if at least 80% of TKE is resolved, and for LES with wall modelling this is possible only away from the surface. For the current simulations this condition is fulfilled for z/δ 0.05, 0.02, 0.01, 0.005, respectively, moving from coarser to finer resolved simulations. Therefore, at the elevation where the plume is released most of the energy is explicitly resolved in all cases. This also means that numerical dissipation plays an important role. Figure 1c shows the mean dissipation rate ( E ) of TKE as obtained from the residual of the resolved TKE balance, Eq. 1,(e.g., Mironov et al. 2000), where horizontal homogeneity and steady state have been used. In this equation, p represents the perturbation pressure (see e.g. Moeng 1984;Maronga et al. 2015, and Appendix 1 for more details). The dissipation rate E , is the last term on the right-hand side (r.h.s), and has been evaluated as a residual. The agreement between this estimation of the dissipation and the experimental value of F&R is extremely good. The experimental values of Nironi et al. (2015) shows a just slightly higher dissipation. We remark again that the values in Fig. 1c pertain to elevations where most of the TKE is resolved (z/δ ≥ 0.05) and therefore Eq. 1 is well representative of the total TKE. We found that other estimates of the mean dissipation are not representative of the actual dissipation. This agrees with Heinze et al. (2015), who used the PALM code to compare three methods of estimation of the TKE dissipation; (i) the parametrized form in the SGS TKE equation where Δ = 3 √ ΔxΔyΔz with Δx, Δy, Δz being the grid spacing in the x, y and z directions, respectively and l is the SGS mixing length; (ii) the computation of the scale-interaction term describing the transfer of energy between resolved and subgrid scales (Brown 1994;Heinze et al. 2015;Mironov and Sullivan 2016); and (iii) the residual of the TKE budget. They found that the only reliable method to compute energy dissipation is from the residual of the TKE budget. This is not surprising, as the latter method includes numerical dissipation, which can be a dominating fraction of dissipation where LES resolves most of the TKE. This also agrees with Glendening and Haack (2001) who compared several numerical schemes and found that in the presence of strong mean advection only pseudo-spectral methods were able to correctly avoid spurious dissipation at high wavenumbers. We note that Heinze et al. (2015) used the total TKE balance, but for the elevation considered here most of the energy is resolved. For values z/δ 0.05, 0.02, 0.01, 0.005, respectively, moving from coarser to finer resolved simulations the SGS dissipation starts to dominate over the resolved scales (not shown here).
Figure 1d-f show the scaled resolved variance of the three velocity components. The variance of the resolved scale along-wind component, σ 2 u , progressively gets closer to the experimental values when increasing the grid resolution, and for the highest resolutions the agreement is very good. The crosswind and vertical velocity variance similarly become closer to the wind-tunnel experimental values when increasing the grid resolution, but the agreement is not as good as for the along-wind velocity component. Similar values of crosswind velocity variance were obtained in the LES of Xie et al. (2004a) using the SGS mixed-scale model of Sagaut (1995). Porté-Agel et al. (2000) obtained similar values of crosswind velocity component using a Smagorinsky model closure, but higher values using a scale dependent dynamic model. The two wind-tunnel experiments have very similar values, except for the crosswind component (v), where F&R is more similar to the LES results compared to Nironi et al. (2015). The Reynolds stress, ū w (not shown here) is a straight line as expected and matches very well the experimental data. In the highest resolution case, the observations are nearly captured, even without the contribution of the SGS stress. Obviously in the lowermost layer close to the wall everything is subgrid.
It can be anticipated that by considering the difference in the crosswind and vertical velocity components variance (Fig. 1e, f), a lower plume spread may be expected for the LES when compared to the wind-tunnel data of Nironi et al. (2015). However, overall the second moments of the high resolution LES are similar to the wind-tunnel experimental values. Plume dispersion statistics in neutral boundary layers are mainly driven by second-order velocity statistics and by velocity length and time scales, which are the footprint of turbulent velocity structures. In the next section, we complete the velocity field analysis using the spectrum, two-point correlation and turbulent length scales across LES resolutions before turning our attention to the scalar field.

Turbulent Flow Structures and Length Scales
To examine the existence and extension of the inertial subrange, the two dimensional spectra of the TKE on the horizontal plane are now considered. As discussed in Sullivan and Patton (2011) the two-dimensional spectra are more representative of the spatial eddy scale than one-dimensional spectra. Figure 2 shows the time-averaged spectrum of the TKE E(k h , z) (e.g., Pope 2000), where k h is here the horizontal dimensionless wave vector module k h = kδ(2π) −1 with k being the horizontal angular wave vector module (rad m −1 ) (e.g., Wyngaard 2010). The spectrum is obtained by integrating circles of increasing radius and it is presented for all grid resolutions both at the source elevation, z = 0.19δ, and at the middle of the boundary layer, z = 0.5δ. The k −5/3 h Kolmogorov inertial subrange scaling is also shown as a straight line. The inertial subrange is visible only in the two higher resolved simulations (1024 and 2048), which therefore can be assumed to converge in this sense. The lowest resolution seems especially under-resolved, lacking the display of an inertial subrange.
Large-scale turbulent structures are fundamental for plume dispersion as they define the persistence of spatial and time correlations. In an inhomogeneous and anisotropic shear flow a variety of integral length scales exists to describe these structures. These length scales can be estimated from the two-point spatial correlation coefficient, which is defined as (e.g., Carlotti and Drobinski 2004;Nironi et al. 2015) where x is a fixed point and r is a generic vector. In vertically inhomogeneous turbulence, the two-point correlation depends on the vertical position. However, the most significant position to characterize plume dispersion for an elevated small (point-like) source is certainly that of the source. Therefore, the analysis here is about the correlations centred at source position and limited to the single velocity component correlation coefficients, ρ ii (no summation implied in repeated index here), which are used to obtain the relevant integral length scales. As extensively discussed in Carlotti and Drobinski (2004) many integral scales can be defined in inhomogeneous anisotropic turbulence by integrating the spatial correlation coefficient over the distance r along all the possible orthogonal directions, where e j is the unit vector in x, y, z directions; x s is the source position; and no summation is implied in repeated index. In the present case of wall bounded flow, the only non homogeneity and possible asymmetry is in the vertical direction (Carlotti and Drobinski 2004). Following Nironi et al. (2015) our analysis is limited to L uu,x = L 11,1 , L vv,y = L 22,2 and L ww,z = L 33,3 , with the latter two being the more relevant for understanding plume dispersion in the crosswind and vertical directions . The computed length scale values are reported in Table 2. The values of the length scales are obtained by fitting the exponential function of the form of exp(∓r /L ii, j ) to the profiles of the correlation coefficient, where the signs depend on the direction being positive or negative, respectively. These profiles must be symmetric except for the vertical correlation coefficient (ρ 33 ) where the wall effect causes asymmetry. Considering the average value of the length scales in Table 2 Estimate of the Eulerian length scales from the two-point spatial correlation coefficient, by fitting the exponential function of the form exp(∓r /L ii, j )e j to the profiles of the correlation coefficient. The first value in L ww,z column corresponds to the positive distance +r e j and the second value is for negative distance −r e j as specified in Eq. 4 Resolution

Turbulent Scalar Field
The transport of the passive scalar in PALM is simply obtained by solving an advection diffusion equation for the resolved scalar field, where the SGS scalar diffusivity is defined as K s = 1 + 2l Δ K m , and l = min(1.8z, Δ) ( Moeng and Wyngaard 1988;Maronga et al. 2015). As remarked in Sect. 2, the Bott advection scheme as modified by Chlond (1994) has been used in this study. Figure 3 shows snapshots of the concentration field for all the resolutions on the (x, y) plane at z = z s = 0.19δ. There are obvious visual differences in the scalar field up to 1024 resolution, while between 1024 and 2048 these differences become more subtle. Interestingly, this visual impression is somewhat confirmed by the quantitative analysis below. Time series of the concentration at the location of maximum mean concentration are shown in Fig. 4 for x/δ = 0.13 and x/δ = 3.75. Similarly to what was observed in the snapshots, the differences are rather obvious up to 1024 resolution and become more subtle between 1024 and 2048.

Mean Concentration Field
Crosswind and vertical profiles of the mean concentration c , through its maximum for different grid resolutions at various downwind distances are plotted in Fig. 5. Angle brackets represent only the time average when the scalar is involved. The concentrations are reported following Nironi et al. (2015) by using the scaling c * = c(u s δ 2 /Q) where Q = 1 kg s −1 is the mass flow rate at the source located at z s . We would expect that the mean concentration is minimally sensitive to the grid resolution in LES, as the mean concentration is influenced by the whole turbulent spectrum and the LES should explicitly capture most of the energy. Contrary to this intuitive picture some difference among grid resolutions is visible. The 2048 simulation stands out with a much higher spread and lower peak concentrations compared with lower grid resolutions for all downwind distances but the first one, and this is more evident for x/δ = 3.75. On the contrary, at the first downwind distance the numerical diffusion dominates and the 256 and 512 resolutions show low concentration values. Figure 6a shows the scaled maximum mean concentration as a function of the equivalent along-wind distance from the source x * defined below. Here we include in the comparison the experimental measurements of Nironi et al. (2015) and F&R. As reported in Table 1, the simulations and the experiments have different mean velocity at the source; this in turn implies that the plume advection time at a certain distance from the source is different. Following Nironi et al. (2015), an appropriate dimensionless advection time can be defined accounting for different mean velocity at the source and different scaling parameter u * and δ as The comparison of plume results from different experiments should be done at the same dimensionless time, i.e., in our case T * (ex p) = T * (L E S) . This equality allows the definition of an equivalent dimensionless advection distance where the experimental results have the same dimensionless advection time as the LES results, Based on this equality, we define the equivalent distance x * as Similarly to the use of the dimensionless advection time, using this distance eliminates the difference originating simply by the mean advection velocity difference and allows also the comparison of plume dispersion from experiments with different u * and δ. The LES results and the data of Nironi et al. (2015) have the same u * and δ, so the scaling only differs through the mean velocities, u s , but is limited since the difference in velocity is small. For u s(L E S) , needed in Eq. 8, the average value between the two most resolved simulations was used. We note that the LES results have not been corrected for the slightly different advection velocities as we want to underline the difference originating from grid resolution in the same exact simulation settings (the correction would be anyway minimal). In Fig. 6a the measurements of Nironi et al. (2015) and F&R have a source size that is smaller than (about half of) that used in the simulations, but source size has a very limited effect on mean concentration and can be neglected here. The model results are in quite good agreement with the measurements, especially those of Nironi et al. (2015), but the maximum mean concentration is generally overestimated and therefore plume spread must be underestimated. This is consistent with the lower simulated variances in the crosswind and vertical velocity components as shown in Fig. 1e, f. Figure 6b, c show the calculated crosswind and vertical plume spread standard deviation as a function of the distance from the source. These standard deviations are defined as the square root of the following variances with and σ 2 z are also commonly called the absolute plume dispersion (e.g., Arya 1999; Dosio and de Arellano 2006). In the main panels of Fig. 6b, c, it is evident that the 2048 simulation has a larger averaged plume spread compared to the other simulations for σ i /δ larger than about 0.06, and that the simulations have generally lower averaged plume spread compared to the wind-tunnel measurements apart from the closest location. The insets in 6b, c show the near source behaviour. In this case the two less resolved simulations have a higher spread and this is due to numerical diffusion as further discussed in Sect. 4.3. Table 3 Estimate of the Lagrangian time scales (T lv , T lw ) from a fit of Taylor's dispersion relations (Eq. 10) to the plume spread in crosswind and vertical directions for all grid resolutions. The estimates of the Lagrangian time scales as the ratio of length scales and standard deviations of velocity components following Tennekes and Lumley (1972)  To explain the mean concentration behaviour it is useful to recall an approximation of plume-spread standard deviation commonly used in conjunction with Gaussian dispersion models. Using the slender plume approximation (e.g., Arya 1999) and adapting Taylor's dispersion theory (Taylor 1922) to anisotropic turbulence the evolution of the plume-spread standard deviation, in crosswind (σ y ) and vertical (σ z ) directions can be approximately written as where σ s is the source standard deviation, σ 2 v and σ 2 w are the variances of the crosswind and vertical velocity components, and T Lv and T Lw are the crosswind and vertical Lagrangian integral time scales, respectively (e.g., Arya 1999). For an elevated plume this approximation is commonly used to describe plume dispersion by taking all the values as those at the source elevation and replacing time and downwind distance according to Taylor's hypothesis t = x/ u . These formulations remind us (e.g., Arya 1999) that for t << T L turbulent dispersion is insensitive to T L , while for t >> T L standard deviation of plume dispersion increases proportionally to (T L t) 1/2 . Close to the source the mean concentration is mainly influenced by the velocity variance, the mean plume increases as σ 2 y,z = σ 2 s +σ 2 v,w t 2 (e.g., Arya 1999). As discussed above, the 2048 simulation has slightly greater variances (σ 2 v , σ 2 w ) but considerably larger structures (see Table 2) compared to lower resolutions. Larger structures result in longer correlation time scales including Lagrangian autocorrelations. For a neutral boundary-layer flow following Tennekes and Lumley (1972) the two-point length scales can be directly related to the Lagrangian integral time scales as T Lv ≈ L vv,y /σ v , T Lw ≈ L ww,z /σ w , T Lu ≈ L uu,x /σ u . Table 3 shows the Lagrangian time scale estimate according to Tennekes and Lumley (1972) and the estimates obtained by fitting Taylor's dispersion relation (Eq. 10) to the plume spread in the crosswind and vertical directions. These are all rather crude estimates, especially in the vertical direction, but are adequate for cross comparing grid resolutions. The 2048 simulation has larger structures and accordingly longer Lagrangian time scales, especially in the crosswind direction, compared to the lower-resolution simulations, which leads to the lower maximum mean concentration of the scalar.  Figure 7 shows the crosswind and vertical profiles, passing through the position of maximum mean concentration, for the scaled standard deviation of concentration σ * c = σ c (u s δ 2 /Q), σ c = ( c 2 ) (1/2) . The effect of grid resolution on scalar fluctuations is evident in the standard deviation (Fig. 7). The change in the mesh resolution between the 256 and 1024 simulations does not produce significant change in the mean field while the scalar fluctuations are more than doubled. The increase is also substantial between the 512 and 1024 grid specifications. In contrast, 2048 and 1024 grid configurations seem to converge to quite similar results and are the only simulations able to capture the presence of the double peak at x * /δ = 0.13. The behaviour in Fig. 7 can be summarized stating that the standard deviation converges at the 1024 grid configuration. This behaviour is in striking contrast with what was observed for the mean concentration where all the simulations but the 2048 one gave quite similar results. This is due to the fact that the mechanisms driving the evolution of the mean and the variance of the concentration are different and are differently influenced by the numerical resolution.

Concentration Fluctuations Variance and Budget
The budget analysis of the resolved scalar variance is introduced to highlight the relevant terms defining the observed concentration variance. The budget is calculated at the position of maximum mean concentration where most of the TKE is resolved for all resolutions. In a stationary state, the budget equation for the resolved scale mean scalar variance is written as where the first term on the r.h.s. corresponds to the advection (Adv.), the second term corresponds to the production (Prod.), the third one corresponds to the turbulent transport (TT) and ξ res is the mean scalar dissipation from the residual. The definition of the scalar dissipation residual follows the same logic applied above to the resolved scale TKE budget and it is comprehensive of numerical dissipation. Two additional estimates of the mean scalar dissipation are possible: one is based on the equilibrium approximation for the SGS Kaul et al. 2009;Heinze et al. 2015) and reads, Alternatively the transfer of resolved-scale scalar variance to the SGS may be defined as Mironov and Sullivan 2016), where τ ci = u i c is the SGS scalar flux computed by the SGS model. These two estimates have been calculated and considered in the analysis. However, they do not account for numerical dissipation. Therefore, similarly to what was discussed above for the TKE, they cannot be considered reliable. Figure 8a-c show the budget at the elevation of maximum mean concentration in the crosswind direction for the 2048 simulation. All the terms are indicated as φ, and are reported normalized on the ordinate as φ * = φ δ u * u s δ 2 Q 2 . The production term, turbulent transport and advection are all important close to the source (x/δ = 0.13) but already at x/δ = 0.725 the budget is dominated by advection and turbulent transport. At x/δ = 1.92 advection dominates and dissipation (ξ res ) becomes as large as turbulent transport. This analysis agrees and completes what was found by F&R, that for distances x/δ > 1.92 advection and dissipation dominates the balance. Figure 8d-f show the budget at x/δ = 0.13 but for the grid resolutions 256, 512 and 1024. Production is the dominant term but for 256 and 512 simulations it has much lower values compared to the 2048 one (Fig. 8a). The shape of the terms is quite different in the 256 and 512 simulations compared to the 2048 one. In contrast, the 1024 simulation shows quite similar values and shape to those of the 2048 simulation. However, the production is still lower. This confirms that a tendency to converge is observed between 1024 and 2048 resolutions for the concentration variance, but the convergence is not yet perfect.
This analysis demonstrates that production close to the source is the critical phase for plume concentration fluctuations variance and that only an appropriate resolution can capture it. As the production of concentration variance occurs close to the source while afterward advection is the dominant term, we may conclude that significant converge, albeit not total in the concentration variance, is ensured by the source being resolved by at least 4 3 grid nodes with the current numerical schemes. The next section analyzes the mechanism-generating fluctuations and how the insufficient resolution alters it.

Absolute and Relative Dispersion, Meandering Motions, and Production of Fluctuations
Eddies of a size larger than the plume generate meandering fluctuations (Gifford 1959;Csanady 1973;Franzese and Cassiani 2007;Cassiani et al. 2009), i.e. large-scale undulations of the plume as a whole and the sweeping of the plume centreline with respect to the source location. By contrast, eddies of a size comparable to, but smaller than, the plume generate relative dispersion, that is the spreading of the plume in a coordinate system relative to the instantaneous centre of mass position (Richardson 1926;Batchelor 1952;Monin and Yaglom 1975;Sawford 2001;Franzese and Cassiani 2007). Eddies much smaller than the plume contribute little to relative dispersion (e.g., Mikkelsen et al. 1987). When assuming that the associated spatial scales are separated, the absolute dispersion of the plume can be partitioned between two statistically independent components (Gifford 1959;Csanady 1973;Franzese and Cassiani 2007). Therefore the absolute dispersion defined in (9a), can also be written as The meandering variance can be readily defined from the local centre-of-mass definition, as The relative dispersion can then be calculated from Eq. 14. Absolute dispersion, relative dispersion and meandering, as a function of along-wind distance from the source, are reported in Fig. 9 for the 2048 simulation. Figure 10 represents the comparison of the absolute and relative dispersion between different grid resolutions. At a constant absolute dispersion, the smaller the relative dispersion is, the higher the fluctuations are since the meandering related production of fluctuations is enhanced by two mechanisms, the increased flapping of the Fig. 9 Absolute and relative dispersion, and the meandering motion of the plume, as a function of downwind distance for the 2048 simulation for, a crosswind, and b vertical directions plume and the fact that a smaller and more concentrated plume is moving. This is most important very close to the source, as illustrated in Fig. 10c, d, where most of the production of the scalar variance occurs, as extensively discussed above using the variance budget. Therefore we may hypothesize that a coarse grid resolution increases the relative dispersion by numerical diffusion and artificially damps the meandering production of concentration fluctuations in a way similar to an increase in the source size and this decreases the production of scalar variance. After this initial phase the scalar variance is mainly transported away and eventually dissipated far away from the source.
The hypotheses discussed above can be nicely and quantitatively demonstrated by using the meandering plume model of Gifford (1959) (see Appendix 3 for more details). In this model the concentration moments can be obtained from σ 2 m and σ 2 r , which, in the present case, are those directly obtained from the LES concentration field with no parametrization used. The meandering plume model shows that the production of fluctuations by meandering depends on the ratio M = σ 2 m /σ 2 r and on the value of σ 2 r (see Appendix 3), larger σ 2 r and smaller M generate a decrease in fluctuations.
The second moment of the concentration from LES and from the meandering plume model are compared in Fig. 11 as a function of downwind distance from the source. Dashed lines correspond to the meandering plume model and solid lines correspond to the LES results at various grid resolutions and at the position of the peak of mean concentration. The meandering plume model reproduces well the value of the variance peak generated by the LES concentration time series at different resolutions. Increasing values are observed by increasing resolution both for the meandering model and for the LES. The agreement is better at the highest resolution. The position of the peak in the meandering is slightly displaced compared to the LES but this difference decreases with resolution. Further downwind, as expected, a considerable difference arises between the meandering plume model and the LES because in the fluctuating plume formulation used here any contribution of internal fluctuations, in the coordinate system relative to the local plume centre of mass, is neglected.
The meandering results in Fig. 11 prove the physical picture discussed above showing that the production of fluctuations at coarse resolutions decreases because of artificially enhanced relative dispersion by numerical diffusion. It also shows that the meandering production of fluctuations for the 1024 and 2048 simulations is similar, thus supporting the existence of a convergence of LES results for the concentration fluctuations variance between the 1024 and 2048 resolutions.    Nironi et al. (2015) and Fackrell and Robins (1982) experimental data for several source sizes

Scaled Concentration Moments and Probability Density Function
So far the mean concentration and the variance have been analyzed, and we will now focus on the shape of the concentration p.d.f.. A parameter that is often used to characterize the p.d.f. shape is the relative intensity of concentration fluctuations (see e.g. F&R) defined as i c = σ c / c . This is also called the coefficient of variation in statistical textbooks, and describes the spread of the p.d.f. relative to its mean value. As discussed in, e.g., Yee et al. (1993c), Nironi et al. (2015), Marro et al. (2018) i c suffices to define the shape of the Gamma p.d.f., which is often used as an analytical model of the concentration p.d.f. based on experimental evidences and theoretical arguments on scalar mixing in confined domains (Duplat and Villermaux 2008). So far it was observed that the N x = 2048 simulation has different flow structures, with larger T L , which generates larger mean dispersion and lower maximum mean concentration. The concentration variance instead was found to approach convergence only for (N x ) ≥ 1024.
Following F&R and Nironi et al. (2015), Fig. 12 shows the intensity of concentration fluctuations defined as max(σ c (x)) max( c(x) ) as a function of along-wind distance from the source. Note that these two maxima are not necessarily in the same position, and that the ratio increases systematically with increasing of the resolution. However, following the results above on mean concentration and variance, it is clear that the changes in the relative intensity of the concentration fluctuations moving from 256 to 512 resolutions and from 512 to 1024 resolutions is dominated by an increase in the fluctuations, i.e. σ c . On the contrary, the increase in intensity from 1024 to 2048 resolutions is dominated by the change (decrease) of mean concentration ( c ) while the standard deviations of fluctuations are very similar between the two resolutions.
Overall, i c does not tend to converge as nicely as the standard deviation when moving from the 1024 to 2048 resolution, with a relative difference of about 30% over most of the distances. However, the fractional variation between 512 and 1024 simulations is clearly more pronounced than between 1024 and 2048 simulations.
The comparison with the wind-tunnel experiments in Fig. 12 shows that the 2048 simulation is in better agreement with the measurements. We note that the source sizes are not identical but our source size lies in between the Nironi et al. (2015) and F&R configurations. A remarkable point is that the relative intensity of fluctuations in the LES seems to have a slower decay compared to the wind-tunnel observations. This indicates that, despite SGS and numerical diffusivity and despite neglecting SGS fluctuations, the LES has less mixing and lower dissipation.  Figure 13b shows that normalizing the crosswind coordinate by the local standard deviation of plume spread, σ y , the curves become similarly spaced in the centreline and towards the plume edges. Figure 13c shows i c for the 2048 simulation, as a function of the crosswind position and for three downwind distances from the source. The vertical dashed lines mark distances of 2σ y from the plume centreline on each side and indicate the position of the plume edges.
The skewness Sk = (c − c ) 3 /σ 3 c and kurtosis K u = (c − c ) 4 /σ 4 c are commonly used indicators of p.d.f. shape. The skewness measures the degree of asymmetry of the p.d.f. and the kurtosis indicates the significance of the tails of the p.d.f.. However, the concentration is bounded at zero, c ≥ 0, and the p.d.f. is positively skewed over most of the downwind distances. Therefore, Sk and K u are both strongly influenced by the right tail of the p.d.f.. Figures 14a, c show the crosswind variation of Sk and K u at an elevation of maximum mean concentration, for x/δ = 1.25 and all grid resolutions, plotted as a function of y/σ y to avoid overlapping due to the differences in mean plume spread. Both Sk and K u strongly increase while moving towards the plume edges where the concentration is more intermittent. Moving towards the plume edges, it seems that the relative difference found in the centreline between grid resolutions disappears between 256 and 512 resolutions while it is conserved between 1024 and 2048 resolutions. However, despite the long averaging time, these highorder statistics are subject to large statistical fluctuations for both 1024 and 2048 resolutions. The statistical fluctuations make the interpretation of the data difficult, especially in the plume edges. However, according to the Gamma p.d.f. (Eq. 17) model Sk = 2i c ≈ 100 for our dataset in the plume edges, and K u = 1.5Sk 2 + 3 ≈ 15000, this seems to suggest that the values reported are not far from reality.
Figures 14b, d show the along-wind variation of Sk and K u on the position of maximum mean concentration. More precisely, the minimum of the statistics over a small area (Δy×Δz) surrounding the position of maximum mean concentration is reported together with the spatially-averaged value over the same area. The area starts with source size and afterward expands up to Δy = Δz = 2 × 0.037δ to account for the plume expansion. This approach simulations, the two estimates are very similar and only the minimum is reported. See text for more details. c, d As in a, b but for the kurtosis was used because the value exactly in the position of maximum was affected by fluctuations that hindered the detection of the along-wind variation. The two estimates are respectively biased towards low values (the minimum over the area) and towards high values (the average over the area). For the average, the bias arises from the crosswind profile strongly increasing away from the plume centre. Together, the two curves provide an upper and lower bound and give a good indication of the behaviour of the statistics along the plume centreline. For the lower resolutions, i.e., the 256 and 512 simulations, only one line (the minimum) is reported because the two estimates are very similar.
On the plume centreline, a significant increase is visible for both Sk and K u while increasing the resolution, and is similar to the behaviour of i c but more pronounced. This means that the weight of the right tail in the p.d.f. increases with the grid resolution. However, a careful inspection of Fig. 14b reveals that, for Sk, the relative difference between the 2048 and 1024 simulations is about 50% and much lower than between 1024 and 512, about 100%. This suggests that also for Sk a clear tendency to convergence is observed between the 1024 and 2048 simulations. It is remarkable that the difference between 256 and 512 is rather small, highlighting again that LES-grid sensitivity studies must span a wide range of resolutions to be meaningful. Understanding the actual behaviour of K u is more difficult since this statistic is affected by the largest statistical errors. Considering only the lower bound of the estimate (the curve given by the minimum in Fig. 14d) a tendency to converge could be again deduced between the 1024 and 2048 simulations, but this is not the case considering the average over the area, i.e. the upper bounds of the shaded area.
Overall, these results show that the behaviour of the right tail of the p.d.f. is even more sensitive to the changes of resolution than is the variance. This was perhaps expected but could be clearly appreciated only when the grid resolution moved from 512 to 1024 and 2048. In contrast, a change between 256 and 512 resolutions showed very limited change in Sk and K u and would erroneously point towards an almost perfect convergence, while both these resolutions are clearly incapable of correctly capturing Sk and K u and therefore the tail of the p.d.f.. A significant point is that the 1024 and 2048 simulations have very similar standard deviations of concentration over most of the downwind distances. Therefore, similar Sk and K u values imply similar third-and fourth-order centred moments. This implies very similar actual concentration values in the p.d.f. tail considering that the differences are greatly enhanced for the third and fourth moment of concentration. Kurtosis especially highlights even quite limited discrepancies in the p.d.f. tail between the 1024 and 2048 resolutions.
Despite the different source sizes in the wind-tunnel experiment and LES, the comparison with the experimental values of Nironi et al. (2015) (see also Marro et al. 2018) shows that the values of both Sk and K u for the 2048 simulation are reasonable and realistic. However, considering the larger source size in the LES, the upper bound of the LES estimates is somewhat higher than the experimental observations while the lower bound is in a good agreement with the observations. As is evident from the changes in i c , Sk and K u, the shape of the p.d.f. significantly changes with the distance from the source, both in the along-wind and crosswind directions. This is true for any grid resolution.
An analytical model of the concentration p.d.f. that in recent years has gained considerable experimental and theoretical support is a family of one parameter Gamma distributions (Duplat and Villermaux 2008;Yee and Skvortsov 2011;Nironi et al. 2015). The concentration p.d.f. is defined as with Γ (κ) the Gamma function, κ = i −2 c and χ = c/ c . For the Gamma p.d.f. Sk = 2i c and K u = 1.5Sk 2 + 3 (e.g., Nironi et al. 2015).
Complementary to the Gamma p.d.f. model is the meandering plume model, extensively discussed above, which fully describes the concentration p.d.f. in the early phase of dispersion where concentration fluctuations are produced. According to the meandering plume model, the centred moments, skewness and kurtosis can be calculated using a standard relationship between centred and uncentred moments (e.g., Monin and Yaglom 1971). This (see Appendix 3, Eq. 28) highlights that close to the source and along the plume centreline (y = y s , z = z s ), scaled moments such as i c , Sk and K u are only a function of M y = σ 2 my /σ 2 r y and M z = σ 2 mz /σ 2 r z . Therefore plotting Sk and K u as a function of i c gives a unique curve. Figures 15a, b show the minumum values of skewness and kurtosis in the small area surrounding the maximum mean concentration for all grid resolutions as a function of i 2 c . The LES data show a clear path, with data moving from the meandering p.d.f. limit, close to the source (dashed line) with low values of i c progressively increasing toward the peak, to the Gamma p.d.f. far away from the source (solid line) with progressively decreasing value of i c after its peak. The peak in i c thus marks the start of the phase where the Gamma p.d.f. is the optimal model. This was never discussed before to our knowledge. We note that the meandering plume model values are only shown for the increasing value of i c up to its peak (see also Fig. 11). Different resolutions seem to behave slightly different with the higher resolutions more closely following both the meandering model close to the source and the Gamma p.d.f. after the i c peak. This confirms that by increasing the resolution, the LES better captures the physics of dispersion. Figure 15c, d report similar comparison but for the off centreline position y = 2σ y . The meandering plume model values almost perfectly collapse on a single curve, therefore only one line is reported in Fig. 15c, d. It overall confirms the behaviour observed on the centreline but noticeably, at constant i c , the difference between values of Sk and K u in the meandering and dissipative phases is strongly reduced.

Summary and Discussion
Four different grid resolutions have been used to perform LES of plume dispersion and concentration fluctuations in an incompressible half-channel flow at infinite Reynolds number. The analysis of the difference in results for changing resolution was supported by a comparison with wind-tunnel experimental data.
For the velocity field, the most salient points of the comparison among resolutions were that, despite the good consistency of one-point second-order velocity statistics, the two-point correlation analysis showed that the highest resolution simulation (2048) developed larger turbulent structures, characterized by longer length scales and time scales. The TKE spectrum revealed that an inertial subrange was present in the two more resolved simulations.
The comparison with the wind-tunnel data showed that the variance of the along-wind velocity component and TKE were well captured, especially at high resolutions, while the variance of crosswind and vertical components (which are more important for plume dispersion) was lower in the LES than in the measurements. The LES showed also a generally higher mean flow speed at source elevation compared to the experimental data of Nironi et al. (2015) and this could be mainly attributed to the specific wall model used in PALM.
Moving to the scalar field, the comparison of the mean concentration among the different resolutions showed that the highest resolution simulation (2048) had a markedly different mean field characterized by higher plume spread and lower mean concentrations. This was tracked back to a larger T L resulting from the larger turbulent structures generated at this grid resolution. As a result, the absolute dispersion obtained at this resolution was more similar to the wind-tunnel measurements despite still being underestimated.
The analysis of concentration variance among different resolutions showed a clear tendency to converge in the two higher resolved simulations (1024 and 2048). The analysis demonstrated that production close to the source is the most critical phase for plume concentration fluctuations and that only an appropriate resolution can capture it, corresponding to N z , N y ≥ 256 and N x ≥ 1024 here. More precisely, the conservation equation of scalar variance showed that the production of fluctuations mainly occured close to the source while afterward advection was the dominant term. Therefore, significant convergence in the variance was ensured by the source being adequately resolved by at least 4 3 grid nodes.
The mechanism of generating fluctuations was analyzed by separating plume dispersion in meandering and relative dispersion and using the Gifford (1959) meandering model for a quantitative analysis. This analysis demonstrated that inadequate resolution artificially enhanced relative dispersion, therefore suppressing meandering related production of fluctuations. As the meandering is often thought to be associated with the large eddies, one could think that a relatively low resolution LES is sufficient to capture it correctly. This intuitive picture is not correct and, as demonstrated in our analysis, numerical diffusion can suppress fluctuations produced by the meandering close to the source if the source itself is not properly resolved. This early meandering phase turns out to be very important for the production of concentration variance. The artificial diffusion imposed by a coarse grid resolution seems to result in a larger effective source size with the corresponding lower level of fluctuations.
The analysis of the shape of the p.d.f. included the examination of scaled moments up to the fourth, i c , Sk, K u. Regarding i c , by comparing all the grid resolutions, the combined behaviour of mean concentration and standard deviation generated a weaker convergence of this scaled moment although a tendency to converge was again visible between the 1024 and 2048 simulations. However, two distinct reasons of this behaviour could be determined: (i) an increase in concentration variance production acting mainly on all resolutions up to 1024, and (ii) a decrease in the mean concentration acting on the 2048 resolution. Therefore, insufficient scalar source definition acted mainly on the 256 and 512 resolutions totally preventing i c convergence at these resolutions, while a change in the turbulent structure modified the mean value of the 2048 simulation.
Regarding Sk and K u, some tendency to converge was again visible between the 1024 and 2048 simulations, especially for Sk. From 256 to 512 to 1024 resolutions, the production of fluctuations close to the source increases significantly. Here, the third and fourth centred moments increase even more than σ c , because the plume retains a progressively thinner initial structure that afterward is transported by mean advection and turbulence; significantly different internal plume structure is present at different resolutions. Moving from 1024 to 2048 resolution the production of concentration variance is very similar and therefore the internal plume structure is similar over most of the concentration values. However the third and, more markedly, the fourth centred moments are mostly influenced by the higher concentration values in the right tail of the p.d.f. and strongly enhance even limited differences in the actual concentrations.
Irrespective of grid resolution, we found that the different p.d.f. shapes could be very nicely reproduced by a family of Gamma distribution in the decaying phase of i c , i.e. downwind of its peak. While the use and observation of the Gamma p.d.f. is not new in the literature of plume dispersion (e.g Yee et al. 1993c;Yee and Skvortsov 2011;Nironi et al. 2015), this very clear transition of the concentration p.d.f. to a family of Gamma p.d.f. after the peak of i c was not observed before to our knowledge. We emphasize that a similar behaviour was observed here in the plume centre and in the plume periphery. The fact that all LES resolutions showed similar behaviour, irrespective of the quite different turbulent velocity spectrum, indicates that the mechanism generating the Gamma p.d.f. is rather general and well captured by the LES.
Complementary to the Gamma p.d.f. model is the concentration p.d.f. generated by the meandering plume model of Gifford (1959), which was shown to quite accurately capture the concentration p.d.f. generated by the LES in the early phase of plume dispersion, before the peak of i c . This is especially the case for the higher resolutions, 1024 and 2048.
We underline that our study highlighted once more that convergence in the turbulent flow of LES is difficult to achieve and this is exacerbated for the scalar field from a localized small source. Different mechanisms drive the production of fluctuations and the long term evolution of mean plume dispersion and therefore multiple distinct factors could prevent convergence of the concentration p.d.f. shape. However, we systematically investigated the differences in the results among grid resolutions and discussed the mechanisms producing these differences. This provides a useful reference to anybody interested in validating and using LES for plume dispersion and concentration fluctuations at very high Reynolds number, and provides a quantification of the uncertainty involved in this type of studies. A very important point is that the grid dependency quantification should span a wide range of grid refinements, in our case the difference in scalar fluctuations between 256 and 512 resolutions was minimal while there was a large change between 512 and 1024 resolutions. Furthermore, the results clearly show that p.d.f. shapes could be very nicely reproduced by a family of Gamma distributions only in the decaying phase of i c , downwind of its peak. The meandering p.d.f. was a better model in the early phase of dispersion.
Detailed analysis and comparison with wind-tunnel measurements showed that, as expected, the 2048 simulation produced the more realistic turbulent structures and therefore concentration p.d.f. evolution. Simulating a wind-tunnel flow means that some important characteristics of the real atmosphere are neglected in the simulations, like for example the momentum flux and TKE shear production at the boundary-layer top. However, the scalar source considered here is placed in the lower part of the boundary layer where these effects should be negligible compared to the interaction with the wall. It remains a task for future studies to elucidate any difference in concentration fluctuations between the current windtunnel setting and a full atmospheric boundary-layer simulation.
Finally, we remark that our results are somewhat specific to our set-up, in particular the SGS model, wall model and numerical methods. Many models are available and the convergence of the model results and the turbulence statistics may depend on the SGS model (e.g., Salesky et al. 2017) and wall model (e.g., Hultmark et al. 2013). In particular, it is likely that the use of suitably formulated Lagrangian stochastic SGS models for the scalar field may improve the convergence of the scalar field results and will be further investigated in future studies. However, the methodology presented here to validate and understand LES results for plume dispersion and concentration fluctuations is general and should be a useful guidance for anyone interested in using LES for such tasks.  conditions for the passive scalar in along-wind (x) and crosswind (y) directions were modified and tested thoroughly. In the inlet plane (the (y, z) plane corresponding to the first grid point in along-wind direction), Dirichlet boundary condition is used. This means that the scalar quantity is set to zero (c = 0), therefore no entrance of scalar is allowed from the inlet side. From the opposite side (outlet), the scalar can only exit. This is done by setting the scalar in ghost layers equal to that of the last grid point in (x), up to which the prognostic equation is solved (∂c/∂ x = 0). This was done according to what PALM used for the scalar in conjunction with non periodic boundary conditions for velocity (outflow). Note that in PALM, some number of ghost layers are considered on each lateral boundary of the processor subdomains. The exchange of ghost layers is adapted to a dynamic number of ghost layers, depending on the applied advection scheme. In the crosswind (lateral) direction (y), the scalar is set to zero at both sides and also the ghost layers before the first grid point in y and after the last grid point in y (for which the prognostic equations are solved). We also used Neumann boundary condition on crosswind direction (∂c/∂ y = 0), however the effect of the lateral boundaries on scalar dispersion statistics was found to be negligible for the range we explored here (x/δ ≤ 4). At the bottom and top, the default PALM boundary conditions of the form of Neumann are used for scalar (∂c/∂z = 0). Boundary conditions of the velocity and scalar fields are detailed in Table 4.
where Q is the emitted mass at the source located in y s , z s . Using the relative dispersion and meandering generated by the LES, all the concentration moments can be calculated according to the meandering plume model. Equation 27 can be rewritten in the form of c n = Q 2π u s σ r y σ r z n 1 1 + nM y (1/2) 1 1 + nM z (1/2) × exp − n(y − y s ) 2 2σ 2 r y (1 + nM y ) exp − n(z − z s ) 2 2σ 2 r z (1 + nM z ) with M y = σ 2 my /σ 2 r y and M z = σ 2 mz /σ 2 r z . Q is 1 kg s −1 , u s is the LES velocity at source elevation, σ r and σ m are relative dispersion and meandering, respectively as derived from the LES concentration data for crosswind and vertical directions, no parametrization is used. The contribution of internal concentration fluctuations, in the coordinate system relative to the local plume centre of mass, is neglected in this standard meandering formulation. In the simplified case with M y = M z = M and in the plume centreline, Yee and Wilson (2000) demonstrated that Therefore, with a constant absolute dispersion the variance of the concentration fluctuations will always increase with M.