Large-eddy simulation and stochastic modeling of Lagrangian particles for footprint determination in the stable boundary layer

Large-eddy simulation (LES) and Lagrangian stochastic modeling of passive particle dispersion were applied to the scalar flux footprint determination in the stable atmospheric boundary layer. The sensitivity of the LES results to the spatial resolution and to the parameterizations of small-scale turbulence was investigated. It was shown that the resolved and partially resolved (“subfilter-scale”) eddies are mainly responsible for particle dispersion in LES, implying that substantial improvement may be achieved by using recovering of small-scale velocity fluctuations. In LES with the explicit filtering, this recovering consists of the application of the known inverse filter operator. The footprint functions obtained in LES were compared with the functions calculated with the use of first-order single-particle Lagrangian stochastic models (LSMs) and zeroth-order Lagrangian stochastic models – the random displacement models (RDMs). According to the presented LES, the source area and footprints in the stable boundary layer can be substantially more extended than those predicted by the modern LSMs.


Introduction
Micrometeorological measurements of vertical turbulent scalar fluxes in the atmospheric boundary layer (ABL) are usually carried out at altitudes z M ≥ 1.5 m due to technological limitations of the eddy covariance method.The measurement results are often attributed to the exchange of heat, moisture and gases at the surface.This procedure is not justi-fied for inhomogeneous surfaces because of a large area contributing to the flux, and because of variability of the second moments with height.The relationship between the surface flux F s (x, y, 0) and the flux F s (x M , y M , z M ), measured in point x M = (x M , y M , z M ), can be formalized via the footprint function f s : f s (x, y, x M , y M , z M )F s (x, y, 0)dxdy. (1) Traditionally, footprint functions f d s (x d , y d , x M ) = f s (x, y, x M ) are expressed in a local coordinate system with the origin that coincides with the sensor position (here, x d = x M − x is the positive upwind distance from the sensor and y d = y M − y is the crosswind distance; see Fig. 1a).In a horizontally homogenous case these functions do not depend on x M and y M .In the ABL the surface area contributing to the flux is elongated in the wind direction; therefore, the crosswind-integrated footprint function f y s defined as is one of the most required characteristics for the practical use.
The measurements of the scalar flux footprint functions in natural environment are restricted (e.g., Finn et al., 1996;Leclerc et al., 1997Leclerc et al., , 2003;;Nicolini et al., 2015) due to the necessity to conduct the emission and detection of artificial tracers.Besides, such measurements are not available for the stably stratified ABL, where the area of the surface influencing the point of measurements increases.
Modeling approaches used for footprint calculation include stochastic models, such as single-particle first-order Lagrangian stochastic models based on the generalized Langevin equation (LSMs) and zeroth-order stochastic models (also known as the random displacement models, RDMs) (see the reviews listed in the papers (Wilson and Sawford, 1996;Wilson, 2015) and the monograph (Thomson and Wilson, 2013)).Besides, one can use the analytical models (e.g., Horst and Weil, 1992;Kormann and Meixner, 2001) and the parameterizations based on the scaling approach (Kljun et al., 2004(Kljun et al., , 2015)).All of these models should be calibrated against the data considered to be representative of real processes.Their results depend on the choice of universal functions in the ABL or in the surface layer (non-dimensional velocity and scalar gradients, non-dimensional dissipation, dispersion of the velocity components, etc.).Commonly, the applicability of the analytical models is limited by a "constant flux layer" simplification, assuming that the measurement height z M is much less than the thickness of the ABL z i .However, under the strongly stable stratification the thickness z i may be several meters; therefore, the vertical gradients of momen-tum and scalar fluxes near the surface can be large.It can lead to incorrect functioning of the models designed for and tested on the data gathered under different conditions.
Large-eddy simulation (LES), employing the Eulerian approach for the transport of scalars, was applied for the first time for a footprint calculation in Leclerc et al. (1997).Modern computational technologies allow one to combine Eulerian and Lagrangian methods for turbulence simulation and particle transport (e.g., Weil et al., 2004;Steinfeld et al., 2008;Cai et al., 2010;Hellsten et al., 2015) and to perform detailed calculations of averaged two-dimensional footprints under different types of stratifications in the ABL and footprints f s (x, y, x M ) over heterogeneous surfaces (for example, urban surfaces and surfaces with alternating types of vegetation).Some examples of such calculations are given in Steinfeld et al. (2008) and Hellsten et al. (2015).
Lagrangian transport in LES is complicated by the problem of the description of small-scale (unresolved) fluctuations of the particle velocity, which is similar to the problem of subgrid modeling of Eulerian dynamics.A common approach for Lagrangian subgrid modeling in LES is the application of subgrid LSMs (e.g., Weil et al., 2004;Steinfeld et al., 2008;Cai et al., 2010;Shotorban and Mashayek, 2006).This approach requires a number of additional calculations for each particle (e.g., interpolations of subfilter stresses τ ij and subgrid dissipation into particle position x p ). Besides, it is necessary to generate a three-component random noise for each particle, which is a time-consuming computational operation.A numerically stable solution to the generalized Langevin equation (see Sect. 2.3,Eq. 9) in LES requires smaller time steps than the steps to solution of Eulerian equations, because local Lagrangian decorrelation time T L (x p , t) can be very small.
The statistics of simulated turbulence in LES may significantly differ from the statistics of real turbulence.For example, the use of dissipative numerical schemes or low-order finite-difference schemes usually results in a suppression of fluctuations over almost the entire resolved spectral ranges of discrete models (see, e.g., Fig. 16 in Piotrowski et al., 2009).Turbulent fluxes (in the Eulerian representation) associated with these fluctuations are restored by subgrid closure.However, in terms of the Lagrangian transport the effects of distortion of the small-scale part of the spectrum are most often not considered.
Numerical simulations of Lagrangian transport in LES are also limited by the low scalability of parallel algorithms.This is due to the impossibility of uniform loading of processors in a joint solution to the Euler and Lagrangian equations, a large number of interprocessor exchanges and an unstructured distribution of characteristics required for Lagrangian advection in the computer RAM memory.
Thus, all methods of numerical and analytical determination of the functions f s have individual drawbacks.At the same time, due to the lack of a sufficient amount of experi-mental data and due to their low accuracy, there are no clear criteria for evaluation of different models.
According to the need for computational cost reduction, one of the objectives of this study is to establish the role of stochastic subgrid modeling in the correct description of the particle dispersion in LES.Is it possible to simplify the calculation and to avoid the introduction of stochastic terms without the loss of accuracy in some integral characteristics, such as the footprints or the concentration of pollutants emitted from the point sources?The role of subgrid fluctuations is reduced with an increase in spatial LES resolution.Therefore, the independence of results from the mesh size is used as a criterion for checking the quality of Lagrangian transport procedures in LES.It will be demonstrated that the subgrid stochastic modeling in LES can be omitted in most cases.Instead, we propose a "computationally cheap" procedure of inverse filtering supplemented by divergent correction of Eulerian velocity to replace the subgrid stochastic modeling in LES (see the description below).
Subgrid transport is especially significant near the surface and/or under stable stratification -all are the cases associated with small eddy size.That is why the stable ABL was selected as the key test scenario in this study.We slightly modified the setup of the GABLS (Beare et al., 2006) numerical experiment for this purpose.
LES results are used as the input data for the stochastic models (LSMs and RDMs).These data are pre-adjusted using known universal dependencies and taking into account an incomplete representation of turbulent energy in LES.The comparison of results of different stochastic models and the results from LES allows one to specify the parameters for the LSMs and permits one to identify the differences between LSMs and RDMs under the conditions that have not been tested previously.
The paper is organized as follows.Section 2 contains the description of some common features of approaches: the implemented numerical algorithm for footprint estimation in the LES and LS models (Sect.2.1); LES governing equations and the definitions of some terminology used for the small-scale modeling description and for the testing of particle transport (Sect.2.2); the definitions of stochastic models (LSMs and RDMs) and pointing to some problems connected with uncertainty in the choice of turbulent statistics for them (Sects.2.3 and 2.4).Section 3 contains a short description of the numerical algorithms, the turbulent closure for the LES model used in this study (Sect.3.1) and the description of the different approaches for the Lagrangian particle transport in the LES tested here (Sect.3.2).Section 4 is mainly devoted to the testing of the ability of the LES model with rough spatial resolution to reproduce particle dispersion correctly.To this end, we implemented a special setup of the numerical experiment (see Sect. 4.1), permitting one to compare Lagrangian and Eulerian statistics (see Sect. 4.2.2).The focus was placed on the approaches with the limited use of subgrid stochastic modeling (see Sect. 4.2.1,where the sensitivity of the computed footprints to the spatial resolution was investigated).The footprints computed with the LES model with simple subgrid LSMs and RDMs (traditional approach) are presented in Sects. 4.2.3 and 4.2.4.Two-dimensional footprints are shown in Sect.4.3.Due to high sensitivity of LSMs to the turbulent statistics, we emphasize data preparation for them using LES results, measurement data and similarity laws in Sect.5.1.Section 5 contains the results of footprint modeling with the use of the set of different RDMs and LSMs (specified in Sect.5.2) in comparison with LES results (see Sect. 5.3).Section 6 summarizes the results.
In addition to the basic calculation, we carried out a series of tests (see Supplement Sect.S1) under unstable stratification in the ABL with different grid steps in the LES model.This allows us to compare the results presented here with similar results obtained in previous studies (e.g., Steinfeld et al., 2008;Weil et al., 2004) and to verify the performance of our LES model in footprint evaluation.Furthermore, we demonstrate the results of footprint calculations above the inhomogeneous surface (Supplement Sect.S2) with a huge number of particles involved in calculations simultaneously.Computational aspects of technology are discussed as well.

