AirWare User Manual

AERMOD: A Dispersion Model
for Industrial Source Applications

This documentation is extracted from Perry et al. (1998).

Algorithms for the Stable Boundary Layer

The form of the AERMOD concentration expression (C(x,y,z) in eq. 1) that is appropriate for stable conditions (L > 0) in flat terrain (and is similar to that of ISC2) is the Gaussian expression:

where

and

In equations 15-17, plume dilution and plume spread are calculated using effective boundary layer variables (see section on inhomogeneous boundary layer), h_p is the plume height, s _v and s _w are discussed below, and h_a is the level where vertical mixing is limited (= MAX(h_p,h)) where h is the stable (mechanical) mixed layer height. h_p is equal to the stack height plus plume rise (D h), where plume rise is given by Weil⁰ as

where x is downwind distance, N is the Brunt-Vaisala frequency, given by

and q is potential temperature. Wind speed, u, and N are evaluated at stack height (subject to change after development evaluation). The plume buoyancy flux, F_b, and momentum flux, F_m, (given by Weil⁰) are functions of stack parameters and ambient temperature. Equation (18) applies while the plume is rising. Weil⁰ defines the distance, x_max, to final rise.

The plume dispersion parameters, s _y and s _z, are given by the familiar formulas:

and

where t is plume travel time (x/U) and T_Lz is the Lagrangian integral time scale given by Venkatram⁰. Note that the values of s _v, s _w, U, and T_Lz used in equations (20) and (21) are effective values representative of the layer over which the plume is diffusing (see section on inhomogeneity).

Inhomogeneity in the Boundary Layer

The algorithms in AERMOD function under the assumption that the atmospheric boundary layer is vertically homogeneous (single values of meteorological variables or parameters represent the layer). Therefore, we designed a method to "convert" the inhomogeneous values (as measured or estimated) into equivalent (representative) homogeneous values. The approach basically interpolates (with downwind distance) between two observations of plume behavior: 1) near the source, plume dispersion is dominated by meteorological variables near the release point; 2) when the plume later disperses through the depth of the mixed layer, it is reasonable to assume that plume behavior is governed by meteorological variables averaged through the layer.

We will denote these two states for a generic meteorological variable, a (z), as a (h_p) for the value of the variable at release height (near the source) and

for the value of the variable averaged through the layer from the surface to h_a. h_a, the averaging height, is the maximum of the plume height, h_p, and the mixed layer height, z_i or h. To interpolate between a (h_p) and , AERMOD uses the following expression to find the "effective" value:

where x is downwind distance, is the arithematic average between u(h_p) and , and _Lz is the arithematic average of T_Lz(h_p) and _Lz. The appropriate formulations for T_Lz in stable⁰ and convective⁰ conditions are used for this averaging. Equation (23) approaches the required limits of a (h_p) and as x approaches zero and a large value, respectively.

When h_p is less than the mixing height, the relevant plume mixing is through that layer so the variables are averaged through the layer. When the plume is above the mixed layer, the variables are averaged through the mixed layer and into the layer above (up to h_p) where there is relatively little turbulence. In AERMOD, plumes that are in the elevated stable layer are allowed to disperse downwards and become entrained into the mixed layer. In these cases, near the source, the rate of spread is controlled by the small value of turbulence at the source height. AERMOD uses equation (23) to calculate the appropriate effective values of wind speed, s _v, and s _w as needed throughout the computations.

Meteorological Interface

Within AERMOD are routines that are designed to interface between the information provided by the meteorological preprocessor (see following section on AERMET) and the meteorological data needed by the model. Specifically, the interface routines take the computed (by AERMET) boundary layer parameters, the measured meteorological data, and other site-specific information to compute vertical profiles of:

wind direction,
wind speed,
temperature,
vertical potential temperature gradient,
vertical turbulence, s _w, and
6) horizontal turbulence, s _v.

For any one of these six variables (or parameters), the interface (in constructing the profile) takes each height above ground and if it is below the lowest measurement or above the highest measurement (or in some cases there is no data at all), the interface computes an appropriate value from selected profile relationships based on boundary layer parameterizations that are adjusted to the nearest data point, if any. If data are available both above and below a given height, an interpolation (based on both the measured data and the shape of the computed profile) is performed. This overall profiling approach simultaneously takes advantage of the information contained in both the measurements and parameterizations. As will be discussed, at least one level of measured wind speed, wind direction, and temperature is needed for the model to function. Turbulence can be parameterized without any direct turbulence measurements.

