Lagrangian Stochastic Simulation of Atmospheric Dispersion from Sources near the Ground

A simple Lagrangian stochastic (LS) dispersion model was used to investigate atmospheric dispersion in the lower atmospheric boundary layer, for near-surface sources. The sensitivity of ground-level concentration to numerical errors was revealed. Numerical treatments of inhomogeneous atmospheric turbulence and surface reflection in Lagrangian trajectory calculations may cause non-negligible overestimates of ground-level concentration for elevated sources, although these errors become trivial for surface sources and for concentrations above the ground-level. Using a numerical error-correction scheme, the LS model was evaluated against an analytical model. Close agreement was found between the two models in predicting ground-level concentrations for dispersions in the surface layer. LS simulations for an elevated source were also conducted. A scaling scheme was proposed to normalize the dispersion results, by including the source height as a scaling length. The relationships between the normalized surface concentration and downwind distance were distinguished by atmospheric stabilities. In the near-field, the distance of peak ground-level concentration was 0.5, 1.0, and 2.0 times zs U/u*, where zs is the source height, U is the mean wind speed at height 10 m, and u* is the friction velocity, for unstable, neutral and stable atmospheric stability conditions, respectively. In the far-field, the concentration approached approximately the “-3/2”, “-1” and “-2/3” law for unstable, neutral and stable atmospheres respectively.


INTRODUCTION
Turbulence in the atmospheric boundary layer (ABL) is vertically inhomogeneous, especially in the surface layer where strong wind shear occurs because of the influence of the Earth's surface.As a result, modeling near-surface dispersion is still challenging, especially if the complexities of the source, the near-field effects, or variations in surface conditions, need to be considered.
Lagrangian stochastic (LS) modeling has proved to be simple and flexible for describing atmospheric dispersion (Shadwick et al. 2007), and is also able to model inhomogeneous surface and complex meteorological conditions (Schmid 2002).After the remarkable work of Thomson (1987), this method became more robust theoretically.Wilson and Sawford (1996) reviewed the theoretical basis and advantages of the LS method and surveyed its applications for dispersion in different types of atmospheric turbulence.
In recent years, the LS method has been extensively used in practical applications.For example, Anfossi et al. (2010) simulated dense gas dispersion using a LS model by considering the processes of "negative buoyancy" and "gravity spreading".Wilson et al. (2009) simulated dispersion with a 3-dimensional LS model in a very complex, urbanlike environment.In addition to treating forward dispersion problems, the LS model has also been used for the inverse dispersion problem to estimate source strength.Hsieh et al. (2003) used a 2-dimensional LS dispersion model to estimate the sources and sinks of sensible heat and carbon dioxide in a forest, and agreements were found between the predicted model and experimental measurements.Flesch et al. (2004) inferred the emission strength of a homogeneous source of well-defined area from a single concentration measurement by using a backward LS dispersion model, and found satisfactory accuracy for determining the source strength.
Another application of LS modeling involves footprint analysis of flux measurements (Kljun et al. 2002;Leclerc et al. 2003;Kljun et al. 2004;Vesala et al. 2008).Besides research on the practical application of LS modeling, efforts to improve the LS method are ongoing.For example, many studies addressed the issue of the lower boundary conditions in LS simulations and improvements were made to keep the well-mixed constraint on it.Wilson and Flesch (1993) showed that a perfect reflection scheme at the surface can satisfy the well-mixed constraint for Gaussian homogeneous turbulence.Thomson and Montgomery (1994) proposed a reflection scheme for non-Gaussian turbulence, but the reflection velocity could not be solved analytically until Anfossi et al. (1997) gave two different approximate analytical solutions.Furthermore, studies were conducted on the analysis of factors affecting trajectory calculations, such as spinning of the stochastic velocity vector (Wilson and Flesch 1997), surface delay effect (Wilson et al. 2001), instability of the LS modeling (Yee and Wilson 2007) and so on.There is clearly interest both in LS model applications and improvements in methodology.
Here we examine the performance of a simple LS model (Cai et al. 2008).First, the numerical error in the estimation of ground-level concentration was addressed.This error arises mostly from the temporal discretization of the LS model and the inhomogeneity of turbulence near the surface.Although this issue was noted by Wilson and Flesch (1997), they did not give a solution except to suggest shortening the time step.Cai et al. (2008) found that using a small time step in LS simulations could not sufficiently reduce this systematic error, and they proposed a numerical scheme to correct it.However, their work focused on the footprint concept and backward-in-time LS simulations results.The discretization error in the modeling of ground-level concentrations has not yet been investigated explicitly.Second, we compared the results of the corrected LS model to those of an analytical model, thereby enabling the performance of the corrected LS model to be assessed.Finally, we used this model to evaluate the near-field influence of an elevated source on ground-level concentration, which could not be reasonably resolved by the current analytical model.