Numerical evaluation of footprints
Computational methods for determination of footprints often reduce to the implementation of Lagrangian transport of marked particles.Each particle can contain a number of attributes, including its initial coordinate x p 0 and time t p 0 .Choose two small horizontal plates δ S and δ M for averaging in the neighborhood of zero with the areas S S and S M , respectively.Define the time interval T p = [t 0 , t 2 ], during which new particles are ejected near the ground with the intensity H (here H is the mathematical expectation of the new particle number emitted per unit area per unit time) and the interval T a = [t 1 , t 2 ] (t 1 > t 0 ), when particles are detected near the point of measurement.If t 1 is sufficiently large for the ensemble averaged flux to attain constant value in time, and T a is quite large for statistically significant averaging, then the footprint f s can be evaluated by the formula where n SM is the number of intersections of the plane z = z M by the particle trajectories at horizontal coordinates x  3) and the description above, the particle crossing the test area δ M brings the impact into the value f s (x S , y S , x M ) if the beginning of its trajectory belongs to the test area δ S after trajectory modification such that the point x p 1 coincides with sensor position x M .For example (see Fig. 1b), when the footprint value is calculated at the point (x S , y S ), only the red particle is counted, but not the blue particle.Such an algorithm of averaging was selected because it permits one to refine the footprint resolution in the vicinity of the sensor independently of the area of δ M using the assumption of some spatial homogeneity.
In the horizontally homogeneous case, one can calculate the footprint f d s (x d , y d , z M ) by performing averaging over statistically equivalent coordinates of the sensor position.For this averaging in LES with a periodic domain, one can prescribe the coordinates (x M , y M ) to the domain center and select the area δ S to be equal to the whole domain size.Analogical methods can be applied when using LSMs or RDMs, whereas in the case of RDMs, particle displacement should be used in Eq. ( 3) instead of velocity.
The nonuniform Cartesian grid x d ij = (x d i , y d j ) (where −20 ≤ i ≤ 160; −120 ≤ j ≤ 120), stretched with the distance from the sensor position, was selected for the footprint function accumulation in the following sections of this paper.The grid was prescribed as (x d 0 , y d 0 ) = (0, 0); y j/|j | if i = 0 and j = 0; x0 = y0 = 2 m; and γ x = γ y = 1.05.This grid is independent of the LES model resolution and coincides with the footprint grids selected for all runs with LSMs and RDMs.

Lagrangian particles embedded in LES
Lagrangian particle velocity u p and the particle position x p can be computed in LES models as follows: i is the interpolation of the resolved Eulerian velocity into the particle position; u p i are the small-scale unresolved Lagrangian velocity fluctuations associated with Eulerian velocity fluctuations belonging to "subgrid" and "subfilter" scales.Here and later we shall use the designation "subfilter" to denote the fluctuations that belong to the resolved spectral range of the discrete model, but are not reproduced numerically, and the designation "subgrid" for the fluctuations, which can not be represented on the grid due to smallness of the scales.LES governing equations for filtered velocity u are where F e i comprises Coriolis and buoyancy forces, p is normalized pressure and τ ij = u i u j − u i u j denotes the modeled "subgrid/subfilter" stress tensor.A system of equations ( 5) can be supplemented by the Eulerian equations of scalar transport: where Q s denotes source intensity; ϑ s i = su i − u i s are the parameterized "subgrid/subfilter" fluxes.Usually, the fluctuations u p are defined as dependent on some random function ξ , introduced in order to provide the missing part of mixing.The particular approaches for computing the unresolved part of particle velocity will be discussed and tested in the following sections.
There is a great practical interest in the calculation of footprints, as well as of spatial and temporal characteristics of pollution transport from localized sources above heterogeneous surfaces and in the areas with complex geometry (in the urban environment, over the surfaces with complex terrain or over the alternating types of vegetation).LES of such flows becomes a routine procedure with increasing performance of computers.However, the calculation of statistical characteristics of Lagrangian trajectories is complicated in this case by the need of transport of huge number of tracers (e.g., Hellsten et al., 2015).For example, it is necessary to calculate the trajectories of about 10 9 particles (see Supplement Sect.S2) to obtain the footprints above the inhomogeneous surface with the explicitly prescribed obstacles (the task similar to that presented in Glazunov and Stepanenko, 2015).
On the other hand, a large number of particles (see, e.g., Supplement Fig. S2.1b) allows one to estimate the local instantaneous spatially filtered concentration of the scalar: where G is the function that coincides with the convolution kernel of the LES filter operator and N is the total number of particles in the domain.If the mathematical expectation Q p of a number of new particles ejected in a unit volume during the unit time interval is proportional to the Eulerian concentration source strength Q p (x, t) = CQ s (x, t), then s P (x, t) ≈ Cs(x, t).One can perform the same operations with the "Lagrangian" concentration s P (x, t) as the operations with the Eulerian scalar s.Below, we will compare the averaged values of s P and s and their spatial variability.Besides, we will use the estimation of concentration s P (x, t) to correct the particle velocities (see Sect. 3.2.1,Eqs. 34 and 35), in order to approximate the effect of subgrid turbulence.