Wind Direction

At least one wind direction measurement is required for each simulation (hour) with AERMOD. Three methods are currently being considered for computing the wind direction change with height. The first is given by van Ulden and Holtslag⁰ as an exponential turning with height:

where WD(z) is the angle change from the surface wind direction, WD(h) is the angle change at h, and d₁ and d₂ are empirical constants equal to 1.58 and 1.0, respectively (NOTE: in convective conditions, z_i is substituted for h in equation 24). Van Ulden and Holtslag tabulated WD(h) as a function of Monin-Obukhov length, L. AERMOD interpolates between these values.

A second wind direction profiling method under consideration is that of Tennekes⁰ (currently implemented in CTDMPLUS) where the turning of the wind within the mixed layer is computed as a function of friction velocity (u_*), h, L, roughness length (z_o), and the Coriolis parameter. In either of methods one and two, the wind direction above the mixed layer is set equal to that at the top of the layer. The third option being considered for AERMOD is to have no wind direction oturning with height.

Wind Speed. At least one wind speed measurement is required for each simulation with AERMOD. As suggested by van Ulden and Holtslag⁰, the AERMOD profile equation for wind speed is:

where k is the von Karman constant, d is the displacement height, z_o is the roughness length, and y _m is a correction to the wind profile for stability other than neutral.

Temperature Gradient

. This section deals with the determination of vertical profiles of potential temperature gradient (dq /dz) from which AERMOD can also construct ambient temperature profiles. At least one air temperature measurement is required for each simulation with AERMOD.

In unstable conditions (L < 0), dq /dz is assumed to be zero from near the surface up to the top of the mixed layer, at z_i. Above z_i, dq /dz is assumed constant at 0.005 ° K/m.

In stable conditions (L > 0), AERMOD uses two profiling equations, one very near the surface and another above. Near the surface (below the lowest measurement of dq /dz or below 10 meters if no measurements) AERMOD uses an expression of Businger⁰ et al.

where q _*, the temperature scale, is a simple function of cloud cover⁰ (and is computed by AERMET). If one or more observations of dq /dz exist, then the value at the lowest level is used with eq. (26) to determine a value for q _*. Above the highest measurement level (or 10 meters with no measurements), AERMOD scales dq /dz with an expression by Stull⁰

where a is constant of 0.44. AERMOD uses the higher of eq. (27) and 0.002 ° K/m, wherever eq. (27) applies.

Turbulence in Convective Conditions. For L < 0, the profile of s _v is assumed constant from the surface to z_i. In the absence of measurements, s _v is calculated from⁰

where w_* is the convective scaling velocity. If observations of s _v are available, then AERMOD persists the value at the lowest level down to the surface and persists the value at the highest level in the mixed layer up to z_i. Between z_i and 1.2z_i, s _v is linearly decreased to 0.5 of its value at z_i.

In the absence of measurements, the profile of s _w is given by⁰

These formulations are complimented by any observations of s _w and the profile is adjusted as discussed earlier. As with s _v, s _w is linearly interpolated between z_i and 1.2z_i except that s _w reduces to 0.1 of its value at z_i.

Turbulence in Stable Conditions

. For L > 0, the vertical profile of s _v is given by⁰

where

The vertical profile of s _w in stable conditions has the same form as equations (31) and (32) with the substitution of s _wo for s _vo and with

If observations of s _v or s _w are available, they are used to adjust the profiles as discussed above.

AERMOD PREPROCESSORS

AERMOD has two data preprocessors that prepare meteorological data (on-site and off-site), terrain information, and receptor information for input to AERMOD. The meteorological preprocessor is named AERMET; the terrain/receptor preprocessor is yet unnamed.

AERMOD's Meteorological Preprocessor (AERMET)

As is the case with AERMOD, its meteorological preprocessor, AERMET, is based upon an existing regulatory model structure, that of the Meteorological Processor for Regulatory Models⁰(MPRM). AERMET is discussed by Weil⁰. Its overall treatment of the processing of meteorological data is similar to that done for the CTDMPLUS⁰ and HPDM⁰ models. The important surface parameters provided by AERMET are the Monin-Obukhov Length, L, surface friction velocity, u_*, surface roughness length, z_o, surface heat flux, H, and the convective scaling velocity, w_*. AERMET also provides estimates of the convective and mechanical mixed layer heights, z_i and h. In addition, AERMET passes through any measurements of wind, temperature, and turbulence (provided to it) in a form AERMOD expects.

