Comparison of conventional Lagrangian stochastic footprint models against LES driven footprint estimates

In this study we introduce a comparison method for footprint model results by evaluating the performance of conventional Lagrangian stochastic (LS) footprint models that use parameterised flow field characteristics with results of a Lagrangian trajectory model embedded in a large eddy simulation (LES) framework. The two conventional models follow the particles backward and forward in time while the trajectories in LES only evolve forward in time. We assess their performance in two unstably stratified boundary layers at observation levels covering the whole depth of the atmospheric boundary layer. We present a concept for footprint model comparison that can be applied for 2-D footprints and demonstrate that comparison of only cross wind integrated footprints is not sufficient for purposes facilitating two dimensional footprint information. Because the flow field description among the three models is most realistic in LES we use those results as the reference in the comparison. We found that the agreement of the two conventional models against the LES is generally better for intermediate measurement heights and for the more unstable case, whereas the two conventional flux footprint models agree best under less unstable conditions. The model comparison in 2-D was found quite sensitive to the grid resolution.


Introduction
Footprint modelling aims at determining the areas of highest influence on concentrations or fluxes of atmospheric con-Correspondence to: T. Markkanen (tiina.markkanen@fmi.fi)stituents at a certain location.This is necessary in interpreting the results of measurements especially when those are performed over a landscape of varying source strengths.While the term footprint refers to the field of view of a measurement device, the outcome of footprint modelling typically is the footprint function or source weight function, which provides information about the relative weights of individual point sources (see Schmid, 2002, and the references therein).Analytical footprint models facilitate an analytic solution for the diffusion equation and are usually applied only for flow within well developed, stationary atmospheric surface layers (ASLs) where Monin-Obukhov scaling is valid.In the Lagrangian stochastic (LS) approach a large number of particles are followed as they traverse between their sources and the observation point.Stochastic models are more time consuming but they can also be applied for measurements over tall canopies where the within canopy flow and alongwind diffusion are of crucial importance and roughness layer effects modifies the characteristics of turbulence (e.g.Baldocchi, 1997;Rannik et al., 2000Rannik et al., , 2003)).LS simulations can be run either backward (e.g., Kljun et al., 2002) or forward (e.g.Leclerc and Thurtell, 1990;Rannik et al., 2000) in time.The heterogeneities in the flow field and even nonstationarity are more straightforward to take into account in backward modelling, whereas forward models are less sensitive for stochastic noise than backward models in which the flux contributions are dependent on the ratio of initial and touchdown velocities (Flesch, 1996).The conventional way to simulate LS dispersion requires pre-determined flow characteristics as input.Large eddy simulations (LES) instead determine the statistics of turbulence themselves according to the given initial and boundary conditions.Those flow field characteristics are consequently either used as input for T. Markkanen et al.: Lagrangian footprint model comparison in 2-D Lagrangian particles (Cai and Leclerc, 2007) or the particles are released during LES run in which case the particles experience the flow field development in real time (Steinfeld et al., 2008).
In the present comparison of footprint models we use the LS model LPDM-B (Kljun et al., 2002) for backward simulations (BW) and the LS model by Univ. of Helsinki (Rannik et al., 2000;Markkanen et al., 2003;Rannik et al., 2003) for forward simulations (FW).As a reference we use the LES model PALM (Raasch and Etling, 1998;Raasch and Schröter, 2001) which simulates trajectories of a large number of particles simultaneously with general flow field calculations (Steinfeld et al., 2008).From this data the footprints are determined in a manner similar to that used in conventional forward LS models (Kurbanmuradov et al., 1999).
In previous footprint model comparison studies the most common measures used are peak position and certain statistical measures of cumulative footprint function.Schmid (1997) defined for 2-D source areas the effect level of a given percentage of the total source weight function, Markkanen et al. (2003) used the distances of cumulative 1-D footprint function reaching a given percentage; the former cumulates the values from the maximum toward areas of smaller contribution up to the selected percentage, and the latter cumulates the crosswind integrated footprint function from infinity downwind from the measurement point towards increasing distances upwind.Both in model comparison studies (Horst and Weil, 1992;Rannik et al., 2000;Kljun et al., 2003) and in studies aiming at model evaluation based on tracer experiments (Finn et al., 1996;Leclerc et al., 2003b,a;Mölder et al., 2004) as well as in wind tunnel experiments (Kljun et al., 2004) crosswind dimension is often neglected, and the results are shown in 1-D -that is, as functions of along mean wind distance from the observation position.The reason for neglecting crosswind dimension in the analysis of LS simulations is often the excessively high stochastic noise of 2-D patterns.On the other hand, any comparison is more straightforward as a function of one Cartesian dimension than of two dimensions.Quite often the crosswind direction is not considered in the simulations at all.In the tracer experiment evaluation case, the reason for exclusion of crosswind direction may also be that the experimental set up only provides information for the vertical dimension.For practical purposes, however, the source areas have to be considered in 2-D.For instance, quality assessment tools (e.g., Göckede et al., 2007Göckede et al., , 2008) make use of 2-D footprints to determine the areas of importance surrounding the measurement point.
Nowadays, as computing power has increased, the run times for simulations with high enough particle numbers have decreased and even the 2-D data can be used as such without smoothing the data.This is especially important when the 2-D footprint function is not symmetrical over the mean wind axis, as in the case when Ekman layer wind direction turning is taken into account (Steinfeld et al., 2008).
Moreover, in a comparison the tendency of the cumulated footprint function of each model has to be considered to make the simulated results truly comparable.According to Horst and Weil (1994), in the case of the uniform surface flux the integrated flux footprint function is required to tend to unity under homogeneous flow within the surface layer.Higher within the atmospheric boundary layer (ABL) (in fact in the convective boundary layer) it tends to 1−z m z −1 i (Horst and Weil, 1992), where z m stands for measurement height and z i is ABL depth.While many members of the family of analytical flux footprint models do not satisfy the above condition, Haenel and Grünhage (1999) developed a 1-D model that tends to unity once integrated to infinity.
As the upper boundary of particle dispersal is often neglected when footprints are determined for low measurement heights, the footprints derived from LS simulations are usually normalised by their cumulative values at certain horizontal distances from the observation point.Nevertheless, when the whole depth of the ABL is considered, the plume reflection from the boundary layer top influences the shape of the flux footprint function strongly, which often does not cumulatively tend to any well defined value within reasonable distances from the observation point.This is due to locally negative flux footprint functions which are to be expected under complicated flow situations as was demonstrated by Finnigan (2004) in a case of convergent flow near a hill top, whereas under simple shear flow situations the integrated flux footprint function is bound by 0 and 1.
Among LS approaches, the cumulative crosswind integrated flux footprint function presented by Leclerc and Thurtell (1990) tends to unity.The versions of the forward flux footprint model by Univ. of Helsinki (Rannik et al., 2000;Markkanen et al., 2003;Rannik et al., 2003) which only consider the surface layer flux footprint functions tend to unity.This is due to the fact that the flow field fulfils the simple shear flow condition discussed by Finnigan (2004).In Kljun et al. (2004) the parameterisation of LPDM-B model satisfies the integral condition of unity, as well.
According to Finnigan (2004) cumulated concentration footprint functions are expected to be bound by zero and one as they can be interpreted as Green's functions of Eulerian mass conservation equation or transition probabilities in a Lagrangian framework.In work by Cai and Leclerc (2007), that combines LES derived flow characteristics with LS simulations, concentration footprints tend to unity at long distances, showing, however at certain distances values well over unity.On the contrary, in the footprint model by Univ. of Helsinki (Rannik et al., 2000) which derives the footprints from dispersion data according to Kurbanmuradov et al. (1999) the concentration footprint functions tend to infinity when normalised only by the simulated particle number (which corresponds to normalisation by source strength).
Because there is a large demand for utilising outcome of footprint analysis in 2-D, but for various reasons listed above, an obvious deficit in methods for their validation, we present in this work a concept to compare the 2-D Lagrangian footprints.The results are presented for observation heights extending through the whole depth of the ABL.For that aim we present a way to normalise the stochastic particle data that is suitable for the data sets in hand.We quantify model agreement both for the equality of their sizes and for the degree of overlapping source area.We asses the impact of grid resolution on the comparison of gridded footprint data sets and finally we present a classification for 2-D model agreement.