Single-particle first-order Lagrangian stochastic models (LSMs)
Another approach (more widespread due to a lower computational cost) is the replacement of the entire turbulent component of velocity by a random process (Lagrangian stochastic models (LSMs)): Here u (p) i is the ensemble-averaged Eulerian velocity at point x p .Note that LSMs are assumed to also be applicable under the temporal evolution of turbulence statistics.In this paper we shall consider the ABL as it approaches a quasisteady state.Therefore, due to the assumption of ergodicity, ensemble averaging can be replaced by averaging in time and in the directions of spatial homogeneity: ϕ ≈ ϕ x,y,t .
A single-particle first-order LSM is formulated as follows.Velocity u p i is described by the stochastic differential equation: where ξ stays for the delta-correlated (usually Gaussian) random noise with the variance dt and with the zero average ξ p i = 0; a i and b ij are the functions depending on the Eulerian characteristics of turbulence and on the Lagrangian velocity of the particle.Typically b ij is calculated by the formula where denotes the energy dissipation rate, averaged for a fixed coordinate, and C 0 is the Kolmogorov constant.This kind of random term (arguments are given in Thomson, 1987 andSawford, 1993) is defined by a Lagrangian velocity structure function in the inertial range (see Monin and Yaglom, 1975) if τ η t T E (here, τ η = (ν/ ) 1/2 is the Kolmogorov microscale, T E = E 2 / is the energy-containing turbulent timescale and E is the turbulent kinetic energy).
The function a i (drift term) determines the behavior of particles at large times t ∼ T L ∼ T E (here T L is the Lagrangian decorrelation timescale).For spatially inhomogeneous and statistically non-stationary turbulent flows, including the ABL, the choice of a i is usually made according to the well-mixed condition (WMC; Thomson, 1987).In general WMC does not lead to a unique solution for a i .Different LSMs are constructed by introducing the additional physical assumptions, and can lead to inequivalent results.
Lagrangian models are very sensitive to the choice of universal functions that define the normalized root mean square (RMS) of the vertical velocity σ w = w 2 1/2 /U * and nondimensional dissipation = z/U 3 * (here U * is the friction velocity).Besides, the simulation results are affected by the choice of value of C 0 .It can be shown (e.g., Durbin, 1984;Wilson and Yee, 2007) that for one-dimensional LSM, these parameters determine the eddy diffusivity K s for the scalar in the diffusion limit (when t T L , i.e., at large distances from the source): The data of measurements in the ABL demonstrate large variation.For example, the values of σ 2 w range from 1.0 to 3.1 (see Table 1 in Banta et al., 2006).According to Eq. ( 13) it implies the change in K s by more than 9 times.
There is no consensus on the value of C 0 either.Formally, C 0 has the meaning of a universal Kolmogorov constant in Eq. ( 11).The estimation of this constant for an isotropic turbulence using the data of laboratory measurements and DNS provides an interval C 0 = 6.± 0.5 (see Lien and D'Asaro, 2002).However, the values C 0 ∼ 3-4 are often used for the LSM of particle transport in the ABL, independently of the type of stratification.These values have been obtained by the different methods.For instance, the value C 0 = 3.1 for a onedimensional LSM corresponds to a calibration performed in Wilson et al. (1981) according to observation data (Barad, 1958;Haugen, 1959).This calibration (see Wilson, 2015) assumes that the turbulent Schmidt number Sc = K m /K s = 0.64 is near the surface (here K m is the eddy viscosity).It is known that determination of the turbulent Prandtl number P r = K m /K h (K h -heat transfer eddy diffusivity) and the Schmidt number based on observation data is complicated by large statistical errors associated with the problem of self-correlation (Anderson, 2009;Grachev et al., 2007).Therefore, this method of estimation of C 0 can not be considered final, and should be confirmed by future studies.In Rizza et al. (2010) the values of C 0 were determined using the LES-based evaluations of the velocity structure functions and the Lagrangian spectra in convective and neutrally stratified ABLs.In this study the LES model had a relatively low resolution, which can be insufficient for accurate determination of this constant in the inertial subrange (see the discussion on the resolution requirements in Lien and D 'Asaro, 2002).Nevertheless, the value C 0 ∼ 3 in the paper by Rizza et al. (2010) is relevant for LSMs applied to the convective ABL; in that case, the value of C 0 is also responsible for the energy containing timescales that are well resolved in LES.The detailed overview of the methods of determination of the constant C 0 can be found in Poggi et al. (2008), where the discussion on the disagreements of the different approaches is also included.The results of the LSMs are very sensitive to the choice for C 0 , as was shown earlier by Du et al. (1995), Rotach et al. (1996), Wilson (2015) and many others.Below we show that the commonly used value of C 0 ∼ 3-4 can be greatly underestimated for use as a parameter in LSMs applied to the stably stratified ABL.
2.4 Zeroth-order Lagrangian stochastic models or random displacement models (RDMs) The simplest approach for development of the models of particle dispersion entails replacement of the Eulerian advection-diffusion equation by the stochastic equation for particle position (random displacement models -RDMs): Probability density of particle position P is connected with scalar field concentration s as follows: Using the Fokker-Planck equation, it can be shown that Eq. ( 15) is equivalent to Eq. ( 14) from the point of view of concentration transport when the time step dt tends to zero (Durbin, 1983;Boughton et al., 1987).
An RDM has some major disadvantages.First, it shares the limitation of Eulerian eddy-diffusion treatment of turbulent dispersion, i.e., "K-theory".Correspondingly, it is not able to describe the non-diffusive near field of a source.Also, an RDM can not be applied for the convective ABL, where the counter-gradient transport is observed.Besides, it requires the exact values of diffusion coefficient K s , which can not be measured directly.

Numerical algorithms and turbulent closure
A system of equations (5-6) is discretized using an explicit finite-difference scheme with the second-order temporal approximation (Adams-Bashforth method) and fourth-order (fully conserved for advective terms) spatial approximation of velocity and scalars on a staggered grid (Morinishi et al., 1998).
A mixed model (Bardina et al., 1980), expressed as the sum of the Smagorinsky and scale-similarity models, is used for calculation of the turbulent stress tensor: where S ij is the filtered strain rate tensor, and C s is the dynamically determined (Germano et al., 1991) dimensionless coefficient that depends on time and spatial coordinates.The a priori tests using the data of laboratory measurements show that scale-similarity models with Gaussian or box filters provide correlation typically as high as 80 % between real and modeled stresses (see the overview in Meneveau and Katz, 2000).The significant part of this correlation can be attributed to non-ideality of the spatial filter and use of common information for computing both the real and modeled stresses (Liu et al., 1994).The discrete spatial filter used in this study has a smooth transfer function in spectral space, so it can be supposed that the scale-similarity part of Eq. ( 17) is mainly responsible for the influence of velocity fluctuations belonging to "subfilter" scales.
The procedure of calculation of the coefficients X(x, t) = (C s ) 2 reduces to minimization of the functional , where is the model domain and ε ij (x) is the residual of the overdefined system of equations obtained by substitution of the mixed model (Eq.17) into the Germano identity as Here T ij are subgrid/subfilter stresses for the smoothed velocity u, obtained by successive application of basic F and test F spatial filters; α = / is the ratio of the filter widths.Tensors M T ij , M τ ij , L ij and H ij are calculated as follows: The generalized solution to the discrete analog of Eq. ( 18) is searched using the iterative conjugate gradients (CG) method with a diagonal preconditioner.To do this, the problem is reduced to a linear system of equations where X is the desired solution (a vector of dimension N = N x N y N z with the values defined in the center of the grid cells); A and R = L −H are the discrete analogs of the operator and the right-hand side of Eq. ( 18) correspondingly; A * is the transpose matrix.The diagonal preconditioner P for the CG method was selected as follows: where µ = const ∼ 1 is the empirical coefficient independent on time and spatial position.The solution X contains negative values (unconditional minimization of the functional is used), however, mixed model (Eq.17) reduces their relative number compared with the dynamic Smagorinsky model.
In the algorithm, negative values are replaced by zeroes.In fact, this dynamic procedure is close to approach proposed in Ghosal et al. (1995), with the difference that the mixed model was applied here and iterative method was replaced by a faster CG method.
Eddy-diffusion models are used for subgrid heat and concentration transfer: |S| is the eddy diffusivity, which is independent of the type of scalar.Subgrid turbulent Schmidt and Prandtl numbers are fixed: Sc subgr = P r subgr = 0.8.
A distinctive feature of this model is that the discrete spatial filter operator F = F x F y F z is explicitly involved in calculation of stresses.The following discrete basic filter is selected: here, i, j, k denote a grid cell number.ϕ is any variable.Similar filtering is applied along the coordinates y and z.It is reasonable to expect that we get the velocity u, smoothed according to the specified filtering operator as a solution to Eq. ( 5) supplemented by the mixed closure (Eqs.17-21).
Since the discrete filtering operator is invertible, we can find the following velocity at any point and time: which better reflects the small-scale spatial variability.The approximate inverse filter is calculated as a series (Van Cittert, 1931) where I is a unity operator; in the calculations presented below we used n = 5.Spatial spectra of "defiltered" velocity u * under the neutral, unstable and stable stratifications were obtained earlier (Glazunov, 2009;Glazunov and Dymnikov, 2013;Glazunov, 2014b).It was found in all cases that this procedure improves the small-scale parts of the spectra according to dependence S ∼ k −5/3 , provides better agreement of spectra calculated with the different spatial resolution, and improves the convergence of non-dimensional spectra if proper length scales are used for normalization.
3.2 Methods for Lagrangian particle transport in LES