The growth and structure of the atmospheric boundary layer is driven by the fluxes of heat and momentum which in turn depend upon surface effects. The depth of this layer and the dispersion of pollutants within it are influenced on a local scale by surface characteristics, such as the roughness of the underlying surface, the reflectivity (albedo), and the availability of surface moisture. For ISC2, default assumptions for these characteristics are built into the RAMMET/MPRM pre-processors, reflecting "open country" vegetation cover (a mixture of grassland with occasional forest).

For AERMET, the user can specify monthly variations of three surface characteristics for up to 12 upwind direction sectors. These include:

the albedo, which is the fraction of radiation reflected by the surface;
the Bowen ratio (B_o), which is the ratio of the sensible heat flux to the evaporation heat flux; and
the surface roughness length, z_o, which is the height above the ground at which the mean horizontal wind velocity is theoretically zero (z_o is also thought of as the average height of the surface roughness elements).

The user will be guided by look-up tables (in the AERMET user's guide) of typical values for these three variables for a variety of seasons and land use types.

Besides the specification of surface characteristics, AERMOD accepts meteorological data from the following sources:

standard hourly National Weather Service (NWS) surface data (winds, temperature, cloud cover, etc.) from the most representative site;
2) twice-daily soundings of winds, temperature, and dew point from the most representative NWS upper air station; and
3) on-site wind, temperature, turbulence, pressure, and radiation measurements.

Determination of the representativeness of data will vary with modeling situation (e.g. complex terrain applications may require on-site measurements).

Unstable Conditions

. When the net radiation is positive for a given hour (net radiation is either measured or estimated as a function of cloud cover, solar angle, albedo, and surface temperature), an energy balance technique is used to determine the sensible heat flux, H, as

where R_n is the net radiation, and B_o is the Bowen ratio. Eq. (35) assumes that the soil heat flux is 10% of the net radiation.

With this estimate of H, AERMET computes the surface friction velocity, u_*, and Monin-Obukhov length with an iterative approach (since each are functions of each other) similar to that of CTDMPLUS⁰. Friction velocity is found from eq. (25) and Monin-Obukhov length from eq. (6), described in previous sections. AERMOD initializes u_* assuming neutral conditions, calculates L, then proceeds with subsequent estimates of u_* and L until convergence (less than 1% change).

The boundary layer height (z_i) for the CBL is calculated with a simple one-dimensional energy balance model^0,0 using the early morning potential temperature sounding (prior to radiation sunrise) and the surface heat flux as a function of time. This model applies to an arbitrary temperature distribution which is obtained with early morning radiosonde soundings. Good agreement between predictions and observations of z_i have been found⁰.

Stable Conditions

. During stable conditions (negative net radiation) AERMET calculates the friction velocity with a method of Venkatram⁰ that is a counterpart to eq. (25) except that the stability corrections to the wind profile are not specifically functions of L. The energy balance approach is not used at night since the uncertainty in the ground heating term can be as large as the sensible and latent heat fluxes. Instead, the temperature scale⁰, q _*, is used by Venkatram in estimating u_*. The temperature scale is given by

where N_c is the opaque cloud cover. Once u_* and q _* are computed, the surface heat flux is given by

where r is air density. Finally, the Monin-Obukhov length is calculated from eq. (6) using u_* and H.

The boundary layer height (h) is calculated by AERMET in stable conditions with the diagnostic expression⁰

where f_c is the Coriolis parameter.

The Terrain Preprocessor

The terrain preprocessor for AERMOD is designed to provide the terrain information necessary to calculate the dividing streamline height, H_c, appropriate for each receptor location. The processor uses an objective method to estimate a "height scale", h_c, that is related to the terrain height in the vicinity of the receptor but does not correspond to the height of a particular hill. Instead, h_c can be thought of as the terrain height of the surrounding terrain that will influence the flow in stable conditions. Use of the height scale (instead of actual terrain feature height as in CTDMPLUS) to calculate H_c, provides a reasonable and objective method to calculate the weighting factor, f, in eq. (1).