The LS model embedded into the LES model PALM
The LS particles are embedded into the LES model PALM whose range of applicability covers boundary layers from convective to weakly stably stratified (Beare et al., 2006;Steinfeld et al., 2008) including neutral stratification (Letzel et al., 2006).The method for particle inclusion is based on Weil et al. (2004) who separate the velocity scales of particles into two scales consisting of deterministic and stochastic parts, corresponding to the division of the turbulent flow field into grid scale and sub-grid scales, respectively.For stochastic transport Weil et al. (2004) adopted the Thomson (1987) model which assumes isotropy and Gaussianity of turbulence (see Weil et al., 2004, for more details).The grid scale flow characteristics are interpolated -linearly in the vertical and bilinearly in the horizontal -to sub-grid scale particle positions.Following Kim et al. (2005), no boundary condition has been used at the ABL top, whereas at the top of the model domain a reflection condition has been applied.Nevertheless the latter condition had no effect as no particle reached the top of the domain.
Importantly, coupled LS simulation of the particles is used in PALM, whereas LES driven LS simulations by Weil et al. (2004), Cai and Leclerc (2007) and Kim et al. (2005) used pre-calculated LES data for subsequent, separate LS simulations.The latter approach is costly concerning the disc space, and limited by the writing and reading rates of the data.Furthermore, the PALM embedded LS calculations are fully parallelised which facilitates the release of exceptionally high number of particles.For more discussion on benefits of our approach see Steinfeld et al. (2008).
In this work PALM was driven in its dry mode and cyclic lateral boundaries were applied both for the flow field and the particles.Furthermore, Monin-Obukhov similarity was applied between the surface and the first computational grid point level.The horizontal distance between two sources of particles was 8 m and alltogether 8 particles were released at each source.The horizontal extension of the computational domain was 2560 m and the vertical extension was 960 m.The simulated time in the LES was 5 h.The particles were released after 10 800 s and particle trajectories were evaluated over the following 7200 s.As Weil's method of coupled footprint calculation has basically problems in cases in that the ABL height changes considerably with time, a strongly stable stratification above the neutrally stratified layer was prescribed at the beginning of the simulation.In fact the ABL height did not change much with time.In the case 1 (cf.Table 1), after 10 800 s when the particles were released the boundary layer height was 530 m, 2 h later the ABL height had not increased above 550 m.