Subgrid and subfilter modeling
Below, the subgrid and subfilter modeling methods used for the simulations in the current study are listed.These methods will also be used in combinations as defined in Sect.4.2.
(1) Improvement of Lagrangian transport using inverse filtering of Eulerian velocity field First, we will use the recovering of "subfilter" fluctuations (Eqs.25 and 26) in order to transport Lagrangian particles more precisely: Note that for the use of such a procedure, LES models should exhibit the properties of a model with an explicit filtering.A similar approach was recently applied by Michalek et al. (2013) in LES with an approximate deconvolution subgrid model (ADM; see Stolz et al., 2001), which can also be considered as the model with explicit filtering.In most cases, the suppression of small-scale fluctuations in LES (particularly in those that use a low-order numerical scheme) occurs as a result of the combined effect of approximation errors and the subgrid closure.Therefore, the shapes of effective spatial filters of most models can only be determined by a posteriori analysis of the calculation results.
(2) Lagrangian stochastic subgrid/subfilter model Second, we will apply the subgrid stochastic model proposed in Shotorban and Mashayek (2006): The parameter C 0 was specified to be equal to 6, because the stochastic part of the model (Eq.28) is mainly responsible for spatial scales and timescales in an isotropic inertial subrange of the turbulence.When using the dynamic mixed model (Eqs.17-21), a value of is not calculated directly, and then it is assumed that the dissipation is locally balanced by shear production and buoyancy production or sink.In addition, since this model can produce a local generation of kinetic energy, the averaging in a horizontal plane was performed to avoid negative values of dissipation: where ϑ 3 is the vertical subgrid flux of potential temperature and g/ 0 is the buoyancy parameter.Timescale T L was evaluated as www.geosci-model-dev.net/9/2925/2016/Geosci.Model Dev., 9, 2925-2949, 2016 Thus, the total unresolved kinetic energy was calculated as the sum of "subfilter" energy and "subgrid" energy: To evaluate the value E subgr it was supposed that "subgrid" fluctuations belong to quite a wide inertial range with the component-wise velocity spectra , and that the minimal wavenumbers for these fluctuations kmin i = π/ gi correspond to wavelengths in two grid steps.Here, gi is the grid step in the appropriate direction and C K = 18 55 C K = 0.5 is the Kolmogorov constant (here, C K ≈ 1.5 is the Kolmogorov constant associated with three-dimensional wavenumbers).
All the values required for a application of this model were linearly interpolated into the particle position everywhere except at heights z < g /2, where we use the constant values T L ( g /2) and ( g /2).This procedure is rather arbitrary, but it does not have large impact on the results due to the small decorrelation time T L ( g /2).Besides, there are no physically grounded reasons for the justification of such interpolations in LES because the resolved velocity in the vicinity of surface is greatly corrupted by the approximation errors.Such procedures should be considered as an adjustments depending on the numerical scheme and on the subgrid closure.
(3) Random displacement subgrid/subfilter model Third, the RDM specified in Sect.2.4 will be adopted for the Lagrangian particles subgrid dispersion.In this case we shall use the same subgrid diffusivity K h subgr both for the Eulerian scalars (Eq.23) and for the particles displacement calculations: This model does not contains the arbitrary specified parameters except those which were already used in the Eulerian LES.The coefficient K subgr s was linearly interpolated into the particle positions at heights z ≥ z 0 with the assumption that K (4) Divergent correction of the Eulerian velocity field Finally, in order to find out whether the subgrid mixing is one of the key processes in the dispersion of Lagrangian tracers, we introduced an additional correction to the particle velocities: where u div is the deterministic divergent additive to the velocity field u: with the imposed restriction u div,i = 0 if s P = 0. Here, the "subgrid" flux ϑ sp i is calculated using the same closure as the closure for Eulerian scalars s, with the only difference that the concentration s P , estimated by the number of particles in a grid cell, is used in Eq. ( 23).The applicability of this procedure justified because of the large number of particles involved in simulation (in all the cases described below we have at least several dozens of particles in each grid cell).
Correction given by Eqs. ( 34), ( 35) does not provide true small-scale mixing, but only introduces an additional "stretching" or "compression" of the small volumes filled with particles and provides concentration fluxes across the borders of grid cells close to "subgrid" fluxes in Eulerian model.Using this correction, we are guaranteed to get a high correlation between the "Eulerian" and "Lagrangian" concentrations (in all our preliminary tests s s p xy s 2 s 2 p ≈ 0.9).The idea of such a correction was based on the assumption that details of the mechanism of subgrid mixing have a little influence on the statistics of trajectories at sufficiently large distances from the source and at long enough time t.It was assumed that the quick mixing on small spatial scales can be implicitly substituted by the approximation errors arising in the procedures of interpolation and by the errors of discrete solution to the advection equation.Correction brings an additional systematic effect to reduce incorrect particle transport by the large eddies.