In defining h_c for a given receptor, all terrain elevations within the user defined modeling domain and the distances of those elevations from the receptor are considered. Therefore each receptor may have a unique height scale. Consider a domain of interest, and a receptor at (x,y,z) for which we need an associated terrain height scale. The inherent assumption in this objective scheme is that 1) the effect of surrounding terrain on the flow near the receptor decreases with increasing distance and 2) the effect increases with increasing elevation of that terrain. In other words, the "effective elevation", h_eff, of surrounding terrain is a function of its actual elevation and its distance from the receptor. This is expressed quantitatively as:

where z_t is the terrain elevation at (x_t,y_t), r is the horizontal distance from the receptor to the terrain such that:

f_t(r/r_o) is the functional form of the distance weighting, and r_o represents the radius of influence and is taken to be:

where h_max is the elevation of the highest terrain in the modeling domain, and a is a constant set to 10. The distance weighting function is currently specified as:

After evaluation (and with further experience) we may find it appropriate to refine either r_o or the distance weighting function.

For a given receptor, h_eff is calculated for all terrain points within the modeling domain. This is why it is imperative to have the terrain information already digitized. The height scale for receptor is then related to the maximum effective value. Specifically, h_c is the actual terrain elevation at the location with the maximum h_eff. Quantitatively, the height scale is given as:

This calculation is repeated for each receptor location requested by the user. The terrain preprocessor is designed to include a receptor generator in which the user can specify receptor grids (polar, rectangular, discrete), as is currently done in ISC2, and the preprocessor will assign the appropriate terrain elevation to each receptor and the calculate the height scale for each receptor.

This objective approach to terrain preprocessing is designed for use with digitized terrain data such as the Digital Elevation Model (DEM) data available from the U. S. Geological Survey. As digitized data becomes more and more readily available, the AERMOD objective approach will eliminate the subjectivity and inconsistency in terrain specification. In addition, it is very user friendly in that receptor and terrain information can be provided to the dispersion model without extracting detailed information from terrain maps.

AERMOD MODEL EVALUATION

A critical element of any model development effort is the evaluation of the model algorithms with data from high-grade measurement studies. AERMIC is approaching model evaluation in two distinct phases: 1) model development evaluation where we are interested primarily in whether the model is based on sound physical principles and yields good results for the "right" reasons, and 2) an operational performance evaluation where the focus is on how well the model predicts concentrations on the high end of the distribution (which is important to regulatory decision making). As the model coding is nearing completion (as of this writing), AERMIC is beginning to undertake the model development evaluation (phase 1).

The are several reasons to perform a thorough developmental evaluation rather than skipping directly to "how well does the model perform?" Establishing that the model is based on good model physics greatly increases a user's confidence in applying the model in situations (sources, terrain, climatology) for which there is no existing evaluation data. In addition, the developmental evaluation will also aid in flushing out pathologies in the model algorithms.

AERMIC has gathered field (meteorological and concentration) data from eight separate studies for the phase 1 (developmental) evaluation of AERMOD and AERMET. These include tracer studies of rural, flat-terrain, surface releases, rural, flat-terrain, elevated-buoyant releases, and urban, flat-terrain, elevated releases; and sulfur dioxide studies of rural, flat-terrain, elevated releases, and rural, complex-terrain, elevated releases. With this variety of data, we intend to challenge the AERMOD algorithms in as many situations as our time will allow.

Following the lead of Weil et al.⁰ and Hanna⁰ , the phase 1 analyses will focus on model residuals (predicted/observed) as functions of key model parameters. This will help to identify both systematic problems in model formulations (signified by trends in the residuals with parameter variation) and sensitivity to a particular parameter (signified by large variations in residuals). Residual plots are done with a time and downwind-distance (modified space) pairing of the data (using maximum predicted and observed on a sampling arc for each time period).

Particularly for the full-year sulfur-dioxide data bases, phase 1 will also include some statistical summaries of performance measures (e.g. fractional bias, mean square error) as suggested by Cox and Tikvart⁰ and also by Hanna⁰. Correlation plots of predictions to observations of concentration will also be helpful in comparing AERMOD to other models such as ISC2 and CTDMPLUS. Finally, it is AERMIC's intention to share the results of the phase one evaluation (as they exist at that time) at both the conference this summer and in future publications about this model.

Model Sensitivity

We can gain additional understanding of model behavior by performing sensitivity analyses on AERMOD. Two types of sensitivity analyses will be performed: 1) to determine the effects (on calculated concentrations) of incremental changes in various input variables, and 2) to determine the effects of extreme values of various input variables on concentrations. The first type will be helpful both in defining the needed accuracy for various input values and in identifying spurious behavior of the model (determining model robustness). The second type of analysis will be helpful in defining the limits of applicability of the model.