Backward model LPDM-B
In LPDM-B the dispersion is based on a model by Rotach et al. (1996) and de Haan and Rotach (1998) which satisfies the well mixed condition by Thomson (1987) from convective to stable stratifications and over the whole depth of the atmospheric boundary layer.We used the model in its most parameterized form in which only surface roughness length, friction velocity, Obukhov length, convective velocity scale and boundary layer height are required as input (see Table 1 for the values of the parameters).Parameterizations of mean wind speed and standard deviations are given in Rotach et al. (1996).Calculation of backward trajectories and the method of deriving the flux and concentration footprint function out of the release and touchdown velocity data are given in Flesch et al. (1995) and Flesch (1996).The particles are reflected near the surface and at the top of the ABL as described in Rotach et al. (1996).For detailed description and sensitivity analysis of the model as whole the reader is referred to Kljun et al. (2002).In this study the number of simulated particles varied between 262 000 and 441 000.

Forward model
The forward model by Univ. of Helsinki (Rannik et al., 2000) was used in this work as a version presented in Rannik et al. (2003).The model simulates transition of stochastic particles according to Thomson's (1987) 3-D model forward in time.The model only considers dispersion within the ASL where Monin-Obukhov similarity is assumed.Universal functions accounting for stabilities for wind statistics, for wind speed and for dissipation rate of turbulence are given in Rannik et al. (2003).As the model version Rannik et al. (2003) actually considers the roughness sublayer, including the canopy sublayer, the leaf area density was adjusted so that roughness length of 0.14 m was attained (see Table 1).Because the canopy height (0.75 m) was low compared to the lowest measurement height (10 m) and the particles were released from the top of the canopy, the inclusion of the canopy layer is not expected to have a strong influence on the particle dispersion.Furthermore, release from the canopy top is considered to be representative of the surface sources of the ABL flow.In the simulations the particles that nevertheless reach the actual surface of the canopy layer, are perfectly reflected.Because Table 1.Parameters characterising the site, flow situation and the domain used in both runs.Parameter set for case 1 is adopted from Leclerc et al. (1997) flow only in the ASL instead of the whole ABL is parameterised, the reflection at the ABL top is not considered.Altogether 300 000 particles were released in the forward model simulations.