The LS Model
The LS model adopted in this study is the forward-intime version in Cai et al. (2008), which is simplified from the models of Guo and Cai (2005), Cai et al. (2006), andCai andLeclerc (2007).This version of the model is actually three-dimensional, but assumes a passive scalar transported within an idealized ABL in a stationary flow along the xdirection only, i.e., U i = (U 1 , U 2 , U 3 ) = (U, V, W) = (U, 0, 0), where subscript i takes values of 1, 2, and 3 to represent the three components in Cartesian coordinates.Separating the deterministic part of the passive particle's motion from its stochastic part and using additional assumptions simplifies the model greatly, as described below.
The passive particle's mean displacement caused by the mean flow over a time step dt is written as: where U i is the mean wind profile.The stochastic movement of the particle caused by the turbulent fluctuation can also be written as: where u i is the stochastic velocity written in the stochastic differential equation as: where j g is a random number with zero mean, Gaussian distribution, and variance dt.The coefficients a i and b ij are to be determined, and subscript j also takes values of 1, 2, and 3 to represent the three components in Cartesian coordinates.Kolmogorov's similarity theory for the statistics of velocity increments over a small time interval dt suggests that: where f is the dissipation rate of turbulent kinetic energy, C 0 is a constant (C 0 = 3 in this study), and ij d is the Kronecker function.The coefficient a i satisfies the Fokker-Planck equation under the well-mixed constraint, but is not uniquely defined under multidimensional conditions (Thomson 1987).By assuming that the probability density function of the Eulerian turbulent velocity is Gaussian, a particular solution for a i was given by Thomson (1987) and later adopted by Flesch et al. (1995).We further assume horizontal homogeneity and independence of each of the three components of turbulent velocities from one another, yielding: where ( , , ) is the standard deviation of the turbulent velocities.Note that the subscript i does not sum up in Eq. ( 5).This form of a i treats the inhomogeneity only for the vertical component of velocity, and only in the vertical direction.

Parameters for the Surface and Upper Layers
This LS model can work over the whole depth of the ABL with various atmospheric stabilities, from stable to unstable stratification.However, this model focuses on surface-layer dispersion in this paper.The mean wind profile was prescribed following the flux-profile relationship given by Businger et al. (1971).Turbulent velocity variances were adopted following Hanna (1982) in the depth of the ABL.The turbulence dissipation rate in the ABL is parameterized as described by Pasquill and Smith (1983).Details of the formulae are provided in the Appendix.
Input parameters included the mean wind profile, profiles of turbulent velocity variances, turbulent dissipation rate, ABL height and source terms.In addition, the surface roughness also needed to be prescribed.To be consistent with an analytical model used for comparison, the mean wind profile and turbulent dissipation rate took the surface layer form, extrapolated to the entire ABL depth.However, a realistic parameterization of turbulent velocity variances in the ABL was employed to correctly incorporate the influence of turbulence in the upper layers, because Lagrangian particles may spread out of the surface layer, particularly for stable cases with a shallow surface layer.Thus, except for the ABL depth, which should be given as an external parameter, all other parameters for the LS simulation could be determined using the surface-layer parameters as the friction velocity u * and Obukhov length L.