Simplified velocity interpolation
In preliminary tests it became clear that trilinear interpolation of each velocity component provides no advantages for footprint calculation in comparison with the following simplified linear interpolation on a staggered grid: where position (i, j, k) is the center of a grid cell containing the particle.Trilinear interpolation and interpolation given by Eq. ( 36) provide nearly the same concentration fluxes across the borders of a grid cell, but the latter does not result in additional substantial smoothing of velocity.An exception was made for the grid layer closest to the surface (z p < g ) where the mean velocity components were adjusted according to the Monin-Obukhov similarity theory with the dimensionless functions taken from Businger et al. (1971).The calculations were performed at the equidistant grids with steps g = 2.0, 3.125, 6.25 and 12.5 m.The size of the horizontally periodic computational domain was equal to 400 × 400 × 400 m 3 .The last hour of numerical experiments was used for averaging the results and subsequent analysis.This setup is based on the observation data (see Kosović and Curry, 2000).As was shown in Beare et al. (2006), the LES results obtained under the same conditions with the different models converged with the higher grid resolutions.Later, this case was used for testing the LES models, e.g., in Zhou and Chow (2012) and Bhaganagar and Debnath (2015) and many others, and for the improvement of subgrid modeling, e.g., in Basu and Porté-Agel (2006), Zhou andChow (2011), andKitamura (2010).The LES model presented here was tested earlier under the non-modified setup of GABLS in Glazunov (2014a), where the turbulent statistics above a flat surface and above an urban-like surface were investigated.In all of these studies, LES results were in agreement with the known similarity relationships for the stable ABL.This allows one to consider the LES data for GABLS as a reference case for testing of the approaches utilizing the statistical averaging of the turbulence (e.g., see Cuxart et al., 2006, where the intercomparison of single-column models was performed).Several of the non-dimensional relationships in the stable ABL were collected and presented in Zilitinkevich et al. (2013).The considered case is also included in the LES database for this study and fits well with the different stability regimes after the appropriate normalization.Therefore, the results obtained in this particular case can be generalized for many cases due to the similarity of the stable ABLs.Besides, the presented simulations are easily reproducible, and they can be repeated using any LES model that contains the Lagrangian particle transport routines.
The mean wind velocity and the potential temperature, calculated with the different spatial steps g , are shown in Fig. 2. The model slightly overestimates the height of the boundary layer at coarse grids; however, the wind velocity near the surface is approximately the same in all runs.As one can see from Fig. 2, the results of the simulation are in good agreement with the results from other LES presented in Beare et al. (2006) (see http://gabls.metoffice.com for more information).The mean wind profile computed in accordance with Högström (1996) is shown in Fig. 2 by the vertical dashes; in the surface layer part of the domain this "standard" profile for the stable conditions almost coincides with the longitudinal velocity obtained in LES.
Passive Lagrangian tracers were transported simultaneously with the calculations of dynamics.Each particle, when reaching a lateral boundary of domain, is returned from the opposite boundary in accordance with periodic conditions.The reflection condition is used at the ground.The particles are ejected at the height z 0 = 0.1 m (one particle per each grid cell adjacent to surface) with regular time intervals t ej = 1 s.The position of the new particle within a grid cell is set randomly with uniform probability.The ejection of particles takes place continuously from the seventh to the ninth hour of the experiment.
To limit the number of particles involved in the calculation, the absorption condition is applied at the height of 100 m within the ABL.It was verified previously that the upper boundary condition does not have a large impact on the results of calculations of footprints for the heights z M up to 60 m and for the distances x − x M considered in this paper (see Appendix A and the test with the LSM shown by the orange curves in Fig. 11a, c, e).This formulation of the numerical experiment allows direct comparison of the concentration of particles s P , estimated by Eq. ( 7), and the scalar concentration s, calculated by the Eulerian approach (Eq.6).For this purpose, an additional scalar s is calculated from the seventh till the ninth hour, with a constant surface flux F s = const = 1, zero initial condition and the Dirichlet condition s = 0 at altitude 100 m.
In the last hour of simulation, the averaged number of particles in each cell of the grid near the surface was approximately equal to 700-800, 350-400, 180-200 and 110-130 for grid steps g = 12.5, 6.25, 3.125 and 2.0 m, respectively.Having such a number of particles, one can estimate the concentration s p (x i,j,k , t m ) at each time step, where x i,j,k is the center of a grid cell.It was assumed that each particle contributes to the concentration s P (x i,j,k ) with the weight r p i,j,k = (V p V i,j,k )/V i,j,k , where V p is the rectangular neighborhood of its position with the side g , (V p V i,j,k ) is the volume of intersection with a grid cell, and V i,j,k is www.geosci-model-dev.net/9/2925/2016/Geosci.Model Dev., 9, 2925-2949, 2016 the cell volume.This averaging is close to the filtering of an Eulerian scalar (Eq.24).The additional normalization is performed as follows: s P = s P t ej / z .The concentration s P corresponds to the number of particles in one cubic meter under the condition that one particle per square meter per second is ejected near the surface.Concentration s P is numerically equal (excluding errors, determined by different methods of transport) to the concentration of the scalar field s if scalar surface flux F s = 1.
Figure 3 shows the resolved and the parameterized components of flux w s in runs with different grid steps.It is seen that the calculation time is not large enough to reach a steady state (the total flux is not constant with the hight, so the average concentration continues to grow during the last hour).However, it was checked that the flux footprint close to the sensor is not affected by nonstationarity.Besides, we can compare the values of s and s P , because the boundary and initial conditions are identical for them.
The unresolved fraction of the flux F sbg s = ϑ s 3 is an essential part of the total flux F tot s = s w + ϑ s 3 .Accordingly, the vertical transport of Lagrangian particles by resolved velocity u may be significantly underestimated.Thus, we have a "hard" enough test to verify Lagrangian transport in LES with a poorly resolved velocity field.

Footprint calculation with limited application of subgrid stochastic modeling in LES
Figure 4 shows the scalar flux footprints averaged in crosswind direction f y s (x M − x, z M ) computed by different methods and with different grid steps.In all cases, we have avoided using the subgrid-scale stochastic modeling, except for calculating the velocity of the particles located within the first grid layer z p < g .For the curves marked "st_1l", the resultant velocity of the particles near the surface was calcu-lated as follows: where the function r(z p ) is defined as r(z p ) = (1 − z p / g ) if z p < g and r(z p ) = 0 if z p ≥ g ; u p is the random velocity component, calculated using the stochastic subgrid model (Eq.28).To take into account the memory effects in Langevin equation, the stochastic model was implemented inside the layer z p < 3 g , so (because of the smallness of scale T L ) this procedure does not lead to significant distortions in the random component of the velocity.
If the particles are advected by the filtered velocity u without any correction then the vertical mixing is too weak and the maxima of footprints f y s are strongly underestimated and shifted at the large distances from the sensor position.Divergent correction of Eulerian velocity (Eqs.34, 35) partially improves the results (squares in Fig. 4a, b).For example, maximum of footprint f y s for the sensor height z M = 30 m (near the fifth computational level) occurs to be close to the maxima of footprints, computed at fine grids, but it is still shifted.Thus, the correction (Eqs.34, 35) alone is not sufficient.Primarily this is due to the weak mixing below the first computational level, where the contribution of the subgrid velocity is crucial.
The inclusion of stochastics within the first layer improves the result (dashed curves in Fig. 4a, b).However, it is not enough to determine footprints at altitudes comparable to the grid spacing.
The advection of particles by the velocity u * leads to close matching of functions f y s , calculated with different grid steps (solid lines of different thickness in Fig. 4c, d).The differences between these footprints are not significant from a practical point of view, and can be equally explained by means of the incorrect Lagrangian particles transport, as well as by means of the insufficiently accurate solution to the Eulerian equations on the coarse grid.

Spatial variability of scalar concentration inferred by Eulerian and Lagrangian methods
While the particles were advected by the "defiltered" flow, we have also used the correction (Eqs.34 and 35).In this case the subgrid diffusion coefficient was reduced twice: , c = 0.5 (coefficient c = 0.5 was chosen because about half of the subgrid flux can be restored using "defiltering": sw * − s w ≈ 0.5 ϑ s 3 ).We note that when the particles are advected by velocity u * (p) , then the presence or absence (crosses in Fig. 4c, d) of correction has no significant effect on the function f y s .Nevertheless, this procedure may be useful for the following reasons.
In the inertial range of three-dimensional turbulence along with the kinetic energy the variance of a passive scalar concentration is transferred from large scales to small scales with the formation of the spatial spectrum S s ∼ s −1/3 k −5/3 (see Obukhov, 1949) (here s is the dissipation rate of the variance of concentration, caused by molecular diffusion).Lagrangian transport of particles by a divergence-free velocity field u * with the truncated small-scale spectrum is equivalent to Eulerian advection of concentration s without any dissipation.The absence of a subgrid-scale part of the velocity spectrum will lead to reduction of the forward cascade and to the accumulation of variance σ 2 sp in the vicinity of the smallest resolved scales.
Figure 5a shows the variances of "Eulerian" concentration σ 2 s (z) = s div is used (filled circles), the values of σ 2 sp and σ 2 s become closer to each other.Besides, the correction (Eqs.34 and 35) increases the correlation corr(s, s P ) = s s P xyt /(σ sp σ s ) of two fields calculated by means of "Eulerian" and "Lagrangian" approaches (see Fig. 5b).
One can expect that in more complicated cases (e.g., the turbulent flow around geometric objects and the formation of quasi-periodic eddies), the accumulation of small-scale noise in the concentration field may lead to the incorrect advection of concentration by the resolved eddies.This effect may also be important for inertial particles when the nonphysical variance of concentration can directly affect dynamics.In additional tests it was found that the correction given by Eqs. ( 34) and ( 35) prevents particle stagnation in zones with unresolved turbulence during the modeling of urban-like environments.Thus, this correction is desirable for a number of reasons as a practical replacement of subgrid stochastics, which requires large computer resources.

Particle advection and footprint determination in LES with subgrid LSM
One can obtain footprints close to those presented in Fig. 4 by means of application of the stochastic subgrid model .The calculations for this model have been carried out on the grids with steps 3.125, 6.25 and 12.5 m (solid lines in Fig. 6a, b).One can note the shortcoming of this stochastic subgrid modeling in LES, which can not be detected by study of the mean characteristics.In the previous subsection, the recovered "subfilter" part of velocity u = u * − u and so the subfilter Lagrangian velocity u (p) were highly correlated with the resolved velocity u in time and space.This is due to the specifics of the spatial filter (Eq.24) used for the recovering given by Eqs. ( 25) and ( 26).This filter has a smooth transfer function in spectral space.The analogous effects of non-ideal filters in LES that lead to the high correlations between modeled and measured turbulent stresses were obtained and discussed earlier in Liu et al. (1994) and Meneveau and Katz (2000), where the laboratory data of turbulent flows were studied.By contrast, additional mixing in the stochastic model (Eqs.28-32) is due to random fluctuations, which are not related to u strictly.When one uses coarse grids, the energy of these Lagrangian fluctuations should be large enough to restore mixing in the vertical direction.This is accompanied by an excessive suppression of the variability of concentration s P near the surface, where the contribution of subgrid mixing is large (stars in Fig. 5a).The correlation between "Eulerian" and "Lagrangian" concentrations is reduced simultaneously (see Fig. 5b).Probably, this defect of the employed Lagrangian stochastic model is connected to the horizontal averaging in the evaluation of "subgrid" dissipation and energy.Nevertheless, this result shows that in some cases the stochastic subgrid modeling can prevent correct reproduction of the resolved spatial variability of particle concentrations in LES along with improvement of the mean transport.

Footprints in LES with subgrid RDMs and the comparison of different methods
In Fig. 7 footprints obtained in LES with intermediate resolution g = 6.25 m are shown.We choose this resolution because LES dynamics is still reproduced sufficiently well, but the effects from the subgrid/subfilter Lagrangian parameterization are already clearly visible.In addition to the approaches that were already discussed above, we applied the subgrid RDM (Eq.33) and the subgrid RDM in combination with the velocity recovering (Eqs.25 and 26) and the correction (Eqs.34 and 35).In the former case we restricted the activity of the subgrid RDM by the multiplying of the diffusivity coefficient K subgr(p) h in Eq. ( 33) on the following ramp function: r(z p ) = (1 − z p / g ) if z p ≤ g and r(z p ) = 0 if z p > g .Generally, results are in close agreement with the results of LES with the fine grid, except for some details.One can see the intrinsic defect of the RDM when it is applied to the dispersion of particles in a near field of a source.That is, as the RDM is the approximation of the diffusion process with the infinite speed of the signal prorogation, this model overestimates values of f y s in the vicinity of the measurement point location (see Fig. 7d, where this effect is highlighted in the logarithmic scale).Nearly the same effect was obtained in Wilson (2015) (see Figs. 1-3 in that paper, where the footprints from the RDM are also shifted left in comparison with the other models).It was also observed that, along with the overestimated vertical mixing, a subgrid RDM leads to the propagation of some portion of the particles in the upwind direction (the function f y s (x M − x, z M = 10) has small but positive values if x M − x < 0).In LES with the intermediate resolution the mentioned overestimated mixing exceeds the similar effect in RDM standing alone (see Sect. 5.3), be-cause the coefficient K subgr h is highly variable in time and space, and it can achieve even larger local values than the magnitude of the averaged turbulent diffusivity K h .At the higher levels of z M = 30 m and z M = 60 m, the footprints are formed as a result of averaging of the turbulent motions over the large spatial distances and over long temporal intervals, and the diffusion approximation becomes acceptable.As will be shown in Sect.5.3, an RDM applied alone gives very close results to the results of LSMs in this particular case of the stable ABL.
In contrast to the subgrid LSM and to the methods of velocity correction proposed above, the advantage of the subgrid RDM consists in the absence of the arbitrary prescribed parameters and in the absence of the need to involve the additional suppositions.In terms of Eulerian statistics, this model is identical to Eq. ( 6) (in the limit dt → 0 and with the precision defined by the spatial approximations).From this point of view, subgrid RDM can be considered the "ideal" model,  3,0x10 -4 3,5x10 -4 4,0x10 -4 x M -x (m) x M -x (m) 0 1000 2000 3000 4000 5000 6000 7000 8000 0,0 2,0x10 -5 4,0x10 -5 6,0x10 -5 8,0x10 -5 x M -x (m) Figure 7. Crosswind-integrated scalar flux footprints f y s , obtained in LES with g = 6.25 m using different stochastic Lagrangian subgrid models RDM (Eq.33) and LSM (Eqs.28-32).The results obtained with these subgrid models applied within the first computational grid layer in combination with velocity recovering u * = F −1 u and correction of velocity (Eqs.34 and 35) are also shown.Black lines are the footprints in LES with g = 2.0 m. because it is determined by the coefficients that are consistent with LES dynamics of the stratified flow (the same subgrid diffusivity is used for the potential temperature that defines the buoyancy and the interchanges between the kinetic and available potential energy).One can see that the variance of "Lagrangian" concentration computed with the use of a subgrid RDM (open circles in Fig. 5a) is very close to the variance of the concentration obtained by the Eulerian method.The correlation between "Eulerian" and "Lagrangian" concentrations (open circles in Fig. 5b) is also large, except for the first computational level; there, the Eulerian non-monotonous numerical advection scheme produces significant numerical noise.Thus, we have one more confirmation of the validity of the results, except for the invariance with respect to the grid steps.
The impact from the subgrid RDM is reduced when it is applied within the first grid layer only.In this case, the footprints are approximately the same as the footprints computed using the other approaches.

Two-dimensional footprints
The trajectories of a large number of particles (∼ 1.8 × 10 8 ) were simultaneously computed in LES with a grid step of 2.0 m.Accordingly, one can get a statistically grounded esti-mation of two-dimensional footprint functions f s (x−x M , y− y M , z M ).These functions, computed for the sensor heights z M = 10 m and z M = 30 m, are shown in Fig. 8a, b.One can see that the area with the negative values of the footprint exists.The negative values of the footprints are typical (e.g., Cai et al., 2010;Steinfeld et al., 2008) of the convective boundary layer due to fast upward advection by the narrow thermal plumes and slow downward advection in the surroundings.Here, the negative values of the function f s are connected to the Ekman spiral and to the mean transport of the particles elevated to large altitudes in the direction perpendicular to the near-surface wind.The negative values of the scalar flux footprint show that the vertical turbulent transport of the scalar emitted in the relevant area is basically directed from the upper levels down to the surface.For example, the positive surface concentration flux in this area will lead to a negative anomaly of the turbulent flux measured in the sensor position.This does not contradict the diffusion approximation of the turbulent mixing, because mean crosswind advection at the upper levels can produce the positive vertical concentration gradient to the right of near-surface wind.
The contribution of the negative part of the flux to the "measured" flux is significant, as shown in Fig. 8c, d, where x M -x (m) cumulative footprints, defined as are separated into positive and negative parts 5 Stochastic modeling and the comparison with LES

Preparation of turbulence data from LES for LSMs and RDMs
The LES results with grid step g = 2.0 m were used for data preparation.To apply an LSM (Eqs.8 and 9), the following Eulerian characteristics are required: the mean wind velocity components u and v , the second moments u i u j and the dissipation .Stochastic models are even more sensitive to some of these characteristics than the advection of particles in LES.For example, the underestimated values of the turbulent kinetic energy in LES are the consequence of the suppression of small eddies.Nevertheless, these eddies exert a relatively small influence on the mixing of scalars, because the effective eddy diffusivity associated with them (K small h ∼ E 1/2 small l small ) is not large due to the small spatial scale.However, the turbulent energy that is substituted into the LSM affects results independently of the scale and has to be evaluated with good accuracy.

Mean velocity
Mean wind velocity at the height z 0 < z ≤ g was computed using log-linear law: and u i = 0 at z < z 0 .Here, U * is the friction velocity, κ = 0.4 denotes the von Karman constant, L is the Obukhov length at the surface where Q s is the kinematic potential temperature flux at the surface, g = 9.81 m s −2 is the acceleration of gravity and A. Glazunov et al.: LES and LS modeling of Lagrangian particles for footprint determination in the SBL Note that the von Karman constant is not included in the definition of the length L here and later (this alternative definition of the Obukhov length is used along with the traditional one; see, e.g., Zilitinkevich et al., 2013, Eq. 41).The linear interpolation of velocity was used if z > g .

Momentum fluxes
The fluxes u i u j = u i u j + τ mix ij (i = j ) were interpolated linearly and additionally smoothed everywhere in the domain.These fluxes are shown in Fig. 9a.

Variances of velocity components
The variances of velocity components σ 2 i = u i 2 were estimated by the formula where E subg is the subgrid energy (Eq.32) and (u * i ) 2 are the variances of recovered velocity components.The vertical velocity variance has the greatest impact on the functions f y s .Figure 9b shows the comparison of evaluated normalized RMS σ w = σ w /|τ | 1/2 (solid line) with the SHEBA data (symbols; see description in (Grachev et al., 2013, Fig. 15b); data kindly provided by Dr. A. Grachev).The data are shown in dependence on non-dimensional stability parameter ξ = κz/ , where is the local Obukhov length, determined using values of fluxes of momentum |τ | and temperature Q at the given height z (local scaling in the stable ABL (Nieuwstadt, 1984)).
The measurements suggest that the mean value of the normalized RMS σ w ≈ 1.33 if the value ξ is small.Figure 9b shows that our estimation of RMS is slightly less than the measured values in the interval 0.03 < ξ < 0.2.Respectively, the final values of vertical velocity variance designed for the substitution in stochastic models were corrected as follows: At the higher levels, the estimation (Eq.41) was applied.
The final estimations of the variances of velocity components are shown in Fig. 9c by the solid lines.Dashed lines are the filtered resolved velocity u i variances.The estimation of the variance σ 2 w using Eq. ( 41) is shown by the circles.One can see that significant parts of variances were not reproduced explicitly in LES and were recovered using the above-mentioned assumptions.

Turbulent energy dissipation rate
Usual interpolation is not applicable to the calculation of dissipation rate near the surface, where ∼ 1/z.Besides, the values of dissipation k computed in LES at the levels z k = (k − 1/2) g are approximately equal to the averaged values inside the layers (k − 1) g < z ≤ k g , but not to the physical dissipation at given altitudes.Under the assumption that |τ | is constant with height and neglecting the stratification inside first layer, one can get the following corrected value of at the height z = g /2: Additional analysis showed that, if z < 0.25z i , then the local balance of turbulent kinetic energy (TKE) is well satisfied: ≈ S + B, where S and B are shear and buoyancy production.Therefore, the non-dimensional dissipation can be approximated by a formula where is the non-dimensional velocity gradient; C m = 5, according to the observation data (e.g., Grachev et al., 2013) and LES results (e.g., Glazunov, 2014b).Here, the assumption is used that the shear ∂ u /∂z and the stress τ are collinear.Previous LES studies of the stable ABL (e.g., Beare et al., 2006) also give negligibly small values of the transport terms in the TKE balance.The experimental confirmation of the validity of Eq. ( 44) can be found in Grachev et al. (2015), where the dissipation in the stable ABL was estimated using the spectral analysis of longitudinal velocity in the inertial range.In accordance with this paper, ≈ φ m , which is almost indistinguishable from Eq. ( 44) within the accuracy of the experimental data and the ambiguity of the method of dissipation evaluation.Discrete values of non-dimensional dissipation k z k /|τ | 3/2 are shown in Fig. 10a by circles.The dashed straight line is the universal function (Eq.44).One can see that the correction (Eq.43) makes the dissipation values closer to the function (Eq.44).Finally, the profile of dissipation cf (z) for the LSM was corrected as follows (see Fig. 10b).The dissipation was set to be constant below some height z e , and was replaced by the universal function = |τ | 3/2 /z up to the level with z/ = 1.The height z e was chosen in such a way to equalize values of the dissipation averaged in a layer 0 ≤ z ≤ g and the dissipation 1 .Figure 10b shows that the corrected dissipation cf (solid line) is very close to "discrete" dissipation k (circles), except for the first computational level.

Diffusion coefficients
A random displacement model (Eq.15) requires the estimation of an eddy-diffusion coefficient K s .Note that, due to + 2/3E subgr ) 1/2 ; symbols -measurements (Grachev et al., 2013)  anisotropy, one should use tensor diffusivity K ij s in a general case.Neglecting this fact, let us assume that the principal axes of the tensor K ij s are aligned with the coordinate axes.The corresponding coefficients K ww s , K uu s and K vv s (see Fig. 9d) can be calculated as follows: The horizontal eddy diffusivities K uu s and K vv s are estimated taking into account Eq. ( 13).
One can see that the formula (Eq.13) provides a good approximation for the coefficient K ww s if one sets the value C 0 = 6.We note that the data of LES were substantially corrected to get this estimation.Very fine grid simulations are needed to verify and justify the given value.There is no guarantee that this constant is actually universal under different stratifications in the ABL.