Footprint calculations
In the following presentation the form of the functions is adopted from Vesala et al. (2008) except for the notation convention, where capital letters in speeds and positions refer to particles, while lower case letters stand for fixed Eulerian reference frames such as positions of the measurement point.
According to Flesch (1996) in the case of the BW model the vertical flux density at the location (x, y, z) is where summations run over the total number of particles (N) and touchdowns (n i ), W i0 is initial speed of the particle and W ij is its speed at touchdown.Q stands for source strength.2-D footprint functions are subsequently calculated as follows: The equation for concentration at the sensor location differs from the respective equation for flux in the sense that initial speed of particles is not considered in the nominator.Thus the equation is as follows: The footprint function, however, is of a form similar to that of flux footprint.
In the FW model case, the respective estimator for flux at position (x, y, z) is as follows: where n i stands for number of interceptions of the particle with the measurement level z m and, consequently, W ij represents the particle velocity at the moment of interception.
Otherwise the notation follows that of backward calculations (Eqs. 1 to 3).Concentration at (x, y, z) is as follows The respective footprints are determined as in the backward model case above (Eq.2).
In this study source strength Q is horizontally homogeneous within the domain thus normalisation by source strength in Eq. 2 simply reduces the variable out of consideration.

Model comparison
In this work we simulated the convective, homogeneous boundary-layer already investigated in the paper of Leclerc et al. (1997) (case 1 with Obukhov length of −32 m, cf.Table 1).Additionally, we simulated a less unstable case 2 with an Obukhov length of −77 m.Since the two conventional LS models in their present form only consider stationary flow conditions, the state of a well developed, very slowly growing ABL situation was considered for the LES as well.The parameters characterising the site and the flow are given in Table 1.
We simulated trajectories for ten observation heights (10, 30, 40, 50, 100, 150, 200, 250, 300 and 350 m) for a mixed layer of approximately 500 m (see Table 1).Nevertheless, FW simulations were only performed up to an observation height of 100 m, which is already beyond the validity of the model as it only considers ASL flow.For each measurement height, we calculated the footprint contributions for evenly spaced grids of resolutions of: 5, 10, 20, 50, 100, 200 and 400 m.A wide range of resolutions was applied to be able to select the most suitable resolutions for each measurement height (see the end of this section for the selection of resolution).Selecting an appropriate resolution is crucial; important patterns in the footprint functions will be levelled out by using a resolution which is too coarse, whereas a resolution that is too high will give rise to outliers due to the statistical noise which is characteristic of stochastic simulations.This noise is particularly problematic in the BW model case, where the touchdown speed occurs in the nominator of flux and concentration functions (Eqs. 1 and 3).The effect of noise can be reduced by increasing the number of simulated particles but here we used particle numbers -on the order of 10 5 -which are typical for conventional models, while the number of particles in the LES was one order of magnitude more.
In the comparison among the three models, we used the LES results as the reference because of that model's most realistic description of the flow field and largest particle number, consequently producing the least noise.For model comparison it is important to compare both, sizes of the areas contributing to the signal as well as their locations.To concentrate on differences in shapes of the footprint patterns instead of differences in total contributions, all footprints were normalised by the total contribution from the selected horizontal footprint domain.Thus the total signal from the domain area was equal to unity.In order to take into account all the particles contributing to the signal, the footprints should be normalised by the contribution from an area that extends to infinity in each direction from the measurement point.As in practice an infinite domain is not computationally possible and moreover, it does not make sense in practice because no infinite fetch occurs in nature, we selected fixed footprint domain sizes for both stratifications.
In crosswind direction it was enough to extend the horizontal domain to 1000 m in each direction from the measurement location, while a distance of 200 m captured all of the contribution from downwind in all of the modelled cases.Because negative contributions from areas where descending particles outnumbered their ascending counterparts were observed in the LES results (see also Steinfeld et al., 2008), the contribution from the total domain is dependent on the domain upwind extension.The areas of negative and positive contribution are distributed symmetrically across mean wind direction, which is partly due to exclusion of the Coriolis force from the simulations (see Sect. 3.1 for effects of this exclusion).Thus, a crosswind integrated footprint shows the relative importance and along mean wind location of the areas of the negative contribution within the domain (Figs.1a  and 1b).Under the convective case 1 (Fig. 1a) the contri- bution of the negative part of the crosswind integrated footprint to the total signal within 5000 m upwind is at its highest at the highest measurement level, being 19% of the contribution of the positive part which dominates within 3000 m upwind from the measurement point.The influence of the negative part decreases monotonously towards low measurement heights and practically vanishes at the lowest height of 10 m.However, according to these results we set the up wind domain extension of case 1 to 3000 m for all measurement heights.Under the less unstable case 2 (Fig. 1b) there is practically no dominance of negative contributions observed within 5000 m upwind.As the crosswind integrated cumulative footprints reach a plateau at this distance the upwind domain extension was set to 5000 m.
In this study we normalise the footprints by the total contribution from the whole domain area of fixed size, i.e. both T. Markkanen et al.: Lagrangian footprint model comparison in 2-D positively and negatively contributing areas are included.Another option, of normalising by the fraction from positively contributing areas alone, would cause bias between the models because of the two fundamentally different sources of negative contributions, stochastic noise and overall flow pattern which dominates negative contributions of the BW and the LES models, respectively.While stochastic noise could be reduced by a coarse grid, the effect of flow pattern is better revealed by using a dense grid.In both cases the fraction of negative areas would decrease, whereas the flux contribution from the whole domain does not depend on the grid cell size and thus serves better as a normalisation parameter.Because the LES model showed most prominent areas of negative contributions to the signal, the selection of extension of total domain area presented above was consequently used for footprint normalisation for all models (Table 1).
We used the following approach to classify and compare the footprint functions.Firstly, we determined the peak locations in along-wind and crosswind directions.Secondly, to quantify the similarities in shapes and extensions of footprint patterns, we determined the smallest areas contributing certain percentages to the LES, BW and FW footprints.
Schmid and Oke (1990) used the term source area of level P ( P ) of the respective measure for analytical two dimensional footprint functions.As in the case of a stochastic model the footprints are given as the probability of each grid cell within the domain to serve as source for the measurements, we cannot determine a well defined single area bounded by an isopleth as is possible in the case of an analytical function.Instead, for some resolution -measurement height combinations, the smallest area contributing a given percentage of the total flux consists of several separate parts.Moreover, for certain combinations such a smallest area cannot be found at all as the contribution from an individual grid cell alone is higher than the given limits.This problem can be avoided by appropriate grid cell size -measurement height combination selection, which was used as the only means to reduce the influence of stochastic noise.In this work we adopted the notation by Schmid and Oke (1990) regardless of the above-mentioned fundamental difference between the presentation of footprint areas in analytical functions and stochastic estimators.
To select the suitable grid cell size -measurement height combinations we firstly ruled out combinations producing obvious outliers in along-wind peak position patterns when plotted as functions of measurement height.This criterion removed outliers from predictions by all the models and set a lower limit of grid cell size ( x) as follows: x>0.4 z m , where z m is measurement height.Secondly, we determined an upper limit such that the shapes of footprint functions would not suffer from too drastic an areal averaging.This limit resulted in the setting of grid cell size as follows: x≤z m .Together these criteria preserved one or two grid cell sizes for each measurement height.
For the selected data we determined the smallest areas contributing 10%, 20%, 50% and 80% to the footprints, that is 10 , 20 , 50 and 80 respectively.Additionally we determined the area that is common for both the validated model (BW or FW) and the reference (LES), that is an intersection of the two P s ( Val P ∩ Ref P ) and we denote it as ∩ P .In order to compare the equality of predicted footprint functions it is possible to determine the signal predicted by both models that originate from ∩ P .When both these values are close or in agreement to the target percentage, the two models agree perfectly for that part.When one of the models is close to the given upper limit but the other one is smaller, the former one is completely disclosed by the latter one.
For practical model comparison ∩ P is of more relevance than equality of the size of the area of level P by the qualified model ( V al P ) and the reference ( Ref P ).Nevertheless, the footprint size also has an influence because a Ref P which is too large would falsely indicate a good agreement by enclosing Val P .Therefore only a combination of both parameters can provide a good measure for comparison.Accordingly, we finally present a classification of the level of model agreement with the reference.The classification is based both on agreement of sizes of the source areas and on the degree of their overlapping.The size agreement between the examined model and the reference is given as follows: P and the degree of overlapping as follows: 1− ∩ P / Ref P The final agreement class ranging from 0 to 3 (no agreement to good agreement) is consequently determined according to the decision table shown in Table 2.This method of classification was principally adopted from Rebmann et al. (2005) and in the updated version Göckede et al. (2008) who had developed a scheme to combine footprints with the land use data.They defined levels of the target land use which must be in the footprint area such as 60 or 80%.Then they checked for how many measurements these levels are fulfilled.Also the quality check of turbulent fluxes is on the first view an arbitrary combination of different parameters (Foken and Wichura, 1996;Foken et al., 2004).