FUTURE PLANS FOR AERMOD

The initial impetus for the development of AERMOD did not come from a desire to develop a new state-of-art model. AERMOD was conceived from a growing sense that EPA's basic regulatory modeling approach could be improved scientifically, without compromising the important need for consistency. As such, the success of AERMOD will be based as much on its regulatory status as on its ability to successfully predict design concentrations. If AERMOD is found to accurately estimate concentrations but does not become a significant part of EPA's regulatory program, we would not have accomplished our primary goal. Consequently, AERMIC's plans for AERMOD go beyond model development and testing. We have planned activities that consider what is needed for both regulatory promulgation and technology transfer.

Promulgation of a regulatory model, i.e. achieving "Appendix A" guideline status⁰, is a long and arduous task. In preparation for the rule-making process, information will need to be developed that identifies the proposed model and how it will affect existing emission limits. Also, the guidance needed for consistent implementation of the model, must be available for public review. In support of these needs, AERMIC has planned for the development of both meteorological guidance and the comparisons of AERMOD with any existing regulatory model appropriate for the same applications.

By using concepts such as similarity scaling and non-Gaussian pdf's, AERMOD represents a significant change in the underlying science of regulatory models. For those environmental engineers and scientists whose modeling experiences have included only regulatory type analyses, the shift from ISC to AERMOD may be substantial. In anticipation of the need for training, we are developing a document, the AERMOD Model Formulation Document(MFD) with this transition in mind. The MFD will include not only a complete description of AERMOD, but also a primer on boundary layer meteorology as it relates to this new model. The MFD will also include explainations designed to make the link between PG dispersion and the approach taken by AERMOD. Also, we will be evaluating the need for other mechanisms (e.g. workshops) of transferring this technology. Finally, AERMIC anticipates that AERMOD will be presented to the public as a regulatory proposal upon completion of model development, testing, and evaluation.

REFERENCES

1. Weil, J. C.,"Updating the ISC model through AERMIC," in Proceedings 85th Annual Meeting of Air and Waste Management Association, 92-100.11, Air and Waste Management Association, Pittsburgh, PA, 1992.

2. U. S. Environmental Protection Agency, User's Guide for the Industrial Source Complex (ISC2) Dispersion Models. Vols 1-3, EPA-450/4-92-008a-c, U. S. Environmental Protection Agency, Research Triangle Park, NC, 1992.

3. Willis, G. E., and J. W. Deardorff,"A laboratory model of diffusion into a convective planetary boundary layer", Quart. J. Roy. Meteor. Soc., 102, 1976, 427-445.

4. Briggs, G. A.,"Some recent analyses of plume rise observations", in Proceedings of the Second International Clean Air Congress, H. M. Englund and W. T. Berry, Eds., Academic Press, New York, 1971, 1029-1032.

5. Briggs, G.A.,"Plume rise and buoyancy effects," in Atmospheric Science and Power Production, D. Randerson, Ed., DOE/TIC-27601, U.S. Dept. of Energy, 1984, 327-366.

6. Berkowicz, R., J. R. Olesen, and U. Torp,"The Danish Gaussian air pollution model (OLM): Description, test and sensitivity analysis, in view of regulatory applications," Air Pollution Modeling and Its Application, V. C. De Wispelaire, F. A. Schiermeier, and N. V. Gillani, Eds. Plemum, New York, 1986, pp 453-481.

7. U. S. Environmental Protection Agency,"Supplement B to the guideline on air quality modeling (revised)," EPA-450/2-78-027R-B, U. S. Environmental Protection Agency, Research Triangle Park, NC, 1993.

8. Paine, R. J. and B. A. Egan,"User's guide to the Rough Terrain Diffusion Model (RTDM)-Rev. 3.20 ERT Document PD-535-585, ENSR, Acton, MA, 1987, 260pp.

9. Perry, S. G.,"CTDMPLUS: A dispersion model for sources in complex topography. Part I: Technical formulations," J. Appl. Meteor., 31:633-645 (1992).

10. Sheppard, P. A.,"Airflow over mountains," Quart. J. Roy. Meteor. Soc., 82:528-529 (1956).