Specification of LSMs and RDMs tested against LES
The following stochastic models were tested using the data prepared as described above.
1. RDM0 is the random displacement model with uncorrelated components.Particle position is computed by the formula similar to Eq. ( 15) but with direction-dependent coefficients (see Eqs. 46 and 47 and Fig. 9d).The components of the Gaussian random noise satisfy Eq. ( 10).
2. RDM1 differs from RDM0 by using the noise with inter-component correlations: where www.geosci-model-dev.net/9/2925/2016/ 3. LSM0 is the Lagrangian stochastic model without a WMC: 4. LSM1 is based on the one-dimensional well-mixed model: supplemented by uncorrelated horizontal mixing similar to Eq. ( 49) with the appropriate variances σ 2 u and σ 2 v .
5. LSMT is a three-dimensional Lagrangian stochastic model satisfying a WMC, which is proposed by Thomson (1987).For the incompressible turbulent fluid in a steady state and under the condition of zero mean vertical velocity, this model (Thomson, 1987, formula 32) reads where τ −1 is the tensor inverse to the stress tensor.
The setups of numerical experiments with RDMs and LSMs were close to particle advection conditions in LES (absorbtion at altitude 100 m, ejection at z 0 = 0.1 m and reflection at z = 0).The particles were generated continuously within 2 h of modeling.The last hour was used for averaging.Models LSM0 and LSM1 use the value C 0 = 6.Threedimensional model LSMT was applied with C 0 = 6, C 0 = 8 and C 0 = 4.