Comparison between LES model versions with and without Coriolis force
As the Coriolis force (CF) was not considered in the basic LES runs, we estimated its influence by making an additional run with the CF.In this run the sub-grid scale resolution was 10 m.The boundary layer characteristics of the additional run also differed from those of the basic LES run (see Table 1).First we evaluate the influence of lower sub-grid scale (SGS) resolution of the runs with the CF.Steinfeld et al. (2008) showed the importance of LES grid resolution at lowest measurement height of 10 m: higher resolution produced footprint function peaks that were located at greater distances and were of lower intensities than those predicted with lower resolution.When all the measurement heights are considered, however, the resolutions of 10 and 2.5 m show high agreement of peak position and intensity, the coefficients of determination being r 2 = 0.94 and r 2 = 0.96, for along mean wind flux footprint peak position and peak intensity respectively.As expected, for total contribution to the flux and concentration footprints from the whole domain area the models agree perfectly, because large flow structures are responsible for transport on that size-scale.Thus, we can conclude that the coarse SGS resolution is acceptable in assessing the effect of the inclusion of the CF.
Next we evaluate the influence of the inclusion of the CF on the footprint predictions.As under the influence of the CF the mean wind direction turns clockwise from the surface towards the top of the ABL, in the comparison the mean wind directions of both model versions were set identical separately at each measurement height investigated, instead of applying a common reference mean wind direction for all through the ABL.
Comparison of flux and concentration footprints predicted by both model parameterisations revealed that, for selected data the flux footprint peak positions along-wind were 10% to 30% closer to the observation point when the CF was included.In crosswind, the peak positions were well centralised implying that the mean wind direction at the observation height mostly determines the position of the footprint peak even when the CF is included.However, in the case of concentration footprints the peak positions were slightly biased towards counter-clockwise positions from the mean wind when the CF was included.The discrepancy between flux and concentration footprints is due to the fact that concentration footprints are positively contributed to by particles travelling both upwards and downwards (Eq.5).Accordingly, the contributing particles originate from greater distances and experience the turning wind in the ABL for longer time than those contributing positively to flux footprints.
The 10 , 20 , 50 and 80 of the model with the CF, are at the lowest observation level up to threefold greater compared to those of the results of the run without it (Fig. 2).The difference gets smaller towards larger heights and is opposite at the highest levels.At the lowest level the overlapping area among the two model versions is nearly equal to the case without the CF, which implies that the areas predicted by this model are enclosed by its counterpart taking the CF into account.This is largely due to the narrower peak resulting from the lower resolution of the version considering the CF, as reported in the beginning of this section and in more detail by Steinfeld et al (2008).At higher levels up to 250 m the footprint areas are quite similar for all effect levels, however the fractions of footprint areas common to both versions ( P s) get lower, which in turn is due to both the above mentioned discrepancy between peak positions and the distortion of symmetry across mean wind axis at far distances when the CF is taken into account.Behaviour of concentration footprint areas are very similar to those of flux footprints (not shown).Moreover, total areas contributing certain fractions to the concentration footprints are from 3 up to 40 -fold compared to those of respective flux footprints (not shown).These fractions are very similar among the two LES versions discussed above and are in agreement with earlier works (Schmid, 1994;Rannik et al., 2000) where the fraction was concluded to be an order of magnitude.
Finally, the overall model agreement classification was applied to this pair of model parameterisations using the version with the CF as reference (Fig. 3). 10 is generally of lowest agreement while agreement gets better and is either high or moderate for larger effect levels.The only exception to the general rule is the lowest measurement level where there is practically no agreement observed, which is obvious because of the reasons discussed above.We would like to point out that as there is no CF, there is no pressure gradient either in the simulations, which leads to decaying flow that in turn may influence the results at large distances.