Simulation and Implementation Method
We performed 21 simulation runs for both surface and elevated sources.The simulation conditions covered six surface types with varying roughness lengths z 0 .There were at least three runs with atmospheric stability from unstable to stable conditions for each surface type.According to Wieringa (1993), these six z 0 values represent surfaces from 'short grass and moss' to 'regularly-built large towns', respectively.Table 1 provides detailed information on all the simulation runs.
In the LS model simulation, particles were released uniformly and continuously from a finite volume 0.01 × 0.01 × 0.01 m at the surface (roughness height z 0 ) for surface sources and at a height of 14.6 m for elevated sources.A random initial velocity was ascribed to each particle with a Gaussian probability density function and turbulent velocity variances of u 2 v , v 2 v , and w 2 v at the initial height.A perfect reflection boundary on the surface was specified.Particles were reflected upward at the surface, and the vertical velocity was reversed (Wilson and Flesch 1993).The touchdown velocity was used to calculate the groundlevel (roughness height z 0 ) concentration normalized by the source strength Q as follows: where w , j i T is the touchdown velocity (Flesch et al. 1995).Through touchdown velocity, the interaction time to the ground surface of a LS particle, which is proportional to surface concentration, is calculated by Eq. ( 6).The subscripts denote the i th touchdown within the surface element (x k -Δx/2 ≤ x ≤ x k + Δx/2 and y l -Δy/2 ≤ y ≤ y l + Δy/2) by the j th particle.The grid position is given as (x k , y l ).N is the total released particle number and n is the total touchdown number of the j th particle within this surface element.
In total, 600000 particles were released under unstable or neutral conditions, and 900000 were released under stable conditions.The simulation domain kept a fixed depth of 1.5 km in the vertical direction for all runs, while the horizontal domain varied with atmospheric stability and z 0 from 12.5 × 12.5 to 80 × 80 km.Particles were released at the center of the horizontal domain.The horizontal resolution of the concentration calculation was 5 m under unstable conditions and 10 m under neutral and stable conditions.Lateral integration was applied to obtain a crosswind-integrated concentration for comparison with analytical results: where M is the lateral (y-direction) grid number of the LS model domain.

Numerical Error Correction
A LS model may introduce systematic errors into the concentration estimate if its trajectory calculation does not correctly deal with inhomogeneous turbulence.The strongest vertical inhomogeneity of turbulence occurs near the ground, and thus simulations of ground-level concentration can be the most contaminated.Wilson and Flesch (1993) recommended small time steps as a way to deal with this problem.Cai et al. (2008) found this error non-negligible in their footprint results with a backward-in-time LS simulation.However, in principle, their result could not quanti-tatively represent a real dispersion problem, because of the asymmetry between the backward and forward dispersion in non-homogeneous turbulence.
These numerical errors occur because temporal discretization is used in the LS simulation: parameters at t are used to calculate the results at the next moment (t t D + ).As long as the turbulence in the vertical direction is inhomogeneous, the upward and downward movements of particles will have different time scales in one time step t D .By discretizing Eqs. ( 2) -( 5) directly, a calculation scheme for the LS particle trajectories (in the vertical direction) can be defined: where + are used to calculate upward and downward motions within a time step t D .As a result, a bias is inevitable.We adopted a 'test-adjustment' scheme to reduce this bias by seeking a better estimate of the time scale within a time step.We used an approximation 2 ( ) x directly at the n th time step.The procedure was as follows: (1) Calculate a set of 'test' increases of velocity Δw and displacements Δz with Eqs.( 8) and ( 9). ( 2) Determine an adjusted time scale xl, according to the new position anticipated by the test displacement Δz: (3) Get the 'adjusted' increase of Δw and Δz by replacing x with xl in Eqs. ( 8) and ( 9).In addition to the problem of inhomogeneous turbulence, the commonly used scheme of surface reflection for trajectory calculations may also produce errors by violating the well-mixed constraint.Wilson and Flesch (1993) showed that no universal reflection scheme is well-mixed for skewed or inhomogeneous turbulence.Only for homogeneous turbulence with a symmetric velocity probability density function can a perfect reflection scheme at the surface satisfy the well-mixed condition.Thus, according to Wilson and Flesch (1993), one way to treat this problem is to put a homogeneous 'buffer' layer on the ground.Here we set this 'buffer' layer in the model by specifying all turbulent parameters below 2z 0 with values as if they were at this height, i.e., ( ) (2 )