Modeling results
Figure 11 shows crosswind-integrated footprints f y s and the corresponding cumulative footprints F , computed by LES (black bold solid lines, g = 2.0 m) and by stochastic models described above.Footprints are shown for sensor heights z M = 10, 30 and 60 m.Models RDM0, RDM1 and LSM1 provide very similar results.Faster mixing is observed in stochastic models below altitude z M = 10 m in comparison to LES.These differences are not crucial and are compensated for in cumulative footprints at the distances x − x M ∼ 1000 m.The differences can be explained either by insufficient subgrid mixing in LES or by an inexact procedure of the data preparation for stochastic modeling.Very weak sensitivity of the models with respect to correlations of particle velocity components is observed as well.Thus, the results close to LES were obtained in stochastic models having the "diffusion limit" with the same or similar vertical diffusion coefficient.The significant advantages of LSMs compared to RDMs were not observed in this particular flow.
The substantial disagreements with LES were obtained using the three-dimensional Thomson model (Eq.51), with C 0 = 4 and C 0 = 6, and the LSM0 model.The last one is designed for the isotropic turbulence and does not satisfy a WMC under the conditions considered here.This model leads to overestimated mixing, and such bias does not vanish at high altitudes.
LSMT (Eq.51) was proposed in Thomson (1987) as one of the possible ways to satisfy a WMC in three dimensions.In our simulations the error of LSMT with C 0 = 6 is substantial and grows with sensor height.This was shown by  (Thomson, 1987) with C 0 = 6 (absorption at z = 100 m); open triangles -LSMT with C 0 = 8; blue dashed lines -LSMT with C 0 = 4. Orange curves -LSMT with C 0 = 6 (absorption at z = 300 m).Short-dashed line -LSM0 (Lagranian stochastic model without a well-mixed condition).Red circles -LSM1 (an LSM with a WMC for vertical mixing).Open green circles -RDM0 (uncorrelated random displacement model).Dash-dot green line -RDM1 (random displacement model with correlation between the displacement components).Sawford and Guest (1988), who derived the diffusion limit of Thomson's multi-dimensional model for Gaussian inhomogeneous turbulence and showed that the implied effective eddy diffusivity for vertical dispersion is Taking into account this expression and Eq. ( 13), which is valid for the one-dimensional LSM, one can estimate the appropriate value of C 0 for LSMT under the conditions considered here: C 0 ≈ 6(1.33 4 + 1)/1.33 4 ≈ 8 (we assume that The results of LSMT with C 0 = 8 are in close agreement with the results of other stochastic models and with the results of LES (open triangles in Fig. 11a, c, e).
One can see that Thomson's multi-dimensional model with C 0 = 4 produces a very short footprint (blue dashed lines in Fig. 11).Similar results can be obtained using LSM1 with the values of C 0 < 6 (not shown here).
Finally, it can be seen from Fig. 11 that the top boundary condition (absorbtion of particles at the height of 100 m) does not affect the footprints considered here.See the orange curves, which are obtained in LSMT (C 0 = 6) with an absorption condition applied at the level above the boundary layer height.
Turbulent Prandtl P r and Schmidt Sc numbers computed using the Eulerian approach are shown in Fig. 12a.These numbers coincide and are approximately equal to 0.8 up to the altitude slightly less than 100 m, where the boundary condition for a scalar is applied.Schmidt numbers Sc were also calculated using the concentrations and the fluxes of Lagrangian particles.Models RDM0 and LSM1 provide the values of Sc close to the results of the Eulerian model.Calculations by LSMT (C 0 = 6) result in Sc ≈ 0.5-0.6,which is also the sign of the overestimated vertical mixing.Two-dimensional footprints f s (x −x M , y −y M , z M ), computed by models RDM0, RDM1 and LSM1 (figures are not shown here), were very close to LES results presented in Fig. 8.In particular, this fact argues for the mechanism of formation of the region with negative values of f s having a simple nature, which can be easily reproduced in the framework of the diffusion approximation.
The crosswind mixing can be characterized by an RMS of transversal coordinates of the particles depending on the mean distance from the source: Y p (X p ) = (y p − Y p ) 2 1/2 , where X p = x p and Y p = y p are the mathematical expectations of the particle position.Functions Y p (X p ) are shown in Fig. 12b.Models RDM0, RDM1, LSM1 and LSMT (with C 0 = 6) result in close horizontal dispersion.All the stochastic models predict slightly less intensive mixing in comparison to LES, which can be a consequence of the inaccurate data preparation algorithm as well.If one neglects the anisotropy of eddy diffusivity, then this dispersion would be substantially underestimated (see the shortdashed line in Fig. 12b, computed by an RDM with coefficients K uu s = K vv s = K ww s ).One can see that choice C 0 = 8 in LSMT (open triangles) does not improve its overall performance because the improved vertical mixing is accompanied by the reduced dispersion of particles in the horizontal direction.
Wind direction rotation leads to widening of a concentration trace from the point source (see the thin dashed line in Fig. 12b, computed with a one-dimensional LSM).At larger distances from the source in the Ekman layer the crosswind dispersion of pollution should be defined by the joint effect of the wind rotation and vertical mixing, but not by the horizontal turbulent mixing.