Comparison of LES, BW and FW models
In Fig. 4 the flux footprint peak positions in along-wind direction as a function of measurement height are shown.Each model predicts relatively linear dependence between the two variables under both stability cases.Both FW and BW LS model results show very good agreement with LES predictions for measurement heights up to 100 m; the peaks predicted by BW model being located slightly further upwind for z m >50 m.This trend is even stronger, i.e., the BW model's peak location clearly further upwind for z m >100 m; while for these measurement heights, the FW model is not valid anymore.In cross wind direction the peak positions are located close to zero as expected (not shown).However, in the less unstable case 2 the LES results show a slight tendency towards positive values, that is, peak locations to the right from the mean wind axis, whereas in the case 1 the peak is located slightly to the left of the axis.This implies non-vanishing crosswind flow of particles, which may be due to discontinuous particle release in puffs of particles, instead of continuous release.
In the concentration footprint case the peak positions increase linearly in the along-wind direction with increas-  ing measurement heights as well, ranging from zero up to 2500 m (not shown).The crosswind positions showed similar asymmetry in the LES case as the flux footprints, while BW and FW model predictions were symmetric across mean wind axis.
In  predicted fluxes close to the respective P 's and mostly larger than those predicted by the LES.This implies that Val P is largely disclosed by Ref P , which in turn is an implication of a more concentrated BW model footprint.The fluxes are closer to each other under convective (Fig. 5) than under less unstable (Fig. 7) stratification.At high measurement levels the fluxes predicted by both models are very close to each other but they deviate from the target percentage more towards higher measurement levels.The ratio between pre- dicted fluxes from ∩ P in the FW model and the LES model comparison (Figs. 6 and 8) is opposite to that from the comparison between the BW and the LES.Generally the fluxes are of closest agreement when effect levels of 80% are observed.