Sensitivity to Numerical Error
By using and deactivating the numerical error correction scheme in the LS model, we simulated the ground-level concentration from an elevated source near the surface.
Two sets of runs with different z 0 were conducted, with three different stabilities in each set (Run10 -Run12 and Run13 -Run15 in Table 1).A comparison of the simulated ground-level concentrations with and without the numerical error correction scheme is shown in Fig. 1.
The effect of the numerical error could be significant.and Run13 -Run15 in Table 1).Source height is 14.6 m.
The simulated concentration using the correction scheme was obviously lower than without it.The difference could be up to 20% at the ground-level concentration peak for the three runs of z 0 = 0.05 m (Panel 1), and up to 14% for cases of z 0 = 0.3 m (Panel 2).It appears that the numerical error diminishes with downwind distance, which is clearer in unstable and neutral conditions.The results shown in Fig. 1 demonstrate that the numerical error is non-negligible in forward LS modeling and should be corrected.However, we also found that for conditions of surface sources, these numerical errors were not significant.Additionally, when the concentrations were calculated at a level just above the ground (1 or 2 m), the numerical errors also diminished quickly (results not shown).In short, the numerical error is only important if a problem involves dispersion results at the ground.This echoes the conclusion of Wilson (2007) that a very small timestep in LS trajectory calculation was needed only when surface concentrations were calculated.All the simulation results described below have employed the error correction scheme.

Evaluation of Model Performance
We used van Ulden's (1978) analytical model as a reference to evaluate the LS model.van Ulden's model is an analytical solution to the two-dimensional advection diffusion equation, for a point source on the ground surface.Power laws for both wind speed and eddy diffusivity in the atmospheric surface layer were assumed in solving this problem.Based on the Monin-Obukhov similarity theory, the spatial distribution of plume concentration can be calculated from the analytical solution.There was only one empirical parameter in the analytical solution, i.e., the shape factor, which was determined by laboratory and field experiments.Since our concern in this study was the ground-level concentration, it was advantageous to use the analytical model to calculate the counterpart result at the exact ground level.Another advantage in comparing our relatively simple LS model with the rather idealized analytical model was that the factors affecting dispersion were relatively isolated, and thus the physical processes influencing the result could be more clearly revealed.
Analytically, the lateral-integrated ground-level concentration C y for a surface point source can be written in dimensionless form [van Ulden 1978, Eq. ( 19)]: where z and U are the mean height and mean wind speed of the dispersion plume, respectively.Following van Ulden's (1978) scaling, the simulation results of all runs in Table 1 were grouped into seven classes according to the stability parameter z 0 /L. Figure 2 compares the LS-simulated and analytical results for each stability class.It shows that the ground-level concentrations simulated by the LS model qualitatively agree with those of the analytical model at most of the dimensionless distances.The results approximately follow a "-3/2" power law falloff with downwind distance under unstable conditions, and a "-2/3" power falloff under stable conditions, as noted by van Ulden (1978).For neutral conditions, a "-1" power falloff was found.However, there were systematic underestimates from the LS model at shorter dimensionless distances, especially for very unstable conditions.The cause of these differences between the LS and analytical models may be that the set of turbulent parameters used in the LS model were not exactly consistent with that used in the analytical model.For example, in vertical dispersion in the surface layer, the vertical velocity variance w v , Lagrangian time scale x , and the turbulent dissipation rate f used in the LS model are as follows: where b w and C 0 are both experimental coefficients, and m { is the dimensionless wind shear, which is a universal function of the dimensionless stability parameter z/L.Combining Eqs. ( 12) -( 14), we obtained the effective diffusivity for far-field dispersion as follows: While in the analytical model for resolving the twodimensional advection diffusion equation: the turbulent diffusivity K z Eu can be written as: if Monin-Obukhov similarity applies in the surface layer.
h { is the dimensionless temperature gradient, which is also adopted to describe other scalars.Different forms of functions and experimental coefficients are used in Eqs. ( 15) and ( 17).Any experimental uncertainty between these two sets of turbulent parameters may cause the differences between dispersion results.There have been previous studies of the influence of relevant experimental coefficients (b w and C 0 ) in LS models, e.g., Du (1997), Cai et al. (2008), Wilson et al. (2009).Detailed discussion of this issue is beyond the scope of this paper.However, despite these uncertainties, the performance of the LS model at shorter distances can still be regarded as reasonable.However, the LS simulation results obviously deviated from those of the analytical model at longer dispersion distances, and under different stability conditions.The LS model gave higher concentrations after a certain dimensionless distance.We determined this critical distance x c by checking the deviation in the figure slopes.Then the mean plume height z c corresponding to x c could be calculated analytically (van Ulden 1978).Table 2 lists these results for all LS simulation runs.All z c cases were close to or greater than 10% of their corresponding ABL heights.Usually the lower 10% of the ABL is considered as the atmospheric surface layer wherein the analytical model of van Ulden (1978) applies.Hence the deviation of the LS simulation results from those of the analytical model could be interpreted as the influence of the practical profiles of turbulence veloc-ity variances, particularly in the w v profile, that were used in the LS model.Figure 3  ) (Flesch et al. 2004) for an unstable case.The "deviation height" ( / z h c ) in Table 2 was roughly indicated in Fig. 3. Clearly, this deviation does correspond to the height where the difference between these two w v profiles becomes significant.On the other hand, the ABL height applied in the LS model, which acts as an upper limit to vertical dispersion, may also influence the far-field dispersion.However, rather simple analytical models are widely used in practice, e.g., in footprint analysis (Schmid 2002).Hence it is useful to note their applicability limits with regard to dispersion distance, as well as their possible underestimates at longer distances.