Conclusions
Scalar dispersion and flux footprint functions within the stable atmospheric boundary layer were studied by means of LES and stochastic particle dispersion modeling.It follows from LES results that the main impact on the particle dispersion can be attributed to the advection of particles by resolved and partially resolved "subfilter-scale" eddies.It ensures the possibility of improving the results of particle advection in discrete LES by the use of recovering of small-scale partially resolved velocity fluctuations.If one uses the LES model with the explicit filtering, then this recovering is straightforward and consists of application of the known inverse filter operator.Apparently, a similar method can be implemented for other LES when the spatial filter is not specified in an explicit form.This would require, however, the prior analysis of the modeled spectra to identify an effective spatial resolution and the actual shape of the implicit filter.For substantial improvement of particle transport statistics, it is enough to use a subgrid Lagrangian stochastic model within the first computational layer only, where the LES model becomes equivalent to the simplified RANS model.
When the particles are advected by a divergence-free turbulent velocity field, then the variance of the particle concentration can be accumulated at small spatial scales.In the considered case, it does not affect directly the particle advection by the large eddies and has no significant influence on the results of footprint calculations.In those cases, when the instantaneous characteristics of the scalar field of a particle concentration are important, additional correction to particle velocities may be required.It can be done both through the introduction of stochastics, resulting in the diffusion of concentration, and through the "computationally inexpensive" divergent correction of the Eulerian velocity field.
Under the stable stratification, to calculate the flux footprint, it is preferable to use stochastic models, which describe the particle dispersion close to the process of scalar concentration diffusion with the effective coefficient K ww s (z) = − w s /( ds /dz) in a vertical direction.RDM and onedimensional "well-mixed" LSM tested in this study are the examples of such stochastic models.The optimal value for the parameter C 0 for LSMs is found to be close to 6 under the conditions considered here.This value coincides with the estimation of Kolmogorov Lagrangian constant in isotropic homogeneous turbulence.It provides additional justification for use of LSMs in stable ABL, due extending their applicability over a wider range of scales including the inertial subrange.Stochastic models that use smaller values C 0 ≈ 3-4 (this choice is widespread now) may produce extra mixing and the shorter footprints, respectively.Note that the estimation C 0 = 6 is based on the LES results combined with the SHEBA data (Grachev et al., 2013), where the non-dimensional vertical velocity RMS was evaluated as σ w ≈ 1.33 (the exact estimation of this value in LES is restricted by the resolution requirements).In the cases when LSMs utilize smaller values of σ w the parameter C 0 should be reduced accordingly (for example, C 0 ≈ 4.7 will be the best suited parameter for LSMs with the widely used value σ w ≈ 1.25 prescribed).
One-dimensional stochastic models can be supplemented by the horizontal particle dispersion in a simple way.Introduction of the correlation between particle displacement components in RDM does not improve or change results substantially.However, the coefficients of horizontal diffusion K uu s and K vv s for RDMs can be evaluated through the vertical diffusion coefficient K ww s multiplied by the square of velocity component variances ratio.
Model LSM1, constructed as a combination of independent stochastic models in each direction (well mixed in the vertical direction only), gives reasonable results, although this model does not satisfy a WMC in general.In contrast, the three-dimensional Thomson model with a WMC and C 0 = 6 provides overestimated vertical mixing, which is manifested in too small Schmidt number values and in reduced lengths of the footprints.The Thomson model with C 0 = 8 produces true mixing in the vertical direction, but underestimates the mixing in the crosswind direction.
Accordingly, one can recommend another well-mixed stochastic model proposed in Kurbanmuradov and Sabelfeld (2000).It was developed under the assumption that the vertical drift term does not depend on the horizontal velocity components, and the vertical component of this model coincides with LSM1.Prior to use, this model should be modified in an appropriate way to take into account the variation of momentum fluxes with height.
According to the presented LES, the source area and footprints in the stable ABL can be substantially more extended than those predicted by the modern LSMs and footprint parameterizations based on their results (e.g., the parameterization by Kljun et al., 2015, which was calibrated with the use of a stochastic model Kljun et al., 2002).The following reasons were identified in this study: (1) too small values of the parameter C 0 are used; and (2) the possible overestimated vertical mixing provided by some stochastic models based on a well-mixed condition.
We emphasize that a very simple case of the moderately stratified stable ABL in almost steady-state conditions was considered here.This setup of numerical experiments permits the detailed intercomparison of different approaches for the particle dispersion modeling, which utilize identical simplifications.On the other hand, in a real environment the scalar flux footprint functions can be greatly influenced by the meteorological non-stationarity, the peculiarities of mixing inside the roughness layer, internal radiative heating or cooling in the ABL, and so on.Also, a wider investigation of different stability regimes from neutrality to strong stratification must be undertaken in future studies to confirm the universality of the findings.

Code and data availability
The code of the LES model is available on request for scientific researches in cooperation with the first author (and.glas@gmail.com).The data from LES are attached to the Supplement.These data were prepared as was discussed in Sect.5.1 and can be used for the stochastic models' evaluation.Besides, the Supplement contains the data for crosswind-integrated footprints and two-dimensional footprints obtained in LES (see Figs. 6 and 9).
The Supplement related to this article is available online at doi:10.5194/gmd-9-2925-2016-supplement.

Figure 1 .
Figure 1.Schematic representation of the footprint evaluation algorithm.(a) Setup of the numerical experiment.(b) Example of two trajectories (red and blue bold curves).Shifted trajectories are shown by the dashed lines.The particle brings the impact into the value f s (x S , y S , x M ) if it intersects the test area δ M in the vicinity of the sensor position x M and the origin of the modified trajectory belongs to the test area δ S .
x M , y p 1 − y M ) ∈ δ M in time interval T a ; I p SM = 1 if the www.geosci-model-dev.net/9/2925/2016/Geosci.Model Dev., 9, 2925-2949, 2016 A. Glazunov et al.: LES and LS modeling of Lagrangian particles for footprint determination in the SBL (y M −y S )) ∈ δ S and I p SM = 0 otherwise.Here, w p is the vertical component of the particle velocity at the moment of crossing the plane z = z M .Schematic representation of the algorithm for the footprint function determination in LES is shown in Fig. 1.In accordance with Eq. ( y, 0) = 0.A constant value K subgr s (x, y, z) = K subgr s (x, y, z 0 ) was used for z < z 0 .

4
LES of stable ABL and footprint calculations 4.1 The setup of the numerical experimentStable boundary layer at the latitude 73 • N in almost steady state conditions was considered.The calculations were carried out according to the GABLS scenario(Beare et al., 2006), with the difference that the geostrophic wind U g has been rotated 35 • clockwise such that the wind direction near the surface approximately coincides with the axis x.The duration of runs is 9 h.The initial wind velocity coincides with geostrophic velocity |U g | = 8 m s −1 .The initial potential temperature is equal to the surface temperature s | t=0 = 265 K up to the height 100 m and increases linearly with the rate d /dz = 0.05 K m −1 if z > 100 m.During the calculations, the surface temperature decreases linearly with time: d s /dt = −0.25 K h −1 .Dynamical and thermal roughness parameters z 0 and z 0 are set to 0.1 m.

Figure 2 .
Figure 2. Mean wind velocity u (a) and temperature (b) in runs with different grid steps (spatial step is pointed in legend).Gray dots are the data from other LES models obtained in Beare et al. (2006) (wind velocity is rotated 35 • clockwise)."Standard" wind profile for stable conditions in accordance withHögström (1996) is shown by the vertical dashes.

Figure 3 .Figure 4 .
Figure 3.Total F tot s = s w + ϑ s 3 (solid lines), resolved F res s = s w (short-dashed lines) and "subgrid" F sbg s = ϑ s 3 (long-dashed lines with shading) scalar fluxes in the runs with different grid steps g .
if particles are advected by the velocity u * (p) (crosses), variance σ 2 sp is much larger than σ 2 s .If the velocity u * (p) + u (p)

Figure 5 .Figure 6 .
Figure 5. (a) Variance σ 2 s = s 2 of the concentration of Eulerian scalar (solid lines) and variance σ 2 sp = s P 2 of concentration s P , determined by Lagrangian particles (symbols); grid steps and the methods of calculations are shown in the legend, and symbolic notations are the same as in Fig. 4; stars -a stochastic model (LSM, Eqs.28-32) is used throughout the domain.Open circles -a subgrid RDM (Eq.33) is applied.(b) Correlation corr(s, s P ) = s s P xyt /(σ sp σ s ) between "Eulerian" and "Lagrangian" concentrations.For remaining notations, see the caption of Fig. 4.
Figure 8. Two-dimensional footprints f s (x − x M , y − y M , z M ) (×10 −6 m −2 ) for sensor height z M = 10 m (a) and z M = 30 m (b) and the corresponding crosswind-integrated cumulative footprints F (x M − x) (c) and (d); long dashed line -F + (impact of the area with positive values of f s ); short dashed line -F − (impact of area with negative values).

Figure 9 .
Figure 9. (a) Total momentum fluxes obtained in LES with g = 2.0 m.(b) Normalized RMS of vertical velocity σ w = σ w /|τ | 1/2 depending on a dimensionless parameter z/ (solid red line -estimation using LES data σ w = ( w * 2+ 2/3E subgr ) 1/2 ; symbols -measurements(Grachev et al., 2013) at different heights).(c) Variances of velocity components (dashed line -resolved fluctuation; solid lines -the final estimation for LSM; bold red lines -vertical component, green curves of medium thickness -crosswind component, blue thin lineslongitudinal component, circles -evaluation of σ 2 w by Eq. 41).(d) Vertical effective eddy diffusivity K ww at different heights).(c) Variances of velocity components (dashed line -resolved fluctuation; solid lines -the final estimation for LSM; bold red lines -vertical component, green curves of medium thickness -crosswind component, blue thin lineslongitudinal component, circles -evaluation of σ 2 w by Eq. 41).(d) Vertical effective eddy diffusivity K ww s (red solid line -coefficient calculated by the gradient and flux of scalar; dashed line -estimation of coefficient using Eq.(13) with C 0 = 6); estimations of diffusion coefficients in crosswind direction K vv s (green dash-dot line) and coefficient in longitudinal direction K uu s (blue dash-dot-dot line).

-
Figure 10.(a) Discrete (LES) non-dimensional dissipation k κz k /|τ | 3/2 (circles), corrected values (solid line), universal function (Eq.44) (dashed straight line).(b) Simulated discrete dissipation k (circles) and corrected dissipation cf (z) for LSM (solid line).Dashed horizontal line denotes the height z e , which was chosen in order to equalize the integral values of the corrected dissipation and the discrete dissipation.

Figure 12
Figure 12.(a) Prandtl number P r (dashed line) and Schmidt number Sc (solid line), computed using Eulerian scalars.Symbols -Schmidt numbers Sc, computed using the Lagrangian particles in LES, LSMs and RDMs.(b) RMS of the crosswind position of particle Y p = (y p − Y p ) 2 1/2 depending on the mean longitudinal position X p = x p .Dashed lines -RDM with K uu s = K yy s = K ww s and