11. Snyder, W. H., R. S. Thompson, R. E. Eskridge, R. E. Lawson, I. P. Castro, J. T. Lee, J. C. R. Hunt, and Y. Ogawa,"The structure of the strongly stratified flow over hills: Dividing streamline concept," J. Fluid Mech., 152:249-288 (1985).

12. Gifford, F.A.,"Statistical properties of a fluctuating plume dispersion model," Advances in Geophysics, 6, Academic Press, 1959, 117-138.

13. Misra, P.K.,"Dispersion of nonbuoyant particles inside a convective boundary layer," Atmos. Environ., 16, 1982, 239-243.

14. Venkatram, A.,"On dispersion in the convective boundary layer," Atmos. Environ., 17, 1983, 529-533.

15. Weil, J.C.,"Dispersion in the convective boundary layer," in Lectures on Air Pollution Modeling, A. Venkatram and J.C. Wyngaard, Eds., Amer. Meteor. Soc., Boston, 1988, 167-227.

16. Weil, J.C., L.A. Corio, and R.P. Brower,"Dispersion of buoyant plumes in the convective boundary layer," Preprints 5th Joint Conference on Applications of Air Pollution Meteorology, Amer. Meteor. Soc., Boston, 1986, 335-338.

17. Weil, J.C., L.A. Corio, and R. P. Brower, 1993: A dispersion model for highly buoyant plumes in convective conditions. Proceedings at the 86th Annual Meeting of the Air and Waste Management Association, 93-TA-26A.06, Air and Waste Management Assoc., Denver, CO, 1993.

18. Baerentsen, J.H., and R. Berkowicz,"Monte Carlo simulation of plume dispersion in the convective boundary layer," Atmos. Environ., 18, 1984, 701-712.

19. Weil, J.C.,"A diagnosis of the asymmetry in top-down and bottom-up diffusion using a Lagrangian stochastic model," J. Atmos. Sci., 47, 1990, 501-515.

20. Venkatram, A.,"Dispersion in the stable boundary layer," in Lectures on Air Pollution Modeling, A. Venkatram and J.C. Wyngaard, Eds., Amer. Meteor. Soc., Boston, 1988, 229-2659.

21. Van Ulden, A. P., and A. M. Holtslag,"Estimation of atmospheric boundary layer parameters for diffusion applications," J. Clim. and Appl. Meteor., 24, 1985, 1196-1207.

22. Tennekes, H.,"Similarity relations, scaling laws, and spectral dynamics," in Atmospheric Turbulence and Air Pollution Modeling, Nieuwstadt and van Ulden, Eds., D. Reidel Publishing Co., Boston, MA., 1981.

23. Stull, R. B., An Introduction to Boundary Layer Meteorology, Kluwer Academic Publishers, Boston, MA, 1988, 666pp.

24. Hanna, S. R., and R. J. Paine,"Hybrid Plume Dispersion Model (HPDM) development and evaluation," J. Appl. Meteor., 28, 1989, 206-224.

25. Irwin, J. S,, J. O. Paumier, and R. W. Brode, Meteorological Processor for Regulatory Models (MPRM) User's Guide, EPA-600/3-88-043, U. S. Environmental Protection Agency, Research Triangle Park, NC, 1993.

26. Carson, D. J.,"The development of a dry inversion-capped convectively unstable boundary layer," Quart. J. Roy. Meteor. Soc., 99, 1973, 450-467.

27. Weil, J. C., and R. P. Brower, Estimating convective boundary layer parameters for diffusion applications, Report No. PPSP-MP-48, Maryland Power Plant Siting Program, Maryland Department of Natural Resources, Baltimore, MD, 1983, 45pp.

28. Venkatram, A.,"Estimating the Monin-Obukhov length in the stable boundary layer for dispersion calculations," Bound. Layer Meteor., 19, 1980, 481-485.

29. Niewstadt, F. T. M.,"The steady-state height and resistance laws of the nocturnal boundary layer: Theory compared with Cabauw observations," Bound. Layer Meteor., 20, 1981, 3-17.

30. Weil, J. C., R. I. Sykes, and A. Venkatram,"Evaluation air-quality models: review and outlook," J. Appl. Meteor., 31, 1992, 1121-1145.

31. Hanna, S. R.,"Uncertainties in air quality simulation model predictions," Bound. Layer Meteor., 62, 1993, 3-20.

32. Cox, W. M., and J. A. Tikvart, "A statistical procedure for determining the best performing air quality simulation model," Atmos. Environ., 21a, 1991, 2387-2395.