Effect of Source Height Z s
A LS model is flexible in application to elevated sources.However, van Ulden's (1978)  restriction, a near-surface source with small elevation can also be treated, if a virtual dispersion distance upstream of the source is coupled to the analytical formulae (van Ulden 1978).Actually, the near-field dispersion for an elevated source cannot be presented correctly by analytical models.We used the LS model to explicitly simulate the near-field dispersion of an elevated source and then analyzed the related variation in ground-level concentration.
All runs in Table 1 were also simulated for an elevated source.The cases in Table 1 cover different surface-roughness heights, atmospheric stabilities, and ABL heights.We believe that the simulations well represent actual conditions of near-surface dispersion, even though a single source height of 14.6 m was used.Figures 4a -c show the results of unstable, stable and neutral runs, respectively, compared directly with their counterpart results for a surface source.
The figure clearly shows that the concentration of an elevated source quickly approaches that of its counterpart for a surface source.The distance where the concentration difference between a surface source and an elevated source tends to vanish is of interest, without considering the virtual dispersion distance added to the surface source as done by van Ulden (1978).In practice, we accepted a relative difference within 5% as the distance where the influence of an elevated source vanishes, and found that the distance varied with atmospheric stability.It could be smaller than 100 m in a very unstable condition (Run04, with L = -2 m) and as large as 8000 m in a very stable condition (Run09, with L = 20 m).Mean wind speed also clearly affects this distance, with stronger wind tending to have a longer distance.In addition, when compared with runs having similar stability conditions (e.g., Run01, Run10, and Run 16; Run03, Run12, and Run18), we found that the roughness length z 0 also has a significant effect on this distance.Larger z 0 was associated with stronger turbulence, and thus resulted in a smaller distance.Table 3 summarizes the distance results estimated for all the runs and their corresponding mean plume height, as estimated by the analytical model of van Ulden (1978).On average, for a downwind distance where z z 2 s $ , the influence of source height can be ignored.This is consistent with the findings of van Ulden (1978).
To better describe the near-field dispersion, proper scaling is needed.For a surface source or a negligible source height, dispersion results were plotted as C y u * z 0 /Q versus x/z 0 (van Ulden 1978, and Fig. 2), or as C y u * |L|/Q versus x/|L|.The former takes only the roughness height into account.The latter considered the effect of atmospheric stability, but it cannot address near-neutral cases with Obukhov length L approaching infinity.Both of them did not take source height into consideration.
If the height of a source cannot be negligible and near-field dispersion is of interest, the source height should be a crucial length scale (Davis 1983, Wilson et al. 2000).By analogy with Willis and Deardorff (1978) who took source height as a scaling parameter, we tentatively propose a dimensionless scaling scheme of C y Uz s /Q versus xu * /z s U, where U indicates the near-surface wind speed at a reference height of 10 m, and z s represents the source height.The dimensionless downwind distance can also be rewritten as a dimensionless travel time of the plume, i.e., (x/U)/(z s /u * ), which can be interpreted as the horizontal travel time (x/U) of the plume scaled by the dispersion time (z s /u * ) of the plume from the source height to the earth's surface.Based on this physical consideration, this scaling scheme is expected to more reasonably apply to near-field dispersion.
Figure 5 shows the scaled results of all 21 simulation runs with an elevated source.In Fig. 5a, all the concentration points can be divided into three general groups: stable, neutral, and unstable, and their peaks almost coincide.The neutral runs show the best performance under this scaling scheme.Much scatter of data could be seen in both the stable and unstable cases.Two extreme runs of very stable and very unstable conditions are also notable (Run09 and Run04, respectively).Their results clearly depart from the majority of results of their respective groups.Figure 5b gives the geometric mean results for stable, neutral, and unstable runs, respectively, noting that the two extreme runs were excluded from the geometric average calculation.Although there is some fluctuation, there is a clear systematic shift of the curves from unstable to stable conditions.The dimensionless concentration peaks are located at dimensionless distances of about 0.5, 1.0, and 2.0 for unstable, neutral, and stable cases, respectively.This result provides a practical estimate for the distance of the peak ground-level concentration for an elevated source: about 0.5 -2 times z s U/u * .
Figure 5b shows the falloff of concentration with an increase in distance.Three ray-lines were plotted as references to the "-2/3", "-1", and "-3/2" power law falloffs, which were regarded as typical dispersion relationships of ground-level concentration and downwind distances (van Ulden 1978;Du and Venkatram 1997).However, the mean result shows slower falloff, regardless of stability type, at farther distances.