Quality classification
From Fig. 9 we see how the features discussed above are reflected by the agreement class qualification scheme.Qualification of the BW against the LES reveals generally quite good agreement.At 80 (Fig. 9b) agreement is high for most of the intermediate measurement heights under convective stratification.However, the same stratification also shows the worst results when agreement of 20 at high measurement levels (Fig. 9a) is observed.This behaviour is expected due to the deviation in peak positions between the BW and the LES models discussed above.The less unstable case 2 (Fig. 9c, d) is mostly of moderate agreement.
The FW model shows mostly moderate agreement against the LES (Fig. 9e to h).Last shown are the FW model predictions qualified against the BW (Fig. 9i to l).Here the agreement turns out to be mostly high and especially good under less unstable stratification.
The agreement class qualification of concentration footprints (not shown) revealed mostly same or better agreement than that of flux footprints.The only cases where the agreement was reduced from their flux footprint counterparts, were the agreement of the FW model against the LES at the lowest measurement level in the L = −32 m and least agreement of concentration footprints was found in the L = −32 m and 20 case at measurement heights z m >200 m in the BW against the LES qualification.Again the obvious reason is the disagreement of peak positions among the two models, which in turn is probably due to the differences in the flow fields in both models.
The importance of comparison of the footprint functions in two dimensions can be clearly seen in Fig. 10 which shows cross wind and along wind integrated flux footprint functions predicted by the three models at the height of 50 m in stability case 2. Cross wind integrated footprint functions (Fig. 10a) are in relatively good agreement: peak positions agree perfectly and the agreement of peak intensities among the three models seems good too.However, the agreement class approach for 80 rates the BW and the FW models against the LES of no and low agreement, respectively.The reason for the bad rating can be seen in Fig. 10b, which reveals that the differences among the three models lie in cross wind distributions of the footprint predictions.