CONCLUSIONS
Dispersion of near-surface sources was simulated by a forward-in-time LS model, and the resulting ground-level concentrations were investigated in detail.Efforts were made to improve model prediction of concentrations at the ground surface.The LS model with a numerical error correction scheme was evaluated by comparison with results from the analytical model of van Ulden (1978).Near-field dispersion of an elevated source was simulated and a scaling scheme was proposed to describe the influence of the source height.
The numerical errors arising from treatment of both inhomogeneous atmospheric turbulence and the surface-reflection were not negligible.These errors led to systematic overestimations of ground-level concentrations, with bias of as much as 20% of the ground-level peak concentration for a source elevated at 14.6 m.However, these errors were significant only for elevated sources and ground-level concentration; for surface sources or for concentration just above the ground of 1 or 2 meters, these errors became trivial.A 'test-adjustment' scheme and a homogeneous buffer layer on the ground were used to amend the numerical errors.The LS model with the correction schemes showed good performance in simulating ground-level concentrations.General agreement was found between the LS model results and those of the analytical model (van Ulden 1978) in 21 runs with different surface roughness and atmospheric stability conditions.Deviation between the results farther downwind demonstrated the different parameterizations in the two models, especially for the w v profile.A scaling method was developed for the near-field dispersion of elevated sources.In this method, the dimensionless downwind distance is xu * /z s U and the dimensionless concentration is C y Uz s /Q, where x, C y , and Q are the downwind distance, crosswind integrated ground-level concentration, and source strength, respectively; u * , z s , and U represent the friction velocity, source height, and the wind speed at height 10m.Using this scheme, all dispersion data from the LS simulation runs for an elevated source were organized into three groups of stable, neutral, and unstable stratification, with similar behavior and almost coinciding concentration peaks.This scaling scheme appears to be useful for describing the general characteristics of near-surface dispersion for elevated sources.
the superscript n indicates the time step, and N(0, 1) obeys the standard normal distribution.The beginning turbulent time scales ( )

Fig. 2 .
Fig. 2. Dimensionless ground-level concentrations for seven classes of stability.The horizontal axis indicates the dimensionless downwind distance (x/z 0 ).Symbols show the LS simulation results and solid lines represent the results by van Ulden's (1978) analytical model.
compares the profile of vertical velocity variance used in our LS model ( _PBL w v ) with that in the atmospheric surface layer ( _ASL w v

Fig. 3 Fig. 4 .
Fig. 3 Comparisons between the profiles of vertical velocity variance used in our LS model _PBL w v and those in the atmospheric surface layer _ASL w v (Flesch et al. 2004).The meteorological condition corresponds to Run16.The solid line roughly indicates / z h c

Table 1 .
Parameters for the simulations.
Note: z 0 , surface roughness length; U 10 , mean wind speed at 10 m; Q s , surface heat flux; h, ABL height.

Table 2 .
analytical model in principle applies only to a surface point source.By relaxing this Critical distance x c and corresponding plume height z c for the deviation between LS simulation and the analytical results in Fig.2.