Conclusions
As routinely used stochastic footprint models lack means of validation (Foken and Leclerc, 2004) we facilitated LES footprint predictions in 2-D to explore the performance of two conventional models under two atmospheric stratifications over homogeneous surface.We present measures to quantify the model agreement against a reference model at different effect levels giving the extension of the area of the most important footprint around its peak position.Our purpose is not to draw conclusions of applicability of conventional models in practical applications such as that of Göckede et al. (2008) combining footprint synthesis with data quality analysis.For absolute evaluation of the performance of the models the covered parameter space of this study is too small.Nevertheless, this study demonstrates that when application of footprint functions facilitates both horizontal dimensions, the footprints need to be evaluated in two dimensions too, while comparison of alongwind agreement only may be misleading.
In comparison of conventional models against the LES according to the agreement classification scheme that takes into account both extension of predicted footprint area and degree of overlapping between the areas, we determined that the performance of two widely used models is rated very good or moderately good for most of the measurement heights.The agreement is generally better for intermediate measurement heights and for the convective case, whereas the two conventional flux footprint models agree best under less unstable conditions.This is not surprising as good agreement among the two models under neutral stratification was already reported by Kljun et al. (2002).Peak positions, that were evaluated additionally, turned out to be very similar at low measurement heights as influence of convective structures is not pronounced.The least agreement among the BW model and the LES was observed at large measurement heights, which implies differences in the flow field predicted by the LES, on the one hand, and the flow field effective in the BW model runs parameterised with four independent variables derived from LES, on the other hand.As the FW model is only parameterised for ASL, it was not compared at all for heights above 100 m.Additionally, the LES model was run including the Coriolis force (CF) to assess the influence of more realistic flow description.Comparison with the case without CF revealed pronounced differences at lowest observation height while at other levels, according to the quality categories set in Table 2, the models were of high or moderate agreement.
In this work we excluded the areas of pronounced negative fluxes from the analysis by restricting the domain within suitable along mean wind distances.However, under neutral and slightly stable stratifications and over a heterogeneous surface the negative areas may be located asymmetrically across mean wind direction (Steinfeld et al., 2008), which makes their exclusion more problematic.For those cases the applicability of the method has to be further considered.

Fig. 1a .Fig. 1b .
Fig.1a.Cumulative crosswind integrated flux footprints in the convective case 1 predicted by the LES model for three observation levels within 5000 m upwind from the measurement point.

Fig. 2 .
Fig. 2. The sizes of smallest source areas contributing (a) 10%, (b) 20%, (c) 50% and (d) 80% to the flux footprints, i.e. 10 , 20 , 50 and 80 .The model with Coriolis force indicated by circles, the model without Coriolis force with squares and the common area of both models ∩ 10 , ∩ 20 , ∩ 50 and ∩ 80 with asterisk.Results are only shown for selected grid resolutions (see Section 2.4 for criteria).

Fig. 3 .
Fig. 3. Agreement classes for validation of the flux footprints predicted by the LES parameterisation without Coriolis force against the parameterisation including Coriolis force in the convective case 1. Validations for (a) 10 , (b) 20 , (c) 50 and (d) 80 .Results shown only for selected grid resolutions (see Section 2.4 for criteria).

Fig. 4 .
Fig. 4. Along-wind peak position of the flux footprint as a function of measurement height for the LES (crosses), backward (circles) and forward (triangles) LS models (a) in the case 1 (convective) and (b) in the case 2 (less unstable).Results shown only for selected grid resolutions (see Sect. 2.4 for criteria).

Fig. 9 .
Fig. 9. Validation of the flux footprints predicted by the BW model against the LES results in the convective case 1 for (a) 20 , (b) 80 and in the less unstable case 2 for (c) 20 , (d) 80 .Validation of the flux footprints predicted by the FW model against the LES results in the convective case 1 for (e) 20 , (f) 80 and in the near less unstable case 2 for (g) 20 , (h) 80 .Validation of the flux footprints predicted by the FW model against BW results in the convective case 1 for (i) 20 , (j) 80 and in the near less unstable case 2 for (k) 20 , (l) 80 .RResults are only shown for selected grid resolutions (see Sect. 2.4 for criteria).

Fig. 10 .
Fig. 10.Flux footprint functions predicted by the three models at the height of 50 m in stability case 2. Footprint functions integrated (a) cross wind and (b) along-wind over the domain.
. (Case 1 with Coriolis force was only used for LES sensitivity run).

Table 2 .
Quality categories for footprint comparison.