AERMOD Model Formulation and Evaluation

EPA-454/ R-18-003

April, 2018

AERMOD Model Formulation and Evaluation

U.S. Environmental Protection Agency

Office of Air Quality Planning and Standards

Air Quality Assessment Division

Research Triangle Park, NC
ii

Notice

The original document, AERMOD: Description of Model Formulation, was published in 

September, 2004, with the EPA document number EPA-454/R-03-004. This version of the 

document has several revisions, including a new title (AERMOD Model Formulation and 

Evaluation), the new model formulations added to the model since the original publication in 

2004, a model evaluation based on the latest version of the model (version 16216 of AERMOD 

and AERMET), and minor changes to the formatting and supplemental content to make it 

appropriate according to current EPA document publication standards. 

This report has been reviewed by the Office of Air Quality Planning and Standards, U.S. 

Environmental Protection Agency, and has been approved for publication. Mention of trade 

names or commercial products does not constitute endorsement or recommendation for use.
iii

Acknowledgments

This project was made possible through the continued support of Mr. Joe Tikvart, of 

EPA’s Office of Air Quality Planning and Standards (OAQPS), and Mr. Frank Schiermeier, of 

NOAA’s Atmospheric Sciences Modeling Division. The authors are particularly grateful to Dr. 

Gary Briggs, and Mr. John Irwin, of NOAA’s Atmospheric Sciences Modeling Division, for 

their thorough and constructive review of an earlier version of this document. Finally, we would 

like to thank the many scientists who participated in peer reviews and beta testing.

This report was compiled primarily by Alan J. Cimorelli, Steven G. Perry, Akula 

Venkatram, Jeffrey C. Weil, Robert J. Paine, Robert B. Wilson, Russell F. Lee, Warren D. 

Peters, Roger W. Brode, James O. Paumier, and James Thurman
iv

Table of Contents

Section Page

Notice.............................................................................................................................................. ii

Acknowledgments..........................................................................................................................iii

Table of Contents........................................................................................................................... iv

Figures............................................................................................................................................ ix

Tables............................................................................................................................................. xi

1. Introduction................................................................................................................................. 1

1.1 Background......................................................................................................................... 1

1.2 The AERMIC focus: a replacement for the ISC3 model.................................................... 3

1.3 Model development process ............................................................................................... 4

1.4 Purpose of the document..................................................................................................... 6

2. Model overview .......................................................................................................................... 7

3. Meteorological preprocessor (AERMET)................................................................................. 11

3.1 Energy balance in the PBL ............................................................................................... 11

3.1.1 Net radiation............................................................................................................. 12

3.1.2 Transition between the CBL and SBL..................................................................... 13

3.2 Derived parameters in the CBL ........................................................................................ 14

3.2.1 Friction velocity (u*) & Monin Obukhov length (L) in the CBL ............................ 14

3.2.2 Convective velocity scale (w*)................................................................................. 15

3.3 Derived parameters in the SBL......................................................................................... 16

3.3.1 Friction velocity (u*) in the SBL.............................................................................. 16

3.3.2 Sensible heat flux (H) in the SBL............................................................................ 19

3.3.3 Monin-Obukhov length (L) in the SBL ................................................................... 19

3.4 Mixing height.................................................................................................................... 19
v

3.4.1 Convective mixing height (zic)................................................................................. 20

3.4.2 Mechanical mixing height (zim) ............................................................................... 20

3.5 Adjustment for the low wind speed/stable conditions in AERMET................................. 22

4. Vertical structure of the PBL - AERMOD’S meteorological interface.................................... 24

4.1 General profiling equations............................................................................................... 25

4.1.1 Wind speed profiling................................................................................................ 25

4.1.2 Wind direction profiles............................................................................................ 27

4.1.3 Profiles of the potential temperature gradient.......................................................... 27

4.1.4 Potential temperature profiling ................................................................................ 30

4.1.5 Vertical turbulence calculated.................................................................................. 31

4.1.6 Lateral turbulence calculated by the interface ......................................................... 34

4.1.6.1 Mechanical portion of the lateral turbulence .................................................. 35

4.1.6.2 Convective portion of the lateral turbulence................................................... 37

4.2 Vertical inhomogeneity in the boundary later as treated by the interface ........................ 38

5. The AMS/EPA regulatory Model: AERMOD.......................................................................... 43

5.1 General structure of AERMOD including terrain............................................................. 45

5.2 Concentration predictions in the CBL .............................................................................. 50

5.2.1 Direct Source contribution to concentration calculations in the CBL..................... 59

5.2.2 Indirect Source contribution to concentration calculations in the CBL................... 61

5.2.3 Penetrated source contribution to concentration calculations in the CBL............... 62

5.3 Concentrations in the SBL................................................................................................ 62

5.4 Treatment of lateral plume meander................................................................................. 64

5.5 Estimation of dispersion coefficients................................................................................ 66

5.5.1 Dispersion from ambient turbulence........................................................................ 68

5.5.1.1 Lateral dispersion from ambient turbulence ................................................... 68
vi

5.5.1.2 Vertical dispersion from ambient turbulence.................................................. 70

5.5.2 Buoyancy induced dispersion (BID) component of σy and σz ................................. 73

5.5.3 Treatment of building downwash ............................................................................ 73

5.6 Plume rise calculations in AERMOD............................................................................... 75

5.6.1 Plume rise in the CBL.............................................................................................. 75

5.6.2 Plume rise in the SBL .............................................................................................. 77

5.7 Source characterization..................................................................................................... 78

5.8 Plume volume molar ration method (PVMRM)............................................................... 79

5.8.1 Definition of plume volume..................................................................................... 80

5.8.1.1 Total vs. relative dispersion ............................................................................ 80

5.8.1.2 Calculation of relative dispersion coefficients................................................ 80

5.8.1.3 Treatment of volume and area sources ........................................................... 82

5.8.1.4 Defining extent of plume ................................................................................ 83

5.8.1.5 Adaption for AERMOD terrain algorithm..................................................... 84

5.8.1.6 Treatment of penetrated plumes...................................................................... 84

5.8.2 Minimum ozone concentration for stable conditions............................................... 85

5.9 Adjustments for the urban boundary layer........................................................................ 85

6. List of symbols.......................................................................................................................... 90

Appendix A: Input / output needs and data usage ........................................................................ 95

A.1 AERMET input data needs.............................................................................................. 95

A.1.1 Meteorology............................................................................................................ 95

A.1.2 Directionally and/or Monthly Varying Surface Characteristics ............................. 96

A.1.3 Other........................................................................................................................ 96

A.1.4 Optional................................................................................................................... 96

A.2 Selection and use of measured winds, temperature, and turbulence in AERMET.......... 97
vii

A.2.1 Threshold Wind Speed............................................................................................ 97

A.2.2 Reference Temperature and Height ........................................................................ 97

A.2.3 Reference Wind Speed and Height ......................................................................... 97

A.2.4 Calculating the Potential Temperature Gradient above the Mixing Height from 

Sounding Data................................................................................................................... 98

A.2.5 Measured Turbulence.............................................................................................. 98

A.2.6 Data Substitution for Missing On-Site Data ........................................................... 98

A.3 Information passed by AERMET to AERMOD.............................................................. 99

A.4 Restrictions on the growth of the PBL height.................................................................. 99

A.5 Initializing the mechanical mixing height smoothing procedure..................................... 99

A.6 Determining the mixing height when the sounding is too shallow.................................. 99

A.7 Input data needs for AERMAP ...................................................................................... 100

A.8 Information passed by AERMAP to AERMOD............................................................ 100

A.9 Wind speed and turbulence limits used in model calculations ...................................... 100

A.10 Using profiles for interpolating between observations................................................ 101

A.11 Using measured mixing heights................................................................................... 103

Appendix B: Model evaluation results ....................................................................................... 104

B.1 Introduction .................................................................................................................... 104

B.2 Database descriptions..................................................................................................... 104

B.2.1 Martin’s Creek....................................................................................................... 106

B.2.2 Tracy...................................................................................................................... 108

B.2.3 Lovett .................................................................................................................... 109

B.2.4 Westvaco ............................................................................................................... 111

B.2.5 DAEC.................................................................................................................... 112

B.2.6 EOCR .................................................................................................................... 113
viii

B.2.7 Alaska.................................................................................................................... 114

B.2.8 Indianapolis........................................................................................................... 116

B.2.9 Kincaid .................................................................................................................. 119

B.2.10 AGA .................................................................................................................... 120

B.2.11 Millstone.............................................................................................................. 121

B.2.12 Bowline ............................................................................................................... 122

B.2.13 Baldwin ............................................................................................................... 122

B.2.14 Clifty Creek......................................................................................................... 124

B.3 Evaluation methodology ................................................................................................ 126

B.3.1 AERMET/AERMOD comparisons....................................................................... 126

B.3.2 Evaluation procedures................................................................................................. 126

B.3.2.1 Robust highest concentrations............................................................................ 126

B.3.2.2 EPA Protocol for determining best performing model ...................................... 127

B.4. Results........................................................................................................................... 130

B.4.1 Turbulence cases................................................................................................... 130

B.4.2 Non-turbulence cases ............................................................................................ 133

B.4.3 Statistical evaluations............................................................................................ 133

B.5 Summary/Conclusions ................................................................................................... 144

Appendix C: Evaluation of modified urban option..................................................................... 145

References....................................................................................................................................... 3
ix

Figures

Figure Page

Figure 1. Data flow in the AERMOD modeling system................................................................. 9

Figure 2. Wind speed profile, for both the CBL and SBL, in the region below 7z0 .................... 26

Figure 3. Wind speed profiling, for both the CBL and SBL, in the region above 7z0................. 27

Figure 4. Profile of potential temperature gradient for the SBL................................................... 30

Figure 5. Convective portion of the vertical turbulence in the CBL ............................................ 32

Figure 6. Mechanical portion of the vertical turbulence in the CBL............................................ 33

Figure 7. Profile of vertical turbulence in the SBL....................................................................... 34

Figure 8. Family of lateral mechanical turbulence profiles over a range of mechanical mixing 

heights........................................................................................................................................... 37

Figure 9. AERMOD’s Treatment of the Inhomogeneous Boundary Layer.................................. 39

Figure 10. AERMOD two state approach. The total concentration predicted by AERMOD is the 

weighted sum of the two extreme possible plume states.............................................................. 46

Figure 11. Treatment of Terrain in AERMOD. Construction of the weighting factor used in 

calculating total concentration ...................................................................................................... 49

Figure 12. Instantaneous and corresponding ensemble-averaged plume in the CBL................... 52

Figure 13. AERMOD’s three plume treatment of the CBL.......................................................... 53

Figure 14. AERMOD’s pdf approach for plume dispersion in the CBL...................................... 56

Figure 15. Probability density functon of the vertical velocity. While the Gaussian curve is 

unskewed, the bi-Gaussian curve has a skewness of S=1............................................................. 57

Figure 16. Lateral spread (Fy) as a function of non-dimensional distance (X). The data is taken 

from the Prairie Grass experiment (Barad, 1958)......................................................................... 69

Figure 17. AERMOD’s construction of a continuous meteorological profile by interpolating 

between observations.................................................................................................................. 102
x
xi

Tables

Table Page

Table 1. AERMOD evaluation databases used for comparisons of AERMOD 15181 and 

AERMOD 16216. ....................................................................................................................... 105

Table 2. Hourly, 3-hour, and 24-hour RHC for turbulence cases. Best performing model 

compared to observed RHC are highlighted in gray................................................................... 132

Table 3. Hourly, 3-hour, and 24-hour RHC for non-turbulence cases. Best performing model 

compared to observed RHC are highlighted in gray................................................................... 133

Table 4. Composite Performance Measure (CPM) for turbulence cases. Scenarios with lowest 

CPM’s for each study location are highlighted in gray. ............................................................. 134

Table 5. Composite Performance Measure (CPM) for non-turbulence databases. Scenarios with 

lowest CPM’s for each study location are highlighted in gray................................................... 134

Table 6. Martins Creek Model Comparison Measure (MCM) results. Confidence intervals 

highlighted in gray are significant at that percent....................................................................... 137

Table 7. Lovett Model Comparison Measure (MCM) results. Confidence intervals highlighted 

in gray are significant at that percent.......................................................................................... 138

Table 8. Westvaco Model Comparison Measure (MCM) results. Confidence intervals 

highlighted in gray are significant at that percent....................................................................... 139

Table 9. Kincaid Model Comparison Measure (MCM) results. Confidence intervals highlighted 

in gray are significant at that percent.......................................................................................... 140

Table 10. Bowline Model Comparison Measure (MCM) results. Confidence intervals 

highlighted in gray are significant at that percent....................................................................... 141

Table 11. Baldwin Model Comparison Measure (MCM) results. Confidence intervals 

highlighted in gray are significant at that percent....................................................................... 142

Table 12. Clifty Creek Model Comparison Measure (MCM) results. Confidence intervals 

highlighted in gray are significant at that percent....................................................................... 143
1

1. Introduction

1.1 Background

In 1991, the American Meteorological Society (AMS) and the U.S. Environmental 

Protection Agency (EPA) initiated a formal collaboration with the designed goal of introducing 

current planetary boundary layer (PBL) concepts into regulatory dispersion models. A working 

group (AMS/EPA Regulatory Model Improvement Committee, AERMIC) comprised of AMS 

and EPA scientists was formed for this collaborative effort.

In most air quality applications one is concerned with dispersion in the PBL, the 

turbulent air layer next to the earth's surface that is controlled by the surface heating and friction 

and the overlying stratification. The PBL typically ranges from a few hundred meters in depth at 

night to 1 - 2 km during the day. Major developments in understanding the PBL began in the 

1970's through numerical modeling, field observations, and laboratory simulations; see 

Wyngaard (1988) for a summary. For the convective boundary layer (CBL), a milestone was 

Deardorff's (1972) numerical simulations which revealed the CBL's vertical structure and 

important turbulence scales. Major insights into dispersion followed from laboratory 

experiments, numerical simulations, and field observations (e.g., see Briggs (1988), Lamb 

(1982) and Weil (1988a) for reviews). For the stable boundary layer (SBL), advancements 

occurred more slowly. However, a sound theoretical/experimental framework for surface layer 

dispersion and approaches for elevated sources emerged by the mid 1980's (e.g., see Briggs 

(1988) and Venkatram (1988)).

During the mid-1980's, researchers began to apply this information to simple dispersion 

models for applications. This consisted of eddy-diffusion techniques for surface releases, 

statistical theory and PBL scaling for dispersion parameter estimation, a new probability density 

function (pdf) approach for the CBL, simple techniques for obtaining meteorological variables 

(e.g., surface heat flux) needed for turbulence parameterizations, etc. Much of this work was 

reviewed and promoted in workshops (Weil 1985), revised texts (Pasquill and Smith 1983), and 

in short courses and monographs (Nieuwstadt and van Dop 1982; Venkatram and Wyngaard 

1988). By the mid 1980's, new applied dispersion models based on this technology had been 
2

developed including PPSP (Weil and Brower 1984), OML (Berkowicz et al. 1986), HPDM 

(Hanna and Paine 1989), TUPOS (Turner et al. 1986), CTDMPLUS (Perry et al. 1989); later, 

ADMS developed in the United Kingdom (see Carruthers et al. (1992)) was added as well as 

SCIPUFF (Sykes et al. 1996). AERMIC members were involved in the development of three of 

these models - PPSP, CTDMPLUS and HPDM.

By the mid-to-late 1980's, a substantial scientific base on the PBL and new dispersion 

approaches existed for revamping regulatory dispersion models, but this did not occur. In a 

review of existing or proposed regulatory models developed prior to 1984, Smith (1984) 

reported that the techniques were many years behind the state-of-the-art and yielded predictions 

that did not agree well with observations. Similar findings were reported by Hayes and Moore 

(1986), who summarized 15 model evaluation studies. The need for a comprehensive overhaul 

of EPA's basic regulatory models was clearly recognized. This need, including a summary of 

background information and recommendations, was the focus of an AMS/EPA Workshop on 

Updating Applied Diffusion Models held 24-27 January 1984 in Clearwater, Florida (see Weil 

(1985) and other review papers in the November 1985 issue of the Journal of Climate and 

Applied Meteorology.

In February 1991, the U.S. EPA in conjunction with the AMS held a workshop for state 

and EPA regional meteorologists on the parameterization of PBL turbulence and state-of-the-art 

dispersion modeling. One of the outcomes of the workshop was the formation of AERMIC. As 

noted above, the expressed purpose of the AERMIC activity was to build upon the earlier model 

developments and to provide a state-of-the-art dispersion model for regulatory applications. 

The early efforts of the AERMIC group are described by Weil (1992). In going through the 

design process and in considering the nature of present regulatory models, AERMIC’s goal 

expanded from its early form. In addition to improved parameterization of PBL turbulence, 

other problems such as plume interaction with terrain, surface releases, building downwash and 

urban dispersion were recognized as needing attention.

The new model developed by AERMIC is aimed at short-range dispersion from 

stationary industrial sources, the same scenario handled by the EPA Industrial Source Complex 

Model, ISC3 (U.S. Environmental Protection Agency, 1995). This work clearly has benefitted 
3

from the model development activities of the 1980's especially in the parameterization of mean 

winds and PBL turbulence, dispersion in the CBL, and the treatment of plume/terrain 

interactions. Techniques used in the new model for PBL parameterizations and CBL dispersion 

are similar to those used in earlier models. Turbulence characterization in the CBL adopts 

"convective scaling" as suggested by Deardorff (1972) as is included in most of the models 

mentioned above (e.g., PPSP, OML, and HPDM). Algorithms used in these earlier models were 

considered along with variants and improvements to them. In addition, the developers of OML 

met with AERMIC to discuss their experiences. Thus, much credit for the AERMIC model 

development is to be given to the pioneering efforts of the 1980s.

1.2 The AERMIC focus: a replacement for the ISC3 model

AERMIC’s initial focus has been on the regulatory models that are designed for 

estimating near-field impacts from a variety of industrial source types. EPA’s regulatory 

platform for near-field modeling, during the past 25 years has, with few exceptions, remained 

fundamentally unchanged. During this period, ISC3 was the workhorse regulatory model (used 

in the construction of most State Implementation Plans, new source permits, risk assessments 

and exposure analysis for toxic air pollutants) with code structure that is conducive to change. 

Therefore, AERMIC selected the EPA’s ISC3 Model for a major overhaul. AERMIC’s 

objective was to develop a complete replacement for ISC3 by: 1) adopting ISC3's input/output 

computer architecture; 2) updating, where practical, antiquated ISC3 model algorithms with 

newly developed or current state-of-the-art modeling techniques; and 3) insuring that the source 

and atmospheric processes presently modeled by ISC3 will continue to be handled by the 

AERMIC Model (AERMOD), albeit in an improved manner. 

The AERMOD modeling system consists of two pre-processors and the dispersion 

model. The AERMIC meteorological preprocessor (AERMET) provides AERMOD with the 

meteorological information it needs to characterize the PBL. The AERMIC terrain preprocessor (AERMAP) both characterizes the terrain, and generates receptor grids for the 

dispersion model (AERMOD).
4

AERMET uses meteorological data and surface characteristics to calculate boundary 

layer parameters (e.g. mixing height, friction velocity, etc.) needed by AERMOD. This data, 

whether measured off-site or on-site, must be representative of the meteorology in the modeling 

domain. AERMAP uses gridded terrain data for the modeling area to calculate a representative 

terrain-influence height associated with each receptor location. The gridded data is supplied to 

AERMAP in the format of the Digital Elevation Model (DEM) data (USGS 1994). The terrain 

preprocessor can also be used to compute elevations for both discrete receptors and receptor 

grids.

In developing AERMOD, AERMIC adopted design criteria to yield a model with 

desirable regulatory attributes. It was felt that the model should: 1) provide reasonable 

concentration estimates under a wide variety of conditions with minimal discontinuities; 2) be 

user friendly and require reasonable input data and computer resources as is the case with the 

ISC3 model; 3) capture the essential physical processes while remaining fundamentally simple; 

and, 4) accommodate modifications with ease as the science evolves.

Relative to ISC3, AERMOD currently contains new or improved algorithms for: 1) 

dispersion in both the convective and stable boundary layers; 2) plume rise and buoyancy; 3) 

plume penetration into elevated inversions; 4) computation of vertical profiles of wind, 

turbulence, and temperature; 5) the urban nighttime boundary layer; 6) the treatment of 

receptors on all types of terrain from the surface up to and above the plume height; 7) the 

treatment of building wake effects; 8) an improved approach for characterizing the fundamental 

boundary layer parameters; and 9) the treatment of plume meander.

1.3 Model development process

A seven step model development process, followed by AERMIC, resulted in the 

promulgation of a regulatory replacement for the ISC3 model, AERMOD. The process 

followed is as follows: 1) initial model formulation; 2) developmental evaluation; 3) internal 

peer review and beta testing; 4) revised model formulation; 5) performance evaluation and 

sensitivity testing; 6) external peer review; and 7) submission to EPA’s Office of Air Quality 

Planning and Standards (OAQPS) for consideration as a regulatory model.
5

The initial formulations of AERMOD are summarized in Perry et al. (1994) and 

Cimorelli et al. (1996). Once formulated, the model was tested (developmental evaluation) 

against a variety of field measurements in order to identify areas needing improvement. The 

developmental evaluation provided a basis for selecting formulation options.

This developmental evaluation was conducted using five data bases. Three consisted of 

event-based tracer releases, while the other two each contain up to a full year of continuous SO2

measurements. These data bases cover elevated and surface releases, complex and simple 

terrain, and rural and urban boundary layers. A description of the early developmental 

evaluation is presented in Lee et al. (1995) and in a later report by Lee et al. (1998). 

Additionally, a comprehensive peer review (U.S. Environmental Protection Agency, 2002) was 

conducted. Many revisions to the original formulation have resulted from this evaluation and 

comments received during the peer review, beta testing, and the public forum at EPA’s Sixth 

Conference on Air Quality Modeling (in 1995). Lee et al. (1998) describe the developmental 

evaluation repeated with the current model (i.e., revisions based on the developmental 

evaluation and peer review).

In addition, AERMOD underwent a comprehensive performance evaluation (Brode 

2002) designed to assess how well AERMOD’s concentration estimates compare against a 

variety of independent data bases and to assess the adequacy of the model for use in regulatory 

decision making. That is, how well does the model predict concentrations at the high end of the 

concentration distribution? AERMOD was evaluated against five independent data bases (two 

in simple terrain and three in complex terrain), each containing one full year of continuous SO2

measurements. Additionally, AERMOD’s performance was compared against the performance 

of four other applied, regulatory models: ISC3 (U.S. Environmental Protection Agency, 1995), 

CTDMPLUS (Perry 1992), RTDM (Paine and Egan 1987), and HPDM (Hanna and Paine 1989; 

Hanna and Chang 1993). The performance of these models against AERMOD has been 

compared using the procedures in EPA’s “Protocol for Determining the Best Performing 

Model” (U.S. Environmental Protection Agency, 1992).
6

On 21 April 2000 EPA proposed1 that AERMOD be adopted as a replacement to ISC3 

in appendix A of the Guideline on Air Quality Models (Code of Federal Regulations 1997). As 

such, upon final action, AERMOD would become EPA’s preferred regulatory model for both 

simple and complex terrain. Furthermore, on 19 May 2000 EPA announced2 its intention to

hold the Seventh Conference on Air Quality Modeling on 28-29 June 2000. The purpose of this 

conference was to receive comments on the April 2000 proposal. At the Seventh Conference, 

results of the performance evaluation and peer review were presented and public comments 

were received. Based on these comments AERMOD was revised to incorporate the PRIME 

algorithms for building downwash, to remove the dependency on modeling domain in 

AERMOD’s complex terrain formulation, and a variety of other less significant issues. A 

description of the fully revised model is presented here and in Cimorelli et al. (2004) and Perry 

et al. (2003). Performance of the final version of AERMOD is documented in Perry et al. 

(2003) and Brode (2002).

1.4 Purpose of the document

The purpose of this document is to provide a comprehensive, detailed description of the 

technical formulation of AERMOD and its preprocessors. This document is intended to provide 

many of the details that are not included in the published journal articles (Cimorelli et al. 2004; 

Perry et al. 2003).

This document does not include information related to model performance. As 

mentioned above, a description of the performance of the model that is described in this 

document can be found in Perry et al. (2003) and Brode (2002). 

 

1

40 CFR Part 51 pages 21506-21546

2

Federal Register on May 19, 2000 (Volume 65, Number 98)
7

2. Model overview

This section provides a general overview of the most important features of AERMOD. 

With the exception of treating pollutant deposition, AERMOD serves as a complete replacement 

for ISC3. However, it is the intention of AERMIC to incorporate both dry and wet particle and 

gaseous deposition as well as source or plume depletion. Once this is accomplished this report 

will be revised to include a description of the deposition formulation. Thus, the AERMOD 

model described here is applicable to rural and urban areas, flat and complex terrain, surface and 

elevated releases, and multiple sources (including, point, area and volume sources). Every 

effort has been made to avoid model formulation discontinuities wherein large changes in 

calculated concentrations result from small changes in input parameters.

AERMOD is a steady-state plume model. In the stable boundary layer (SBL), it 

assumes the concentration distribution to be Gaussian in both the vertical and horizontal. In the 

convective boundary layer (CBL), the horizontal distribution is also assumed to be Gaussian, 

but the vertical distribution is described with a bi-Gaussian probability density function (pdf). 

This behavior of the concentration distributions in the CBL was demonstrated by Willis and 

Deardorff (1981) and Briggs (1993). Additionally, in the CBL, AERMOD treats “plume 

lofting,” whereby a portion of plume mass, released from a buoyant source, rises to and remains 

near the top of the boundary layer before becoming mixed into the CBL. AERMOD also tracks 

any plume mass that penetrates into the elevated stable layer, and then allows it to re-enter the 

boundary layer when and if appropriate. For sources in both the CBL and the SBL AERMOD 

treats the enhancement of lateral dispersion resulting from plume meander.

Using a relatively simple approach, AERMOD incorporates current concepts about flow 

and dispersion in complex terrain. Where appropriate the plume is modeled as either impacting 

and/or following the terrain. This approach has been designed to be physically realistic and 

simple to implement while avoiding the need to distinguish among simple, intermediate and 

complex terrain, as required by other regulatory models. As a result, AERMOD removes the 

need for defining complex terrain regimes. All terrain is handled in a consistent and continuous 

manner while considering the dividing streamline concept (Snyder et al. 1985) in stablystratified conditions.
8

One of the major improvements that AERMOD brings to applied dispersion modeling is 

its ability to characterize the PBL through both surface and mixed layer scaling. AERMOD 

constructs vertical profiles of required meteorological variables based on measurements and 

extrapolations of those measurements using similarity (scaling) relationships. Vertical profiles 

of wind speed, wind direction, turbulence, temperature, and temperature gradient are estimated 

using all available meteorological observations. AERMOD is designed to run with a minimum 

of observed meteorological parameters. As a replacement for the ISC3 model, AERMOD can 

operate using data of a type that is readily available from National Weather Service (NWS) 

stations. AERMOD requires only a single surface measurement of wind speed (measured 

between 7zo and 100m - where zo is the surface roughness height), wind direction and ambient 

temperature. Like ISC3, AERMOD also needs observed cloud cover. However, if cloud cover 

is not available (e.g. from an on-site monitoring program) two vertical measurements of 

temperature (typically at 2 and 10 meters), and a measurement of solar radiation can be 

substituted. A full morning upper air sounding (rawinsonde) is required in order to calculate the 

convective mixing height throughout the day. Surface characteristics (surface roughness, 

Bowen ratio, and albedo) are also needed in order to construct similarity profiles of the relevant 

PBL parameters. 

Unlike existing regulatory models, AERMOD accounts for the vertical inhomogeneity of 

the PBL in its dispersion calculations. This is accomplished by "averaging" the parameters of 

the actual PBL into "effective" parameters of an equivalent homogeneous PBL.

Figure 1 shows the flow and processing of information in AERMOD. The modeling 

system consists of one main program (AERMOD) and two pre-processors (AERMET and 

AERMAP). The major purpose of AERMET is to calculate boundary layer parameters for use 

by AERMOD. The meteorological INTERFACE, internal to AERMOD, uses these parameters 

to generate profiles of the needed meteorological variables. In addition, AERMET passes all 

meteorological observations to AERMOD. 
9

Figure 1. Data flow in the AERMOD modeling system

Surface characteristics in the form of albedo, surface roughness and Bowen ratio, plus 

standard meteorological observations (wind speed, wind direction, temperature, and cloud 

cover), are input to AERMET. AERMET then calculates the PBL parameters: friction velocity 

(u*), Monin-Obukhov length (L), convective velocity scale (w*), temperature scale (θ*), mixing 

height (zi), and surface heat flux (H). These parameters are then passed to the INTERFACE 

(which is within AERMOD) where similarity expressions (in conjunction with measurements) 

are used to calculate vertical profiles of wind speed (u), lateral and vertical turbulent fluctuations 

(σv, σw), potential temperature gradient (dθ/dz), and potential temperature (θ ).

The AERMIC terrain pre-processor AERMAP uses gridded terrain data to calculate a 

representative terrain-influence height (hc), also referred to as the terrain height scale. The 

terrain height scale hc, which is uniquely defined for each receptor location, is used to calculate 

the dividing streamline height. The gridded data needed by AERMAP is selected from Digital 
10

Elevation Model (DEM) data. AERMAP is also used to create receptor grids. The elevation for 

each specified receptor is automatically assigned through AERMAP. For each receptor, 

AERMAP passes the following information to AERMOD: the receptor’s location (xr , yr), its 

height above mean sea level (zr), and the receptor specific terrain height scale (hc).

A comprehensive description of the basic formulation of the AERMOD dispersion 

model including the INTERFACE, AERMET, and AERMAP is presented in this document. 

Included are: 1) a complete description of the AERMET algorithms that provide quantitative 

hourly PBL parameters; 2) the general form of the concentration equation with adjustments for 

terrain; 3) plume rise and dispersion algorithms appropriate for both the convective and stable 

boundary layers; 4) handling of boundary layer inhomogeneity; 5) algorithms for developing 

vertical profiles of the necessary meteorological parameters; 6) a treatment of the nighttime 

urban boundary layer; 7) treatment of building downwash (incorporation of PRIME); and 8) 

enhancement of lateral dispersion due to plume meander. The model described here represents 

the 02222 versions of AERMOD, AERMET and AERMAP. In addition, all of the symbols 

used for the many parameters and variables that are referred to in this document are defined, 

with their appropriate units, in the section titled “List of Symbols.”
11

3. Meteorological preprocessor (AERMET)

The basic purpose of AERMET is to use meteorological measurements, representative of 

the modeling domain, to compute certain boundary layer parameters used to estimate profiles of 

wind, turbulence and temperature. These profiles are estimated by the AERMOD interface 

which is described in Section 0.

While the structure of AERMET is based upon an existing regulatory model 

preprocessor, the Meteorological Processor for Regulatory Models (MPRM) (Irwin et al. 1988),

the actual processing of the meteorological data is similar to that done for the CTDMPLUS 

(Perry 1992) and HPDM (Hanna and Paine 1989; Hanna and Chang 1993) models. 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 surface 

roughness, reflectivity (albedo), and the availability of surface moisture. The surface 

parameters provided by AERMET are the Monin-Obukhov Length (L), surface friction velocity 

(u*), surface roughness length (zo), surface heat flux (H), and the convective scaling velocity 

(w*). AERMET also provides estimates of the convective and mechanical mixed layer heights, 

zic and zim, respectively. AERMET defines the stability of the PBL by the sign of H (convective 

for H > 0 and stable for H < 0). Although AERMOD is capable of estimating meteorological 

profiles with data from as little as one measurement height, it will use as much data as the user 

can provide for defining the vertical structure of the boundary layer. In addition to PBL 

parameters, AERMET passes all measurements of wind, temperature, and turbulence in a form 

AERMOD needs.

3.1 Energy balance in the PBL

The fluxes of heat and momentum drive the growth and structure of the PBL. To 

properly characterize the PBL, one first needs a good estimate of the surface sensible heat flux 

(H) which depends on the net radiation (Rn) and surface characteristics such as the available 

surface moisture (described in the form of the Bowen ratio (Bo)). In the CBL, a simple energy 

balance approach, as in Oke (1978), is used to derive the expression, used in AERMET, to 
12

calculate the sensible heat flux, H. We begin with the following simple characterization of the 

energy balance in the PBL:

HH + λλ EE + GG = RRnn (1)

where H is the sensible heat flux, λE is the latent heat flux, G is the soil heat flux, and Rn is the 

net radiation. To arrive at an estimate of H simple parameterizations are made for the soil and 

latent heat flux terms; that is G=0.1 Rn and λ E = H / Bo, respectively. Substituting these 

expressions into eq. (1) the expression for surface heat flux becomes

HH = 0.9 RRnn

(1 + 1

BBoo � )

. (2)

3.1.1 Net radiation

If measured values for Rn are not available, the net radiation is estimated from the 

insolation and the thermal radiation balance at the ground following the method of Holtslag and 

van Ulden (1983) as

RRnn = (1 − rr{φφ}) RR + cc1 TTrrrrrr

6 − σσSSSS TTrrrrrr

4 − cc2 nn

1 + cc3

, (3)

where c1 = 5.31 x 10-13 W m-2 K-4, c2 = 60 W m-2, c3 = 0.12, σSB is the Stefan Boltzman constant 

(5.67 x 10-8 W m-2 K-4), Tref is the ambient air temperature at the reference height for temperature 

and Rn is the net radiation. The albedo is calculated as 

{φφ} = rr′ + (1 − rr′)exp(aa + bb),

where a = -0.1, b = -0.5(1-r'), and r' = r{φ-90°). Note, braces, {}, are used throughout this report 

to denote the functional form of variables.

Solar radiation, R, corrected for cloud cover, is taken from Kasten and Czeplak (1980) as
13

RR = RRoo (1 − 0.75nn3.4) , (4)

where n is the fractional cloud cover and Ro is the clear sky insolation which is calculated as 

Ro = 990(sin φ) - 30, and φφ(= (φφ�ttpp� + φφ{tt})⁄2) is the solar elevation (tp and t are the previous 

and present hours, respectively(1975). Note that when observations of cloud cover are 

unavailable a value of 0.5 is assumed in eq. (3) and measurements of solar radiation are required.

3.1.2 Transition between the CBL and SBL

When the PBL transitions from convective to stable conditions the heat flux changes 

sign from a positive to a negative value. At the point of transition the heat flux must therefore 

vanish, implying that the net radiation is equal to zero. By setting Ro equal to zero in eq. (3), 

and solving for sin φ, the critical solar elevation angle, φcrit, corresponding to the transition point 

between the CBL and the SBL can be determined from

sin(φφcc cccc) = 1

990 � −cc1TT6 + σσSSSSTT4 − cc2nn

(1 − rr{φφ})(1 − 0.75nn3.4) + 30�. (5)

Therefore, AERMET defines the point of transition between the CBL and SBL (day to 

night) as the point in time when the solar elevation angle φ = φcrit. On average, for clear and 

partly cloudy conditions, the transition from stable to convective conditions occurs when φ

reaches approximately 13°

; for overcast conditions φcrit increases to about 23° (Holtslag and van 

Ulden 1983).

However, if solar radiation measurements are available AERMET determines φcrit from 

an estimate of cloud cover rather than the actual observations themselves. In eq. (5) the cloud 

cover (n) is replaced with an equivalent cloud cover (neq) that is calculated from eq. (4) such that 

nn ee + �1 − RR⁄RRRR

0.75 �.
14

3.2 Derived parameters in the CBL

3.2.1 Friction velocity (u*) & Monin Obukhov length (L) in the CBL

In the CBL, AERMET computes the surface friction velocity, u*, and the MoninObukhov length, L, using the value of H estimated from eq. (2). Since the friction velocity and

the Monin Obukhov length depend on each other, an iterative method, similar to that used in 

CTDMPLUS (Perry 1992), is used. AERMOD initializes u*, and L by assuming neutral 

conditions (i.e., L=∞). The final estimate of u* and L is made once convergence is reached 

through iterative calculations (i.e., there is less than a 1% change between successive iterations). 

The expression for u* (e.g., Panofsky and Dutton (1984)) is

uu∗ = kk uurrrrrr

ln(zzrrrrrr⁄zz0) − ΨΨmm{zzrrrrrr⁄LL} + ΨΨmm{zz0⁄LL}

, (6)

where k is the von Karman constant (= 0.4), uref is the wind speed at reference height, zref is the

reference measurement height for wind in the surface layer, and zo is the roughness length. The 

stability terms (Ψm’s) in eq. (6) are computed as follows:

ΨΨmm �

zzrrrrrr

LL � = 2 ln �

1 + μ

2 � + ln �

1 + μ2

2 � − 2 tan−1 μ + ππ/2

ΨΨmm �

zz0

LL � = 2 ln �

1 + μ0

2 � + ln �

1 + μ0

2

2 � − 2 tan−1 μ0 + ππ/2

(7)

where μ = (1 - 16zref / L)1/4 and μ0 = (1 - 16z0 / L)1/4.

The initial step in the iteration is to solve eq. (6) for u* assuming that ψm = 0 (neutral 

limit) and setting u = uref. Having an initial estimate of u*, L is calculated from the following 

definition (e.g., see Wyngaard (1988)):
15

LL = ρρ ccpp TTrrrrrruu∗

3

kk gg HH (8)

where g is the acceleration of gravity, cp is the specific heat of air at constant pressure, ρ is the 

density of air, and Tref is the ambient temperature representative of the surface layer. Then u*

and L are iteratively recalculated using eqs. (6), (7) and (8) until the value of L changes by less 

than 1%.

The reference heights for wind speed and temperature that are used in determining the 

friction velocity and Monin-Obukhov length are optimally chosen to be representative of the 

surface layer in which the similarity theory has been formulated and tested with experimental 

data. Typically, a 10 m height for winds and a temperature within the range of 2 to 10 m is 

chosen. However, for excessively rough sites (such as urban areas with zo can be in excess of 

1 m), AERMET has a safeguard to accept wind speed reference data that range vertically 

between 7zo and 100 m. Below 7 zo (roughly, the height of obstacles or vegetation), 

measurements are unlikely to be representative of the general area. A similar restriction for 

temperature measurements is imposed, except that temperature measurements as low as zo are 

permitted. Above 100 m, the wind and temperature measurements are likely to be above the 

surface layer, especially during stable conditions. Therefore, AERMET imposes an upper limit 

of 100 meters for reference wind speed and temperature measurements for the purpose of 

computing the similarity theory friction velocity and Monin-Obukhov length each hour. Of 

course, other US EPA guidance for acceptable meteorological siting should be consulted in 

addition to keeping the AERMET restrictions in mind.

3.2.2 Convective velocity scale (w*)

AERMOD utilizes the convective velocity scale to characterize the convective portion of 

the turbulence in the CBL. Field observations, laboratory experiments, and numerical modeling 

studies show that the large turbulent eddies in the CBL have velocities proportional to the 

convective velocity scale (w*) (Wyngaard 1988). Thus in order to estimate turbulence in the 

CBL, an estimate of w* is needed. AERMET calculates the convective velocity scale from its 

definition as: 
16

ww∗ = � gg HH zz

ρρ ccpp TTrrrrrr�

1⁄3

, (9)

where zic is the convective mixing height (see Section 0).

3.3 Derived parameters in the SBL

In this section the parameters used to characterize the SBL are discussed along with their 

estimation methods. During stable conditions the energy budget term associated with the 

ground heating component is highly site-specific. During the day, this component is only about 

10% of the total net radiation, while at night its value is comparable to that of the net radiation 

(Oke 1978). Therefore, errors in the ground heating term can generally be tolerated during the 

daytime, but not at night. To avoid using a nocturnal energy balance approach that relies upon 

an accurate estimate of ground heating, AERMIC has adopted a much simpler semi-empirical 

approach for computing u* and L.

3.3.1 Friction velocity (u*) in the SBL

The computation of u* depends on the empirical observation that the temperature scale, 

θ* defined as

θθ∗ = −HH⁄ρ cp uu∗ , (10)

varies little during the night. Following the logic of Venkatram (1980) we combine the 

definition of L eq. (8) with eq. (10) to express the Monin-Obukhov length in the SBL as

LL = TTrrrrrr

k g θ∗

uu∗

2

.

(11)

From (Panofsky and Dutton 1984) the wind speed profile in stable conditions takes the 

form
17

uu = uu∗

kk �ln �

zz

zz0

� + ββmmzzrrrrrr

LL �, (12)

where βm = 5 and zref is the wind speed reference measurement height. Substituting eq. (11) into 

eq. (12) and defining the drag coefficient, CD, as k / ln(zref /z0) (Garratt 1992), results in

uu

uu∗

= 1

CCDD

+ ββ zzrrrrrrgg θθ∗

TTrrrrrr uu∗

2 . (13)

Multiplying eq. (13) by CD u*

2 and rearranging yields a quadratic of the form

uu∗

2 − CCDDuu uu∗ + CCDDuu0

2 = 0 (14)

where u0

2 = βm zref g θ* / Tref. As is used in HPDM (Hanna and Chang 1993) and CTDMPLUS 

(Perry 1992) this quadratic has a solution of the form

uu∗ = CCDD uurrrrrr

2

⎣

⎢

⎢

⎢

⎡

1 +

⎝

⎛1 − � 2 uu0

CCDD

1

2 uurrrrrr�

2

⎠

⎞

1/2

⎦

⎥

⎥

⎥

⎤

(15)

Equation (15) produces real-valued solutions only when the wind speed is greater than or 

equal to the critical value uucc = �4 ββmm zzrrrrrrgg θθ∗⁄TTrrrrrrCCDD�

1/2

. For the wind speed less than the 

critical value, u* and θ* are parameterized using the following linear expression:

uu∗ = uu∗{uu = uucc } � uu

uucc

� ff uu < uucc

θθ∗ = θθ∗ � uu

uucc

� ff uu < uucc

These expressions approximate the u* verses θ* dependence found by van Ulden and 

Holtslag (1983).
18

In order to calculate u* from eq. (15) an estimate of θ*is needed. If representative cloud

cover observations are available, the temperature scale in the SBL is taken from the empirical

form of van Ulden and Holtslag (1985) as

θθ∗ = 0.09 (1 − 0.5 nn2) , (16)

where n is the fractional cloud cover. However, if cloud cover measurements are not available, 

AERMET can estimate θ* from measurements of temperature at two levels and wind speed at 

one level. This technique, known as the Bulk Richardson approach, starts with the similarity 

expression for potential temperature (Panofsky and Dutton 1984), that is,

θθ{zz} − θθ0 = θθ∗

kk �ll

zz

zz0

+ ββmm

zz

LL

� , (17)

where βm ≈ 5 and k (= 0.4) is the von Karman constant. Applying eq. (17) to the two levels of

temperature measurements and rearranging terms yields

θθ∗ = kk(θθ2 − θθ1)

��ll zz2

zz1

� + ββmm

(zz2 − zz1)

LL �

. (18)

Since both u* (eq.(12)) and θ* (eq. (18)) depend on L, and L (eq. (11)) in turn depends on 

u* and θ*, an iterative approach is needed to estimate u*. First u* and θ* are found by assuming 

an initial value for L and iterating among the expressions for u*, θ* (eq. (18)) and L (eq. (11)) 

until convergence is reached. The expression used for u*, in the iteration, is taken from (Holtslag 

1984) and depends on atmospheric stability. For situations in which z/L < 0.5 is estimated using 

eq. (12), otherwise (for more stable cases) u* is calculated as follows:

( ) ( )

u

ku

z

z

z

o L z L z L

m

*

ln ln . . .

= .

+ + − + − 















7 4 25 0 5

2

1648 2

β (19)
19

3.3.2 Sensible heat flux (H) in the SBL

Having computed u* and θ* for stable conditions, AERMET calculates the surface heat 

flux from eq. (10) as

HH = ρρ ccppuu∗θθ∗ . (20)

AERMET limits the amount of heat that can be lost by the underlying surface to 64 W 

m-2

. This value is based on a restriction that Hanna (1986) placed on the product of θ* and u*. 

That is, for typical conditions Hanna found that

[θθ∗uu∗]mmmmmm = 0.05mmss−1 KK . (21)

When the heat flux, calculated from eq. (20), is such that θ* u* > 0.05 m s-1 K, AERMET 

recalculates u* by substituting 0.05/u* into eq. (15) for θ* (uo in eq. (15) is a function of θ*).

3.3.3 Monin-Obukhov length (L) in the SBL

Using the sensible heat flux of eq. (20) and u* from eq. (15), the Monin-Obukhov 

Length, for the SBL is calculated from eq. (8).

3.4 Mixing height 

The mixing height (zi) in the CBL depends on both mechanical and convective processes 

and is assumed to be the larger of a mechanical mixing height (zim) and a convective mixing 

height (zic). Whereas, in the SBL, the mixing height results exclusively from mechanical (or 

shear induced) turbulence and therefore is identically equal to zim. The same expression for 

calculating zim is used in both the CBL and the SBL. The following two sections describe the 

procedures used to estimate zic and zim, respectively.
20

3.4.1 Convective mixing height (zic)

The height of the CBL is needed to estimate the profiles of important PBL variables and 

to calculate pollutant concentrations. If measurements of the convective boundary layer height 

are available they are selected and used by the model. If measurements are not available, zic is 

calculated with a simple one-dimensional energy balance model (Carson 1973) as modified by 

Weil and Brower (1983). This model uses the early morning potential temperature sounding 

(prior to sunrise), and the time varying surface heat flux to calculate the time evolution of the 

convective boundary layer as

zz θθ{zz } − � θθ{zz}dd

zziiii

0

= (1 + 2AA) � HH{tt′

}

ρρ ccpp

dd

0

,

(22)

where θ is the potential temperature, A is set equal to 0.2 from Deardorff (1980), and t is the hour 

after sunrise. Weil and Brower found good agreement between predictions and observations of 

zic, using this approach.

3.4.2 Mechanical mixing height (zim)

In the early morning when the convective mixed layer is small, the full depth of the PBL 

may be controlled by mechanical turbulence. AERMET estimates the heights of the PBL during 

convective conditions as the maximum of the estimated (or measured if available) convective 

boundary layer height (zic) and the estimated (or measured) mechanical mixing height. 

AERMET uses this procedure to insure that in the early morning, when zic is very small but 

considerable mechanical mixing may exist, the height of the PBL is not underestimated. When 

measurements of the mechanical mixed layer are not available, zim is calculated by assuming 

that it approaches the equilibrium height given by Zilitinkevich (1972) as

zz = 0.4 �

uu∗LL

ff� �, (23)

where zie is the equilibrium mechanical mixing height and f is the Coriolis parameter.
21

Venkatram (1980) has shown that, in mid-latitudes, eq. (23) can be empirically 

represented as 

zz = 2300 uu∗

3/2

, (24)

where zie (calculated from eq. (24)) is the unsmoothed mechanical mixed layer height. When 

measurements of the mechanical mixed layer height are available they are used in lieu of zie.

To avoid estimating sudden and unrealistic drops in the depth of the shear-induced, 

turbulent layer, the time evolution of the mechanical mixed layer height (whether measured or 

estimated) is computed by relaxing the solution toward the equilibrium value appropriate for the 

current hour. Following the approach of Venkatram (1982)

ddzz

dd = (zz − zz )ττ . (25)

The time scale, τ , governs the rate of change in height of the layer and is taken to be 

proportional to the ratio of the turbulent mixed layer depth and the surface friction velocity (i.e. 

τ = zim / βτ u*). AERMOD uses a constant βτ value of 2. For example, if u* is of order 0.2 m s-1, 

and zim is of order 500 m, the time scale is of the order of 1250 s which is related to the time it 

takes for the mechanical mixed layer height to approach its equilibrium value. Notice that when 

zim < zie, the mechanical mixed layer height increases to approach its current equilibrium value; 

conversely, when zim > zie, the mechanical mixed layer height decreases towards its equilibrium 

value.

Because the friction velocity changes with time, the current smoothed value of zim{t+Δt} 

is obtained by numerically integrating eq. (25) such that

zz {tt + ΔΔ } = zz {tt}�ee−ΔΔ ⁄ττ� + zz {tt + ΔΔ }�1 − �ee−ΔΔ ⁄ττ�� . (26)

where zim{t} is the previous hour’s smoothed value. For computing the time scale in eq. (26), zim

is taken from the previous hour’s estimate and u* from the current hour. In this way, the time 
22

scale (and thus relaxation time) will be short if the equilibrium mixing height grows rapidly but 

will be long if it decreases rapidly.

Although eqs. (24) and (26) are designed for application in the SBL, they are used in the 

CBL to ensure a proper estimate of the PBL height during the short transitional period at the 

beginning of the day when mechanical turbulence generally dominates. The procedure, used by 

AERMET, guarantees the use of the convective mixing height once adequate convection has 

been established even though the mechanical mixing height is calculated during all convective 

conditions. Since AERMET uses eq. (26) to estimate the height of the mixed layer in the SBL, 

discontinuities in zi from night to day are avoided.

In AERMOD, the mixing height zi, has an expanded role in comparison to how it is used 

in ISC3. In AERMOD the mixing height is used as an elevated reflecting/penetrating surface, 

an important scaling height, and enters in the w* determination found in eq. (9). The mixing 

height zi for the convective and stable boundary layers is therefore defined as follows:

z MAX[z z ] for L CBL

z z for L SBL

i ic im

i im

= <

= <

; ( )

( )

0

0

(27)

Since algorithms used for profiling differ in the SBL and CBL, the stability of the PBL 

must be determined. For this purpose the sign of L is used by AERMET; if L < 0 then the PBL 

is considered to be convective (CBL) otherwise it is stable (SBL).

3.5 Adjustment for the low wind speed/stable conditions in AERMET

An option has been incorporated in AERMET to address issues associated with model 

overpredictions under low-wind/stable conditions. The ADJ_U* option is available in 

AERMET by specifying ADJ_U* on the METHOD STABLEBL keyword in the Stage 3

AERMET input file.
23

The ADJ_U* option can be specified with or without the Bulk Richardson Number 

option in AERMET, which utilizes delta-T measurements. The Bulk Richardson Number 

option in AERMET is selected by specifying BULKRN on a separate METHOD STABLEBL 

keyword in the Stage 3 AERMET input file. The formulation for the ADJ_U* option without 

the BULKRN option is based on Equation 26 of Qian and Venkatram (2011). The formulation 

for the ADJ_U* option with the BULKRN option is based on Equations 22, 23, and 25 of Luhar 

and Rayner (2009), with a critical value of z/L of 0.7.
24

4. Vertical structure of the PBL - AERMOD’S meteorological interface 

The AERMOD interface, a set of routines within AERMOD, uses similarity 

relationships with the boundary layer parameters, the measured meteorological data, and other 

site-specific information provided by AERMET to compute vertical profiles of: 1) wind 

direction, 2) wind speed, 3) temperature, 4) vertical potential temperature gradient, 5) vertical 

turbulence (σw ) and 6) lateral turbulence (σv ). 

For any one of these six variables (or parameters), the interface (in constructing the 

profile) compares each height at which a meteorological variable must be calculated with the 

heights at which observations were made and if it is below the lowest measurement or above the 

highest measurement (or in some cases data is available at only one height), the interface 

computes an appropriate value from selected PBL similarity profiling relationships. If data are 

available both above and below a given height, an interpolation is performed which is based on 

both the measured data and the shape of the computed profile (see Section 0). Thus the 

approach used for profiling, simultaneously takes advantage of the information contained in 

both the measurements and similarity parameterizations. As will be discussed, at least one level 

of measured wind speed, wind direction, and temperature is required. However, turbulence 

profiles can be parameterized without any direct turbulence measurements.

The following sections provide a comprehensive description of AERMOD’s profiling 

equations and how these estimated profiles are used to extract pertinent layer-averaged 

meteorology for AERMOD’s transport and dispersion calculations. Also, example profiles (one 

typical of the CBL and one typical of the SBL) for the various parameters have been constructed 

for illustration. The CBL case assumes that zi = 1000 m, L = -10 m and zo = 0.1 m (i.e., 

zo = 0.0001zi and L = - 0.01zi). The SBL case assumes that zi = 100 m, L = 10 m and zo = 0.1 

m (i.e., zo = 0.001 zi and L = 0.1 zi). 
25

4.1 General profiling equations

4.1.1 Wind speed profiling

The AERMOD profile equation for wind speed, has the familiar logarithmic form:

(28)

At least one wind speed measurement, that is representative of the surface layer, is 

required for each simulation with AERMOD. Since the logarithmic form does not adequately 

describe the profile below the height of obstacles or vegetation, eq. (28) allows for a linear 

decrease in wind speed from its value at 7zo.

For the CBL, the Ψm’s are evaluated using eq. (7) with zref replaced by z, and during 

stable conditions they are calculated from van Ulden & Holtslag (1985) as

(29)

For small z/L (<<1) and with a series expansion of the exponential term, the first 

equation in (29) reduces to the form given in eq. (12), i.e., ψm = βm z/L with βm = 5. However, 

for large z/L (>1) and heights as great as 200 m in the SBL, the ψm given by eq. (29) is found to 

fit wind observations much better than the ψm given by eq. (12) (van Ulden and Holtslag 1985). 

Using the example case parameter values Figure 2 and Figure 3 were constructed to illustrate 

the form of the wind profiles used by AERMOD in the layers above and below 7zo.

{ }

{ }

u u z

z

z

for z z

u

u

k

z

z

z

L

z

L for z z z

u u z for z z

o

o

o

o

m m

o

o i

i i

= 



 



 <

= 



 



 − 









 +

















 



 ≤ ≤

= >

7

7

7

7 * ln Ψ Ψ

Ψ

Ψ

m

m

o o

z

L

z

L

z

L

z

L











 = − − −



 



 



 















 = − − −





 



 



 





17 1 0 29

17 1 0 29

exp .

exp . .
26

Figure 2. Wind speed profile, for both the CBL and SBL, in the region below 7z0
27

4.1.2 Wind direction profiles

For both the CBL & SBL wind direction is assumed to be constant with height both 

above the highest and below the lowest measurements. For intermediate heights, AERMOD 

linearly interpolates between measurements. At least one wind direction measurement is 

required for each AERMOD simulation.

4.1.3 Profiles of the potential temperature gradient

Above the relatively shallow superadiabatic surface layer, the potential temperature 

gradient in the well mixed CBL is taken to be zero. The gradient in the stable interfacial layer 

just above the mixed layer is taken from the morning temperature sounding. This gradient is an 

important factor in determining the potential for buoyant plume penetration into and above that 

Figure 3. Wind speed profiling, for both the CBL and SBL, in the region above 

7z0
28

layer. Above the interfacial layer, the gradient is typically constant and slightly stable. 

Although the interfacial layer depth varies with time, for the purposes of determining the 

strength of the stable stratification aloft, AERMET uses a fixed layer of 500 m to insure that a 

sufficient layer of the morning sounding is sampled. A 500 m layer is also used by the 

CTDMPLUS model (Perry 1992) for this same calculation. This avoids strong gradients 

(unrealistic kinks) often present in these data. For a typical mixed layer depth of 1000 m an 

interfacial layer depth of 500 m is consistent with that indicated by Deardorff (1979). A 

constant value of 0.005 K m-1 above the interfacial layer is used as suggested by Hanna and 

Chang (1991). Using the morning sounding to compute the interfacial temperature gradient 

assumes that as the mixed layer grows throughout the day, the temperature profile in the layer 

above zi changes little from that of the morning sounding. Of course, this assumes that there is 

neither significant subsidence nor cold or warm air advection occurring in that layer. Field 

measurements (e.g. Clarke et al. (1971)) of observed profiles throughout the day lend support to 

this approach. These data point out the relative invariance of upper level temperature profiles 

even during periods of intense surface heating.

Below 100 m, in the SBL, AERMOD uses the definition of the potential temperature 

gradient suggested by Dyer (1974) as well as Panofsky and Dutton (1984). That is,

(30)

Eq. (30) is similar to that of Businger et al. (1971). Above 100 m the form of the 

potential temperature gradient, taken from Stull (1983) and van Ulden & Holtslag (1985) is 

(31)

( )

∂θ ( )

∂

θ

∂θ

∂

θ

z k L for z m

z k z

z

L for m z m

= +





 



 ≤

= +





 



 < ≤

*

* .

2

1 5

2

2

1 5 2 100

∂θ { } ( )

∂

∂θ

∂ θ z

z

z

z z

z

mx mx

i

= −

 −



 



 exp 0.44
29

where zmx = 100 m, ziθ = max[zim; 100m], and the constant 0.44 within the exponential term of eq. 

(31) is inferred from typical profiles taken during the Wangara experiment (Andre and Mahrt 

1982). For all z, ∂θ/∂z is limited to a minimum of 0.002 K m-1 (Paine and Kendall 1993).

In the SBL if dθ/dz measurements are available below100 m and above zo, then θ* is 

calculated from eq. (30) using the value of ∂θ/∂z at the lowest measurement level and zTref

replaced by the height of the ∂θ/∂z measurements. The upper limit of 100 m for the vertical 

temperature gradient measurements is consistent with that imposed by AERMET for wind speed 

and temperature reference data used to determine similarity theory parameters such as the 

friction velocity and the Monin-Obukhov length. Similarly, the lower limit of zo for the vertical 

temperature gradient measurements is consistent with that imposed for reference temperature 

data. If no measurements of ∂θ/∂z are available, in that height range, then θ* is calculated by 

combining eqs. (8) and (20). θ* is not used in the CBL.

Figure 4 shows the inverse height dependency of ∂θ/∂z in the SBL. To create this curve 

we assumed that: Zim=100 m; and therefore, Zi θ = 100 m; L = 10 m; u*= .124, which is 

consistent with a mixing height of 100 m; Tref = 293 K; and therefore based on eq. (11) θ* =

0.115 K. These parameter values were chosen to represent a strongly stable boundary layer. 

Below 2 m ∂θ/∂z is persisted downward from its value of 0.228 K m-1 at 2m. Above 100 m 

∂θ/∂z is allowed to decay exponentially with height.
30

Figure 4. Profile of potential temperature gradient for the SBL

4.1.4 Potential temperature profiling

For use in plume rise calculations, AERMOD develops the vertical profile of potential 

temperature from its estimate of the temperature gradient profile. First, the model computes the 

potential temperature at the reference height for temperature (i.e., zTref) as

(32)

where zmsl = zref + zbase is the user specified elevation for the base of the temperature profile (i.e., 

meteorological tower). Then for both the CBL and SBL the potential temperature is calculated 

as follows:

θ {z } T , g z

c Tref ref

msl

P

= +
31

(33)

where ∂∂

∂∂

� is the average potential temperature gradient over the layer Δz. Note that for 

z < zTref , Δz is negative.

4.1.5 Vertical turbulence calculated

In the CBL, the vertical velocity variance or turbulence (σσwwww

2 ) is profiled using an 

expression based on a mechanical or neutral stability limit (σwm ∝ u*) and a strongly convective 

limit (σwc ∝ w*). The total vertical turbulence is given as:

(34)

This form is similar to one introduced by Panofsky et al. (1977) and included in other 

dispersion models (e.g., Berkowicz et al. (1986), Hanna and Paine (1989), and Weil (1988a)).

The convective portion (σ 2

wc ) of the total variance is calculated as:

(35)

where the expression for z ≤ 0.1 zic is the free convection limit (Panofsky et al. 1977), for 

0.1zi < z ≤ zic is the mixed-layer value (Hicks 1985), and for z > zic is a parameterization to 

connect the mixed layer σ 2

wc to the assumed near-zero value well above the CBL. An example 

profile of convective vertical turbulence described in eq. (35) is presented in Figure 5. 

θ θ ∂θ

∂ {z z} {z}

z

z

z

+ ∆ = + ∆

θ θ ∂θ

∂ {z z} {z}

z

z

z

+ ∆ = + ∆

( )

σ

σ

σ

wc

ic

ic

wc ic ic

wc

ic

ic

ic

z

z

w for z z

w for z z z

w

z z

z

for z z

2

2 3

2

2 2

2 2

16 01

0 35 01

0 35

6

= 



 



 ⋅ ≤

= < ≤

= −

 −



 



 >

. .

. .

. exp

*

*

*
32

Figure 5. Convective portion of the vertical turbulence in the CBL

The mechanical turbulence (σwm) is assumed to consist of a contribution from the 

boundary layer (σwml) and from a “residual layer” (σwmr) above the boundary layer (z > zi) such 

that,

(36)

This is done to satisfy the assumed decoupling between the turbulence aloft (z > zi) and 

that at the surface in the CBL shear layer, and to maintain a continuous variation of σ 2

wm with z

near z = zi. The expression for σwml following the form of Brost et al. (1982) is

σ wm σ wml σ wmr

2 2 2 = + .
33

(37)

where the σwml =1.3u* at z = 0 is consistent with Panofsky et al. (1977).

Above the mixing height σwmr is set equal to the average of measured values in the 

residual layer above zi. If measurements are not available, then σwmr is taken as the default value 

of 0.02 u{zi}. The constant 0.02 is an assumed turbulence intensity iz ( = σwm / u) for the very 

stable conditions presumed to exist above zi (Briggs 1973). Within the mixed layer the residual 

turbulence (σwmr) is reduced linearly from its value at zi to zero at the surface. Figure 6 presents 

the profile of the mechanical portion of the vertical turbulence in the CBL. The effect of 

combining the residual and boundary layer mechanical turbulence (eq. (36)) can be seen in this

figure.

Figure 6. Mechanical portion of the vertical turbulence in the CBL

σ

σ

wml

i

i

wml i

u

z

z

for z z

for z z

= − 



 



 <

= ≥

13 1

0 0

1 2

.

.

*
34

In the SBL the vertical turbulence contains only a mechanical portion which is given by 

eq. (36). The use of the same σ 2

wm expressions for the SBL and CBL is done to ensure 

continuity of turbulence in the limit of neutral stability. Figure 7 illustrates AERMOD’s 

assumed vertical turbulence profile for the SBL. This is similar to the profile for the CBL 

except for a notable increase in the value of σwmr. Since values for σwmr are based on the 

magnitude of the wind speed at zi, the differences in the two figures stem from setting

zo = 0.0001zi in the CBL example case while for the SBL case zo = 0.001zi.

4.1.6 Lateral turbulence calculated by the interface

In the CBL the total lateral turbulence, σσvv

2 , is computed as a combination of a 

mechanical (σvm) and convective (σvc) portions such that

σ σ σ (38) vT vc vm

2 2 2 = + .

Figure 7. Profile of vertical turbulence in the SBL
35

In the SBL the total lateral turbulence contains only a mechanical portion. AERMOD, 

uses the same σvm expression in the CBL and SBL. This is done to maintain continuity of σvm in 

the limit of neutral stability. A description of mechanical and convective profiles of lateral 

turbulence follows.

4.1.6.1 Mechanical portion of the lateral turbulence

The variation with height of the mechanical portion of the lateral turbulence is bounded 

by its value at the 1surface and an assumed residual value at the top of the mechanical mixed 

layer. The variation between these two limits is assumed to be linear. Based on observations 

from numerous field studies, Panofsky and Dutton (1984) report that, in purely mechanical 

turbulence, the lateral variance near the surface has the form 

(39)

where the constant, C, ranges between 3 and 5. Based on an analysis of the Kansas data, 

Izumi (1971) and Hicks (1985) support the form of eq. (39) with a value of 3.6 for C. 

Between the surface and the top of the mechanically mixed layer, σ2

vm is assumed to 

vary linearly as

(40)

where σσvv

2 {zzmm} = min[σσvv2 ; 0.25 mm2ss−2] and σvo2, the surface value of the lateral turbulence, is 

equal to 3.6 u*

2

. This linear variation of σ 2

vm with z is consistent with field observations (e.g., 

Brost et al. (1982)). In the SBL the total lateral turbulence contains only a mechanical portion 

and it is given by eq. (40).

σ v Cu 0

2 2 = *

{ }

{ }

σ

σ σ

σ

σ σ

vm

vm im vo

im

vo im

vm vm im im

z

z

z for z z

z for z z

2

2 2

2

2 2

=

 −















+ ≤

= > ,
36

Above the mixed layer, lateral turbulence is expected to maintain a modest residual 

level. Hanna (1983) analyzed ambient measurements of lateral turbulence in stable conditions. 

He found that even in the lightest wind conditions, the measurements of σvc were typically 0.5 m 

s-1, but were observed to be as low as 0.2 m s-1. AERMOD adopts the lower limit of 0.2 m s-1

for σvc in near-surface conditions, as discussed below, but uses the more typical value of 0.5 m s1 for the residual lateral turbulence above the mixed layer. Above the height of the CBL, the 

model linearly decreases σvc2 from σvc2{ zic } to 0.25 at 1.2 zic and holds σvc2 constant above

1.2 zic. However, if σvc2{zic } < 0.25 m2 s-2, then σvc2{zic } is persisted upward from zic. 

Furthermore, it was found that a value of the order σvc2 = 0.25 m2 s-2 provided consistently good 

model performance (for plumes commonly above zim) during the developmental evaluation 

(Perry et al. 2005) supporting the presence of residual lateral turbulence in this layer.

Figure 8 shows how the vertical profile of lateral mechanical turbulence changes over a 

range of mechanical mixing heights, and related friction velocities. The values of u* used to 

produce these curves are consistent with the relationship between zim and u* which is found in 

eq. (24). For the SBL Figure 8 represents profiles of the total lateral turbulence. In the CBL 

these curves depict only the mechanical portion of the total lateral variance. Note that for zim = 

300 m and 100 m the values σvo2 are less than 0.25 m2 s-2. Therefore, the profiles are constant 

with height. 
37

Figure 8. Family of lateral mechanical turbulence profiles over a range of mechanical mixing 

heights

4.1.6.2 Convective portion of the lateral turbulence

The convective portion of the lateral turbulence within the mixed layer is constant and 

calculated as:

(41)

This constant value of σσww

2 ww∗2 ⁄ = 0.35 is supported by the Minnesota data (Readings et al. 

1974; Kaimal et al. 1976) and by data collected at Ashchurch England (Caughey and Palmer 

1979).

σ vc w

2 2 = 0.35 *
38

For z > zic , the model linearly decreases σvc2 from σvc2{ zic } to 0.25 at 1.2 zic and holds 

σvc2 constant above 1.2 zic. However, if σvc2{zic } < .25 m2 s-2, then σvc2{zic } is persisted upward 

from zic.

4.2 Vertical inhomogeneity in the boundary later as treated by the interface

AERMOD is designed to treat the effects on dispersion from vertical variations in wind 

and turbulence. Consideration of the vertical variation in meteorology is important for properly 

modeling releases in layers with strong gradients, for capturing the effects of meteorology in 

layers into which the plume may be vertically dispersing, and to provide a mechanism (in the 

CBL) by which sources that are released into or penetrate into an elevated stable layer can 

eventually re-enter the mixed layer. However, AERMOD is a steady-state plume model and 

therefore can use only a single value of each meteorological parameter to represent the layer 

through which these parameters are varying. Thus, the model "converts" the inhomogeneous 

values into equivalent effective or homogeneous values. This technique is applied to u, σvT, σwT, 

∂θ/∂z and the Lagrangian time scale. The effective parameters are denoted by a tilde throughout 

the document (e.g., effective wind speed is denoted by uu�).

Fundamental to this approach is the concept that the primary layer of importance, 

relative to receptor concentration, is the one through which plume material travels directly from 

source to receptor. Figure 9 presents a schematic illustration of the approach AERMOD uses to 

determine these effective parameters (α is used to generically represent these parameters). The 

effective parameters are determined by averaging their values over that portion of the layer that 

contains plume material between the plume centroid height, Hp {x}, (a simplified surrogate for 

the height of the plume’s center of mass) and the receptor height (zr). In other words, the 

averaging layer is determined by the vertical half-depth of the plume (defined as 2.15 σz {xr} 

where xr is the distance from source to receptor) but is bounded by Hp {xr} and zr. The values 
39

used in the averaging process are taken from AERMOD’s vertical profiles. This technique is 

best illustrated with examples.

Consider the two receptors depicted in Figure 9. Both receptors are located at the same 

distance xr from the source but at different heights above ground, i.e., zr1 and zr2. An example 

profile of some parameter α is shown at the far left of the figure. The value of the effective 

parameter used by AERMOD to represent transport and diffusion from source to receptor 

depends on the location of the receptor. For receptor 1 the effective parameter value αα�1 (shown 

in the figure as αeff1) is determined by averaging the values of α {z} between Hp {xr} and zr1. 

Therefore, the layer over which this average is taken is smaller than the plume’s half-depth. 

Whereas, αα�2 (shown in the figure as αeff2) is determined by averaging α {z} over the full layer 

from Hp {xr} down through a depth of 2.15σz {xr} since the receptor is located below the defined 

lower extent of the plume.

Figure 9. AERMOD’s Treatment of the Inhomogeneous Boundary Layer
40

Since σz {xr }depends on the effective values of σwT and u, the plume size is estimated by 

first using the plume height values of σwT {Hp } and u {Hp } to calculate σz {xr}. As illustrated 

in Figure 9 σz{xr} is then used to determine the layer over which σσ�wwww{xxrr} and uu�{xxrr} are 

calculated. Once the averaging layer for a given plume and receptor is established the effective 

values, αα�, are computed as simple averages:

(42)

where hb and ht are the bottom and top, respectively, of the layer of importance such that:

(43)

For all plumes, both limits are bounded by either the zr or Hp. For both the direct and 

indirect sources ht, in eq. (43) is not allowed to exceed zi and if hb ≥ zi then αα� = αα{zz }.

For plumes in stable conditions and for the penetrated source in the CBL, Hp is always 

set equal to the plume centerline height (Δhss + hss) where hs is the stack height corrected for 

stack tip downwash and Δhs is the stable source plume rise. The stable source plume rise Δhs is 

calculated from expressions found in Section 0.

In the CBL, the specification of Hp is somewhat more complicated. Because of limited 

mixing in the CBL the center of mass of the plume will be the plume height close to the source 

and the mid-point of the PBL at the distance where it becomes well mixed. Beyond final plume 

rise, Hp is varied linearly between these limits.

( ) { } ~

α = α − ∫ 1

h h

z dz

t b h

h

b

t

{ } { }

{[ { } { }] } { }

{[ { } { }] } { }

{ } { }

h

H x y if H x y z

MAX H x y x z if H x y z

h

MIN H x y x z if H x y z

H x y if H x y z

b

p r r p r r r

p r r z sr r p r r r

t

p r r z sr r p r r r

p r r p r r r

=

≤

− >





















= + ≤

>





















, , ,

, . , , ,

, . , , ,

, , ,

215

215

σ

σ
41

Prior to plume stabilization, i.e., x < xf (distance to plume stabilization),

HHpp = hss + ΔΔhdd,pp,

where Δhd is the plume rise for the direct source (estimated from eq. (91)), and Δhp (= hep - hs) is 

the plume rise for the penetrated source, where hep (penetrated source plume height) is calculated 

from eq. (94).

The distance to plume stabilization, xf , is determined following Briggs (Briggs 1975; 

Briggs 1971) as

(44)

where the buoyancy flux (Fb) is calculated from eq. (57). 

For Fb = 0 the distance to final rise is calculated from the ISCST3 ((U.S. Environmental 

Protection Agency, 1995)) expression

(45)

where up is the wind speed at source height, rs is the stack radius, and ws is the stack exit gas 

velocity.

Beyond plume stabilization ( x>xf ), Hp varies linearly between the stabilized plume 

height (H{xf}) and the mid-point of the mixed layer (zi /2). This interpolation is performed over 

the distance range xf to xm, where xm is the distance at which pollutants first become uniformity 

mixed throughout the boundary layer.

x F for F

x F for F

f b b

f b b

= <

= ≥

49 55

119 55

5 8

2 5

( ) x

r w u

w u f

s s p

s p

= 8 + 3

2
42

The distance xm is taken to be the product of the average mixed layer wind speed and the 

mixing time scale, zz /σσwwww. That is,

(46)

where the averaging of u and σwT are taken over the depth of the boundary layer.

For distances beyond xf , Hp is assumed to vary linearly between the plume's stabilized 

height, H {xf }, and zi /2 such that:

(47)

Note that in the CBL, both the direct and indirect source will have the same α (effective 

parameter) values. In eq. (43) σz is the average of the updraft σz and the downdraft σz , the 

maximum value of ht is zi , and when hb ≥ zi , α = α {zi}.

As discussed previously, when multiple vertical measurements of wind direction are 

available a profile is constructed by linearly interpolating between measurements and persisting 

the highest and lowest measurements up and down, respectively. The approach taken for 

selecting a transport wind direction from the profile is different from the above. The transport 

wind direction is selected as the mid-point of the range between stack height and the stabilized 

plume height.

x

u z m i

wT

= σ ,

( )

( ) H H x

z H x

x x

x x p f

i

f

f

m f

= + − 



 



 ⋅ −

− { } { } 2
43

5. The AMS/EPA regulatory Model: AERMOD

AERMOD is a steady-state plume model in that it assumes that concentrations at all 

distances during a modeled hour are governed by the temporally averaged meteorology of the 

hour. The steady state assumption yields useful results since the statistics of the concentration 

distribution are of primary concern rather than specific concentrations at particular times and 

locations. AERMOD has been designed to handle the computation of pollutant impacts in both 

flat and complex terrain within the same modeling framework. In fact, with the AERMOD 

structure, there is no need for the specification of terrain type (flat, simple, or complex) relative 

to stack height since receptors at all elevations are handled with the same general methodology. 

To define the form of the AERMOD concentration equations, it is necessary to simultaneously 

discuss the handling of terrain.

In the stable boundary layer (SBL), the concentration distribution is assumed to be 

Gaussian in both the vertical and horizontal. In the convective boundary layer (CBL), the 

horizontal distribution is assumed to be Gaussian, but the vertical distribution is described with 

a bi-Gaussian probability density function (pdf). This behavior of the concentration 

distributions in the CBL was demonstrated by Willis and Deardorff (1981) and Briggs (1993). 

Additionally, in the CBL, AERMOD treats “plume lofting,” whereby a portion of plume mass, 

released from a buoyant source, rises to and remains near the top of the boundary layer before 

becoming vertically mixed throughout the CBL. The model also tracks any plume mass that 

penetrates into an elevated stable layer, and then allows it to re-enter the boundary layer when 

and if appropriate. 

In urban areas, AERMOD accounts for the dispersive nature of the “convective-like” 

boundary layer that forms during nighttime conditions by enhancing the turbulence over that 

which is expected in the adjacent rural, stable boundary layer. The enhanced turbulence is the 

result of the urban heat flux and associated mixed layer which are estimated from the urbanrural temperature difference as suggested by Oke (1978; 1982).

In complex terrain, AERMOD incorporates the concept of the dividing streamline 

(Snyder et al., 1985) for stably-stratified conditions. Where appropriate the plume is modeled as 
44

a combination of two limiting cases: a horizontal plume (terrain impacting) and a terrainfollowing (terrain responding) plume. That is, AERMOD handles the computation of pollutant 

impacts in both flat and complex terrain within the same modeling framework. Generally, in 

stable flows, a two-layer structure develops in which the lower layer remains horizontal while 

the upper layer tends to rise over the terrain. The concept of a two-layer flow, distinguished at 

the dividing streamline height (Hc), was first suggested by theoretical arguments of Sheppard 

(1956) and demonstrated through laboratory experiments, particularly those of Snyder et al. 

(1985). In neutral and unstable conditions Hc = 0. 

A plume embedded in the flow below Hc tends to remain horizontal; it might go around 

the hill or impact on it. A plume above Hc will ride over the hill. Associated with this is a 

tendency for the plume to be depressed toward the terrain surface, for the flow to speed up, and 

for vertical turbulent intensities to increase. These effects in the vertical structure of the flow 

are accounted for in models such as the Complex Terrain Dispersion Model (CTDMPLUS) 

(Perry 1992). However, because of the model complexity, input data demands for CTDMPLUS 

are considerable. EPA policy (Code of Federal Regulations 1997) requires the collection of 

wind and turbulence data at plume height when applying CTDMPLUS in a regulatory 

application. As previously stated, the model development goals for AERMOD include having 

methods that capture the essential physics, provide plausible concentration estimates, and 

demand reasonable model inputs while remaining as simple as possible. Therefore, AERMIC 

arrived at a terrain formulation in AERMOD that considers vertical flow distortion effects in the 

plume, while avoiding much of the complexity of the CTDMPLUS modeling approach. Lateral 

flow channeling effects on the plume are not considered by AERMOD.

AERMOD captures the effect of flow above and below the dividing streamline by 

weighting the plume concentration associated with two possible extreme states of the boundary 

layer (horizontal plume and terrain-following). As is discussed below, the relative weighting of 

the two states depends on: 1) the degree of atmospheric stability; 2) the wind speed; and 3) the 

plume height relative to terrain. In stable conditions, the horizontal plume "dominates" and is 

given greater weight while in neutral and unstable conditions, the plume traveling over the 

terrain is more heavily weighted.
45

5.1 General structure of AERMOD including terrain

In general, AERMOD models a plume as a combination of two limiting cases: a 

horizontal plume (terrain impacting) and a terrain-following plume. Therefore, for all 

situations, the total concentration, at a receptor, is bounded by the concentration predictions 

from these states. In flat terrain the two states are equivalent. By incorporating the concept of 

the dividing streamline height, in elevated terrain, AERMOD’s total concentration is calculated 

as a weighted sum of the concentrations associated with these two limiting cases or plume states

(Venkatram et al. 2001).

The AERMOD terrain pre-processor (AERMAP) uses gridded terrain data to calculate a 

representative terrain-influence height (hc) for each receptor with which AERMOD computes 

receptor specific Hc values. Through this approach, AERMOD handles the computation of 

pollutant impacts in both flat and elevated terrain within the same modeling framework thereby 

obviating the need to differentiate between the formulations for simple and complex terrain (as 

required with previous regulatory models) . 

The general concentration equation, which applies in stable or convective conditions is 

given by

(48)

where CT{xr,yr,zr} is the total concentration Cc,s{zr,yr,zr

} is the contribution from the horizontal 

plume state (subscripts c and s refer to convective and stable conditions, respectively), 

Cc,s{zr,yr,zr

} is the contribution from terrain-following state, f is the plume state weighting 

function, {x , y ,z } r r r is the coordinate representation of a receptor (with zr defined relative to 

stack base elevation) zp = zr - zt is the height of a receptor above local ground, and zt is the terrain 

height at a receptor. Note that in flat terrain, zt = 0, zp = zr and the concentration (eq. (48)) 

reduces to the form for a single horizontal plume. It is important to note that for any 

concentration calculation all heights (z) are referenced to stack base elevation. Figure 10

illustrates the relationship between the actual plume and AERMOD’s characterization of it.

C x y z f C x y z ( f )C x y z T r r r c s r r r c s r r p { , , } { , , } { , , } = , , ⋅ + 1−
46

Figure 10. AERMOD two state approach. The total concentration predicted by AERMOD is the 

weighted sum of the two extreme possible plume states

The formulation of the weighting factor requires the computation of Hc. Using the 

receptor specific terrain height scale (hc) from AERMAP, Hc is calculated from the same 

algorithms found in CTDMPLUS as:

(49)

where u {Hc} is the wind speed at height Hc, and NN = �

gg

θθ

∂∂

∂∂ �

1

2�

is the Brunt-Vaisala frequency. 

The height scale, hc, characterizes the height of the surrounding terrain that most dominates the 

flow in the vicinity of the receptor.

1 2 ( ) 2 2 / ⋅ { } = − . ∫ u H N h z dz c c

H

h

c

c
47

The weighting between the two states of the plume depends on the relationship between 

Hc and the vertical concentration distribution at the receptor location. Assuming that the wind 

speed increases with height, Hc can be thought of as the level in the stable atmosphere where the 

flow has sufficient kinetic energy to overcome the stratification and rise to the height of the 

terrain. However, in determining the amount of plume material in the terrain-following state at 

a receptor, it is only important to know the lowest height in the flow where the kinetic energy is 

sufficient for a streamline to just maintain its height above the surface, i.e. terrain-following. 

Whether it will be deflected further and reach the top of some specified hill is not important for 

determining the amount of plume material in the terrain-following state for this receptor. 

Venkatram et al. (2001) first proposed the idea that for real terrain, often characterized by a 

number of irregularly-shaped hills, Hc should be defined in relation to a terrain-following height 

at each receptor location. This is in contrast to the more classical definition where Hc is defined 

in relation to the top of a single representative hill upon which may reside many receptor 

locations.

In the AERMOD approach, plume height, receptor elevation, and Hc will determine how 

much plume material resides in each plume state. For a receptor at elevation zt and an effective 

plume at height he the height that the streamlines must reach to be in the terrain-following state 

is zt+he. Therefore the terrain height of importance, hc, in determining Hc is simply equal to this 

local terrain-following height. Any actual terrain above hc = zt + he is of no consequence to the 

concentration at the receptor. This receptor and plume dependent approach to computing Hc

assumes that there is sufficient terrain affecting the flow near the receptor to vertically force the 

streamlines to the terrain-following level. If the actual surrounding terrain does not reach the 

height of the terrain-following state, hc is calculated from the highest actual terrain height in the 

vicinity of the receptor. Therefore, for any receptor, hc is defined as the minimum of the highest 

actual terrain and the local terrain-following height. Given hc, the dividing streamline height is 

computed with the same integral formula found in the CTDMPLUS model. 

The fraction of the plume mass below Hc (i.e., φp) is computed as:
48

(50)

where Cs{xr.yr,zr} is the concentration in the absence of the hill for stable conditions. In 

convective conditions Hc = 0 and φp = 0. As described by Venkatram et al. (2001), the plume 

state weighting factor f is given by f = 0.5(1 + φp). When the plume is entirely below Hc 

(φp = 1.0 and f =1.0) the concentration is determined only by the horizontal plume. When the 

plume is entirely above the critical dividing streamline height or when the atmosphere is either 

neutral or convective, (φp = 0 and f =0.5). Therefore, during convective conditions the 

concentration at an elevated receptor is simply the average of the contributions from the two 

states. As plumes above Hc encounter terrain and are deflected vertically, there is also a 

tendency for plume material to approach the terrain surface and to spread out around the sides of 

the terrain. To simulate this the estimated concentration is constrained to always contain a 

component from the horizontal state. Therefore, under no conditions is the plume allowed to 

completely approach the terrain-following state. For flat terrain, the contributions from the two 

states are equal, and are equally weighted. 

Figure 11 illustrates how the weighting factor is constructed and its relationship to the 

estimate of concentration as a weighted sum of two limiting plume states.

φ p

s r r r

H

s r r r

C x y z dz

C x y z dz

c

= ∫

∫

∞

{ , , }

{ , , }

0

0
49

Figure 11. Treatment of Terrain in AERMOD. Construction of the weighting factor used in 

calculating total concentration

The general form of the expressions for concentration in each term of eq. (48) for both 

the CBL and the SBL can be written as follows:

(51)

where Q is the source emission rate, uu� is the effective wind speed, and py and pz are probability 

density functions (pdf) which describe the lateral and vertical concentration distributions, 

respectively. AERMOD assumes a traditional Gaussian pdf for both the lateral and vertical 

distributions in the SBL and for the lateral distribution in the CBL. The CBL’s vertical 

distribution of plume material reflects the distinctly non-Gaussian nature of the vertical velocity 

distribution in convectively mixed layers. The specific form for the concentration distribution in 

the CBL is found in eq. (54) which uses the notation Cc {xr, yr, zr}. Similarly, in the SBL, the 

C x y z (Q u) P y x P z x y z { , , } ~ = { ; } { ; },
50

concentration takes the form of eq. (67) and used the notation Cs {xr, yr, zr}.

AERMOD simulates five different plume types depending on the atmospheric stability 

and on the location in and above the boundary layer: 1) direct, 2) indirect, 3) penetrated, 4) 

injected and 5) stable. All of these plumes will be discussed, in detail, throughout the remainder 

of this document. During stable conditions, plumes are modeled with the familiar horizontal 

and vertical Gaussian formulations. During convective conditions (L<0) the horizontal 

distribution is still Gaussian; the vertical concentration distribution results from a combination 

of three plume types: 1) the direct plume material within the mixed layer that initially does not 

interact with the mixed layer lid; 2) the indirect plume material within the mixed layer that rises 

up and tends to initially loft near the mixed layer top; and 3) the penetrated plume material that 

is released in the mixed layer but, due to its buoyancy, penetrates into the elevated stable layer.

During convective conditions, AERMOD also handles a special case referred to as an 

injected source where the stack top (or release height) is greater than the mixing height. 

Injected sources are modeled as plumes in stable conditions, however the influence of the 

turbulence and the winds within the mixed layer are considered in the inhomogeneity 

calculations as the plume material passes through the mixed layer to reach receptors.

As described above, AERMOD accounts for the vertical variation of meteorology 

through the use of effective values of wind speed, turbulence, and the Lagrangian time scale. 

Being a steady state plume model, AERMOD uses a single value of each meteorological 

variable to represent the state of the dispersive layer for each modeling period (typically one 

hour). Specifically, the effective parameters are determined by averaging values from the 

meteorological profile within the layer between the plume’s center of mass and the receptor. 

Effective variables or parameters are denoted by an overbar tilde (e.g., uu�).

5.2 Concentration predictions in the CBL

In AERMOD, the dispersion formulation for the convective boundary layer (CBL) 

represents one of the more significant model advances by comparison with existing regulatory 

models. One assumes that plume sections are emitted into a traveling train of convective 
51

elements - updrafts and downdrafts - that move with the mean wind. The vertical and lateral 

velocities in each element are assumed to be random variables and characterized by their 

probability density functions (pdf). The mean concentration is found from the pdf of the 

position of source-emitted “particles”; this position pdf in turn is derived from the pdf of the 

lateral and vertical velocities as described by Weil et al. (1997); also see Misra (1982), 

Venkatram (1983), and Weil (1988a).

In the CBL, the pdf of the vertical velocity (w) is positively skewed and results in a nonGaussian vertical concentration distribution, Fz (Lamb 1982). The positive skewness is 

consistent with the higher frequency of occurrence of downdrafts than updrafts; for an elevated 

non-buoyant source the skewness also leads to the decent of the plume centerline, as defined by 

the locus of maximum concentration (Lamb 1982; Weil 1988a). Figure 12 presents a schematic 

representation of an instantaneous plume in a convective boundary layer and its corresponding 

ensemble average. The base concentration prediction in AERMOD is representative of a one 

hour average. Notice that since a larger percentage of the instantaneous plume is effected by 

downdrafts, the ensemble average has a general downward trend. Since downdrafts are more 

prevalent the average velocity of the downdrafts is correspondingly weaker than the average 

updraft velocity to insure that mass is conserved. In AERMOD, a skewed vertical velocity pdf is 

modeled using a bi-Gaussian distribution, which has been shown to be a good approximation to 

laboratory convection tank data (Baerentsen and Berkowicz 1984). In contrast to the vertical 

component, the lateral velocity pdf is approximately Gaussian (Lamb 1982), and this pdf and 

the resulting concentration distribution, Fy, are assumed to be Gaussian.
52

Figure 12. Instantaneous and corresponding ensemble-averaged plume in the CBL

In addition to the non-Gaussian Fz, AERMOD has the following features. For buoyant 

releases, there is no “final” plume rise assumed. Instead, the plume or particle trajectories are 

determined by the addition of a distance-dependent plume rise and the random vertical 

displacement caused by the vertical distribution of w. Ground level concentrations first appear 

when the negative or downdraft velocities are sufficiently large to overcome the plume rise 

velocity and carry plume sections to the surface. The direct transport of plume material to the 

ground is treated by the “direct” source located at the stack. That is, the direct source treats that 

portion of the plume’s mass to first reach the ground, and all subsequent reflections of the mass 

at z = zi and 0 (where zi is the mixed layer height in the CBL (Cimorelli et al., 2004). For plume 

segments or particles initially rising in updrafts, an “indirect” or modified-image source is 

included (above the mixed layer) to address the initial quasi-reflection of plume material at z = zi, 

i.e., for material that does not penetrate the elevated inversion. This source is labeled “indirect” 

because it is not a true image source (i.e., as is found in models such as ISC) - the plume is not 
53

perfectly reflected about zi. Thus, the indirect source treats that portion of the plume’s mass that 

first reaches zi and all subsequent reflections of that particular mass at z = 0 and zi. For the 

indirect source, a plume rise (Δhi) is added to delay the downward dispersion of material from 

the CBL top (see Figure 13); this mimics the plume’s lofting behavior, i.e., the tendency of 

buoyant plumes to remain temporarily near zi and resist downward mixing. For non-buoyant 

sources the indirect source reduces to the first image source (as found in ISCST3) resulting from 

the first reflection at z = zi. Additionally, a “penetrated” source or plume (above the CBL top) is 

included to account for material that initially penetrates the elevated inversion but is 

subsequently reentrained by and disperses in the growing CBL.

Figure 13. AERMOD’s three plume treatment of the CBL

In line with the above concepts there are three main mathematical sources that contribute 

to the modeled concentration field: 1) the direct source (at the stack), 2) the indirect source, and 

3) the penetrated source. The strength of the direct source is fpQ, where Q is the source emission 
54

rate and fp is the calculated fraction of the plume mass trapped in the CBL (0 ≤ fp ≤ 1). Likewise, 

the indirect source strength is fpQ since this (modified image) source is included to satisfy the noflux boundary condition at z = zi for the trapped material. The strength of the penetrated source 

is (1- fp)Q, which is the fraction of the source emission that initially penetrates into the elevated 

stable layer. In addition to the three main sources, other image sources are included to satisfy the 

no-flux conditions at z = 0 and zi.

For material dispersing within a convective layer, the conceptual picture (see Figure 12) 

is a plume embedded within a field of updrafts and downdrafts that are sufficiently large to 

displace the plume section within it. The relationship between the particle (or air parcel) height,

zc and w is found by superposing the plume rise (Δh) and the vertical displacement due to w (i.e., 

wx/u), as

(52)

where hs is the stack height (corrected for stack tip downwash), u is the mean wind speed (a 

vertical average over the convective boundary layer) and x is the downwind distance. The Δh

above includes source momentum and buoyancy effects as given by eq. (91) below (see Briggs 

(1984)). The Fz or pdf of zc is found from the vertical velocity pdf pw as described in Weil et al. 

(1997). In the CBL a good approximation to pw is the superposition of two Gaussian 

distributions (Baerentsen and Berkowicz 1984; Weil 1988a) such that

(53)

where λ1 and λ2 are weighting coefficients for the two distributions with λ1 + λ2 = 1(the subscripts 

1 and 2 refer to the updraft and downdraft distributions, respectively). The parameters of the pdf 

(w1, w2, σw1, σw2, λ1, λ2) are functions of σw (the “total” or overall root mean square vertical 

turbulent velocity), the vertical velocity skewness SS = ww

3

σσww

3 � (where ww

3 is the third moment of 

z h h wx

u c = s + ∆ + ,

( ) ( ) p

w w w w

w

w w w w

= − − 













 + − − 















λ

πσ σ

λ

πσ σ

1

1

1

2

1

2

2

2

2

2

2

2 2 2 2 2 exp exp ,
55

w), and a parameter RR = σσww1⁄ww1 = − σσww2⁄ww2 = 2. An expanded discussion of the pdf 

parameters is given in Weil et al. (1997).

The instantaneous plume is assumed to have a Gaussian concentration distribution about 

its randomly varying centerline. The mean or average concentration is found by summing the 

concentrations due to all of the random centerline displacements. This averaging process results 

in a skewed distribution which AERMOD represents as a bi-Gaussian pdf (i.e., one for updrafts 

and the other for downdrafts). Figure 14 illustrates the bi-Gaussian approach to approximate the 

skewed vertical concentration distribution in the CBL. The figure shows two mean trajectories, 

each representing the average of many individual trajectories of parcels (or particles) released 

into downdrafts (the downdraft plume) or updrafts (the updraft plume). The velocities 

determining these mean trajectories are: 1) the mean horizontal wind speed (u), 2) the vertical 

velocity due to plume buoyancy (vbuoy), and 3) the mean updraft ( w1 ) or downdraft ( w2 ) velocity. 

The mean height of each trajectory, zzcc1 or zzcc2, can be found by averaging eq. (53). These parcel 

(or particle) height distributions are thus related to concentration and are characterized by σz1 

(= σw1x/u) and σz2 (= σw2x/u), the standard deviations of the two concentration distributions 

comprising the bi-Gaussian form as derived in Weil et al. (1997).
56

Figure 14. AERMOD’s pdf approach for plume dispersion in the CBL. AERMOD approximates 

the skewed distribution by superimposing two Gaussian distributions, the updraft and downdraft 

distributions

Figure 15 compares the bi-Gaussian pdf with the Gaussian form, which is symmetric 

about w = 0. As can be seen, for the negative and positive tails of the distributions, the biGaussian pdf is biased towards smaller and larger pw values, respectively, than the Gaussian. In 

addition, for the bi-Gaussian forms, approximately 60% of the area under the pw curve is on the 

negative side of the w axis and approximately 40% on the positive side. This is consistent with 

the results of numerical simulations and field observations (Lamb 1982; Weil 1988a).
57

Figure 15. Probability density functon of the vertical velocity. While the Gaussian curve is 

unskewed, the bi-Gaussian curve has a skewness of S=1

In the pdf approach used here (Weil et al. 1997), there are, as mentioned in the previous 

section, three primary sources that contribute to the modeled concentration field: 1) the “direct” 

or real source at the stack, 2) an “indirect” source that the model locates above the CBL top to 

account for the slow downward dispersion of buoyant plumes that “loft” or remain near, but 

below, zi , and 3) a “penetrated source” that contains the portion of plume material that has 

penetrated into the stable layer above zi. The direct source describes the dispersion of plume 

material that reaches the ground directly from the source via downdrafts. The indirect source is 

included to treat the first interaction of the “updraft” plume with the elevated inversion - that is, 

for plume sections that initially rise to the CBL top in updrafts and return to the ground via 

downdrafts. Image sources are added to treat the subsequent plume interactions with the ground 

and inversion and to satisfy the zero-flux conditions at z = 0 and at z = zi. This source plays the 

same role as the first image source above zi in the standard Gaussian model, but differs in the

treatment of plume buoyancy. For the indirect source, a modified reflection approach is adopted 
58

in which the vertical velocity is reflected at z = zi, but an “indirect” source plume rise Δhi is 

added to delay the downward dispersion of plume material from the CBL top. This is intended 

to mimic the lofting behavior. The penetrated source is included to account for material that 

initially penetrates the elevated inversion but subsequently can reenter the CBL via turbulent 

mixing of the plume and eventual reentrainment into the CBL. Figure 13 illustrates this three 

plume approach; a fundamental feature of AERMOD’s convective model. In AERMOD, the 

total concentration (Cc) in the CBL is found by summing the contribution from the three sources. 

For the horizontal plume state, the Cc is given by

(54)

where Cd, Cr, and Cp are the contributions from the direct, indirect and penetrated sources, 

respectively. The total concentration for the terrain-following state has the form of eq. (54) but 

with zr replaced by zp.

The fraction fp of the source material that remains trapped in the CBL is found from

(55)

where Δhh = zi - hs, and Δheq is the equilibrium plume rise in a stable environment. The Δheq has 

the form Berkowicz et al. (1986)

(56)

where PPss = FFbb uuNNh

2Δhh

3 ⁄ is the penetration parameter, and the stack buoyancy flux (Fb), and 

Brunt-Vaisala wfrequency (Nh) are given respectively by

C x y z C x y z C x y z C x y z c r r r d r r r r r r r p r r r { , , } = { , , }+ { , , }+ { , , },

f if h h

f if h h

f h

h if h h h

p h eq

p h eq

p

h

eq

eq h eq

= <

= >

= − < >

0 0 5

1 15

0 5 0 5 15

∆ ∆

∆ ∆

∆

∆

∆ ∆ ∆

.

.

. . . .

∆h ( P ( ) ) ∆h eq = s + h 2 6 2 3 3 3 1 3

.
59

(57)

and

(58)

Here, u is the wind speed at stack height; g is the gravitational acceleration; ws, rs,, and Ts

are the stack exit velocity, radius, and temperature, respectively; and θ is the ambient potential 

temperature. The Nh in eq. (58) is based on the potential temperature gradient in the elevated 

stable layer, provided by AERMET, capping the CBL. In general this layer is within zi and 

zi + 500 m.

5.2.1 Direct Source contribution to concentration calculations in the CBL

Following Weil et al. (1997), the concentration due to the direct plume is given by:

(59)

where

(60)

u is the wind speed at stack top, FFyy �= 1

√2ππσσyy

exp �

−yy2

2σσyy

2�� the lateral distribution function with 

meander (discussed in Section 0), wwjj = aajjww∗ (aj is defined below in eq. (62), Δhd is the direct 

source plume rise calculated from eq. (91), and z = zr and zp in the horizontal and terrainfollowing states, respectively. Here, Ψdj and σzj are the effective source height and vertical 

F gw r

T

T b s s

s

= 2 ∆

{ } N g

z z

h

i z zi

= 











 θ > 

∂θ

∂

1 2

.

( ) ( ) C x y z

Qf

u

F

z mz z mz

d r r

p

y

j

zj

dj i

zj

dj i

m zj j

{ , , } ~ = ⋅ exp − exp , − − 















+ − + + 





























 =

∞

=

∑ ∑ 2

2

2

2

2

2

2

2

2

1 0

2

π

λ

σ σ σ

Ψ Ψ

Ψdj s ∆ d

j h h

w x

u = + + ,
60

dispersion parameter corresponding to each of the two distributions in eq. (53). The subscript j is 

equal to 1 for updrafts and 2 for downdrafts. The lateral and vertical dispersion parameters (σy

and σzj), resulting from the combined effects of ambient, buoyancy-induced, and buildinginduced turbulence are calculated as discussed in Sections 0 and 0 respectively. Here, σzj (with j

= 1 or 2) is the vertical dispersion parameter corresponding to each of the Gaussian distributions 

used in the bi-Gaussian pdf, (see Section 0) and λj, the weighting coefficient for each distribution 

in eq.(53), is calculated from Weil et al. (1997) as 

(61)

where

(62)

Recall that ~

σ wT is the total effective vertical turbulence and is calculated from eq. (34). 

The parameters appearing in eq. (62) are given by

(63)

where,

λ

λ

1

2

2 1

2

2 1

2

1

2 1

1

2 1

= + − = + −

= − − = − −

w

w w

a

a a

w

w w

a

a a

a

w

S S

a

w

S S

wT

wT

1

2 2

1 2

2

2 2

1 2

2

1

2

4

2

1

2

4

= + +





 





















= − +





 





















~

~

*

*

σ α

α

β

σ α

α

β

w

w

for H x z

w

w

H x

z

for H x z

p i

p

i

p i

3

3

3

3

0125 01

125 01

*

*

. { } .

. { } { } .

= ≥

= <
61

(64)

and R is assumed to be 2.0 (Weil et al., 1997). Likewise, the term wwjj xx⁄uu in eq. (60) follows 

from the Fz derivation and the wj appearing in the bi-Gaussian form (see discussion of eq. (53)). 

The lateral dispersion parameter (σy,) is calculated from eq. (75) (Weil et al., 1997).

In eq. (59), an image plume is used to satisfy the no-flux condition at the ground, i.e., an 

image plume from a source at z = -hs, which results in the exponential terms containing z +Ψdj

on the right-hand side of eq. (59). This image source results in a positive flux of material at z = 

zi, and additional image sources are introduced at z = 2 zi + hs, -2 zi - hs, 4zi + hs, -4zi - hs, etc. to 

satisfy all the subsequent no-flux conditions occurring at z = 0 and zi.

5.2.2 Indirect Source contribution to concentration calculations in the CBL

The concentration due to the indirect source is calculated from:

(65

where Ψrj = Ψdj - Δhi, and z is either zr.(for the horizontal plume state) or zp (for the terrainfollowing state). As shown in Figure 13, the indirect plume is modeled as a reflected version of 

the direct plume with an adjustment (Δhi - calculated from eq. (92)) to the reflected plume height 

to account for the delay in vertical mixing due to plume lofting at the top of the boundary layer. 

( )

α

β

σ

= +

+

= +

= ≡

1

1 3

1

2

2

2

3 3

3

R

R

R

S w w

w

Skewness factor

wT

*

*

~ ,

( ) ( ) C x y z

Qf

u

F

z mz z mz

r r r

p

y

j

zj

rj i

zj

rj i

m zj j

{ , , } ~ = ⋅ ⋅ exp − exp , + − 















+ − − + 





























 =

∞

=

∑ ∑ 2

2

2

2

2

2

2

2

2

1 1

2

π

λ

σ σ σ

Ψ Ψ
62

5.2.3 Penetrated source contribution to concentration calculations in the CBL

For the penetrated source the concentration expression has a Gaussian form in both the 

vertical and lateral directions. The concentration due to this source is given by:

(6

where zieff is the height of the upper reflecting surface in a stable layer (see Section 0) and z is 

either zr for the horizontal plume state or zp for the terrain-following state. The vertical 

dispersion parameters (σzp) are calculated as described in Section 0.

The penetrated plume height, hep, is taken as the height of the plume centroid above the 

mixed layer and is calculated from eq. (94).

5.3 Concentrations in the SBL

For stable conditions, the AERMOD concentration expression (Cs in eq. (48)) has the 

Gaussian form, and is similar to that used in many other steady-state plume models (e.g., HPDM 

(Hanna and Paine 1989)). The Cs is given by

(67)

where zieff is the effective mechanical mixed layer height, σzs is the total vertical dispersion in the 

SBL (see discussion in Section 0), and hes is the plume height (i.e., stack height plus the plume 

rise - see Section 0).

Above the mechanical mixed layer height, zim (eq. (26)), the turbulence level is generally 

expected to be small and thus supports little vertical mixing of the plume. AERMOD is 

( ) ( ) ( ) C x y z

Q f

u

F

z h mz z h mz

p r r

p

zp

y

ep ieff

zp

ep ieff

m zp

{ , , } ~ = exp exp − ⋅ − − + 















+ − + + 





























 =−∞

∞

∑

1

2

2

2

2

2

2

2

2

2 π σ σ σ

( ) ( ) C x y z Q

u

F

z h mz z h mz

s r r

zs

y

es ieff

zs

es ieff

m zs

{ , , } ~ = ⋅ ⋅ exp − exp , − − 















+ − + + 





























 =−∞

∞

∑ 2

2

2

2

2

2

2

2

2 π σ σ σ
63

designed (in the SBL) with an effective mixing lid, zieff, that retards but does not prevent plume 

material from spreading into the region above the estimated mechanical mixed layer. When the 

final plume height is well below zim, the plume does not interact with zim. When the plume is 

below zim yet the “upper edge” (plume height plus 2.15 σzs) of the stabilized plume reaches zim, 

the effective mixing lid is allowed to increase and remain at a level near the upper edge of the 

plume. In this way, AERMOD allows the plume to disperse downwards, but where the 

turbulence aloft is low, vertical plume growth is limited by an effective reflecting surface that is 

folding back only the extreme tail of the vertical plume distribution. There is no strong 

concentration doubling effect as occurs with reflections from an assumed hard lid. Downward 

dispersion is primarily a factor of σw averaged from the receptor to the plume height. If the 

plume height is above the mixed layer height, the calculation of the effective σw will include 

regions in which σw is likely to be small. This, in effect, retards plume growth by an amount 

dependent upon how much of the plume is above zim. Therefore, whether the plume is above or 

below zim, the region of low turbulence above zim will have an appropriate effect on the 

concentration distribution within the mixing layer.

When the plume buoyancy carries the rising plume into the relatively non-turbulent layer 

above zim, the reflecting surface is still placed at 2.15 σzs above the effective plume height 

because there will be plume spread due to plume buoyancy and downward mixing is still 

important. Therefore, in the SBL, plume material is assumed to reflect off an elevated surface 

which is defined as:

(68)

where σzs in eq. (68) is determined from equations found in Section 0 with σwT and u evaluated at 

hes; not as an effective parameter. It is important to note that zieff depends on downwind distance 

since σzs is distance dependent. In fact, as eq. (68) suggests, this effective reflecting surface is 

only folding back the extreme tail of the upward distribution. Also, if the height of the receptor

zr ≥ zieff then the effective reflecting surface is not considered. This approach is also implemented 

for the penetrated source. For the penetrated and injected sources zieff is calculated using eq. (68)

with σzs and hes replaced by σzp and hep respectively.

zieff = MAX[(hes + 2.15σ zs{hes}; zim )].
64

5.4 Treatment of lateral plume meander

In AERMOD we include the effect that lower-frequency, non-diffusing eddies (i.e., 

meander) have on plume concentration. Meander (or the slow lateral back and forth shifting of 

the plume) decreases the likelihood of seeing a coherent plume after long travel times. This 

effect on plume concentration could best be modeled with a particle trajectory model, since 

these models estimate the concentration at a receptor by counting the number of times a particle 

is seen in the receptor volume. However, as a simple steady state model, AERMOD is not 

capable of producing such information. AERMOD accounts for meander by interpolating 

between two concentration limits: the coherent plume limit (which assumes that the wind 

direction is distributed about a well-defined mean direction with variations due solely to lateral 

turbulence) and the random plume limit, (which assumes an equal probability of any wind 

direction). 

For the coherent plume, the horizontal distribution function (FyC) has the familiar 

Gaussian form:

(69)

where σy is the lateral dispersion parameter (see Section 0). For the random plume limit, the 

wind direction (and plume material) is uniformly distributed through an angle of 2π. Therefore, 

the horizontal distribution function FyR takes the simple form:

FFyyyy = 1

2ππxxrr

(70)

where xr is radial distance to the receptor. Although the form of the vertical distribution function 

remains unchanged for the two plumes, its magnitude is based on downwind distance for the 

coherent plume and radial distance for the random plume.

FyC

y

y y

=  −















1

2 2

2

2 πσ σ

exp
65

Once the two concentration limits (CCh - coherent plume; CR - random plume) have been 

calculated, the total concentration for stable or convective conditions (Cc,s) is determined by 

interpolation. Interpolation between the coherent and random plume concentrations is 

accomplished by assuming that the total horizontal “energy” is distributed between the wind’s 

mean and turbulent components. That is,

(71)

where σh

2 is a measure of the total horizontal wind energy and σr2 is a measure of the random 

component of the wind energy. Therefore, the ratio σr

2

/σh

2 is an indicator of the importance of 

the random component and can therefore be used to weight the two concentrations as done in eq. 

(71).

The horizontal wind is composed of a mean component uu, and random components σu

and σv. Thus, a measure of the total horizontal wind “energy” (given that the alongwind and 

crosswind fluctuations are assumed equal i.e., σu = σv), can be represented as 

(72)

where uu = (uu�2 − 2σσ�vv

2)1/2. The random energy component is initially 2σσ�vv2 and becomes equal to 

σh

2 at large travel times from the source when information on the mean wind at the source 

becomes irrelevant to the predictions of the plume’s position. The evolution of the random 

component of the horizontal wind energy can be expressed as

σσrr

2 = 2σσ�vv2 + uu

2 �1 − exp �

−xxrr

uu�TTrr � �� (73)

where Tr is a time scale (= 24 hrs) at which mean wind information at the source is no longer 

correlated with the location of plume material at a downwind receptor. Analyses involving 

autocorrelation of wind statistics (Brett and Tuller 1991) suggest that after a period of 

approximately one complete diurnal cycle, plume transport is “randomized.” Equation (73)

Cc,s = CCh (1− r h ) + CR( r h ) 2 2 2 2 σ σ σ σ

σ h σ v u

2 2 2 = 2 + ~
66

shows that at small travel times, σσrr

2 = 2σσ�vv2, while at large times (or distances) σσrr2 = 2σσ�vv2 + uu

2

, 

which is the total horizontal kinetic energy (σh

2

) of the fluid. Therefore, the relative 

contributions of the coherent and random horizontal distribution functions (eq. (71)) are based on 

the fraction of random energy contained in the system (i.e., σσrr

2 σσh

2 ⁄ ).

The application of eq. (71) is relatively straight forward in the SBL. Since 

concentrations in the SBL are represented as a single plume, Cs can be calculated directly from 

eq. (71). By contrast for convective conditions the situation is complicated by the inclusion of 

plume penetration. Since σr

2 depends on the effective parameters (eq. (73)), the concentration 

weighting factors found in eq. (71) will be different for the non-penetrated and penetrated 

plumes of the CBL. This is handled by combining the penetrated and non-penetrated weighting 

factors (σr

2

/σh

2

|P and σr

2/σh2

|NP) into a single effective factor (σr

2

/σh

2

|CBL). That is,

(74)

where fp (see eq. (55)) is the fraction of the source material that remains trapped in the CBL. 

Using eq. (74), concentrations in the CBL (Cc) are calculated from eq. (71) with (σr

2

/σh

2

) 

replaced by (σr

2

/σh

2

|CBL).

5.5 Estimation of dispersion coefficients

The overall standard deviations (σy,z) of the lateral and vertical concentration 

distributions are a combination of the dispersion (represented by σya, σza) resulting from ambient 

turbulence, and dispersion (σb) from turbulence induced by plume buoyancy. Building induced 

dispersion is not included here since a separate approach (see Section 0) is taken for situations 

in which building wake effects contribute to the total dispersion. Dispersion induced by 

ambient turbulence is known to vary significantly with height, having its strongest variation 

near the earth’s surface. Unlike present regulatory models, AERMOD has been designed to 

account for the effect of variations of turbulence with height on dispersion through its use of 

“effective parameters” (see Section 0), which are denoted by an overscript tilde, e.g., σσ�wwww.

( ) σ

σ

σ

σ

σ

σ r

h CBL

p r

h P

p r

h NP

f f

2

2

2

2

2

2 = ⋅ + 1− ⋅
67

AERMOD treats vertical dispersion from ambient turbulence (σza) as a combination of a 

specific treatment for surface dispersion and the more traditional approach based on Taylor 

(1921) for elevated dispersion. Using this approach good agreement with observations was 

achieved in the SBL. However, the results in the CBL indicated that the treatment of lateral

dispersion near the surface was problematic. This problem was corrected through the 

development of an empirical relationship for σya near the surface using the full (CBL and SBL) 

Prairie Grass data set. A description of the resulting formulations for σya & σza is presented in 

the next section.

The approach used to combine the above contributions to dispersion assumes that the 

effects are independent of one another. Thus, the total dispersion coefficients, for situations that 

do not include building downwash effects, are calculated from the following general expression 

(Pasquill and Smith 1983): 

(75)

where the subscripts y and z are deleted from σb because σyb is assumed equal to σzb. With the 

exception of the CBL’s penetrated source the form of eq. (75) applies to all source dispersion in 

both the CBL and SBL such that σy,z becomes σys,zs and σyjs,zj and σya,za becomes σyas,zas and σyajs,zaj

for the SBL and CBL, respectively. For the penetrated source, the total dispersion is assumed to 

include ambient and buoyancy induced turbulence only; building wakes are assumed to have 

little influence. For the injected source, the total dispersion is calculated as if the source were in 

the SBL.

A comment on notation: eq. (75) applies for both lateral and vertical dispersion in the 

SBL and CBL. In references to the SBL, σz appears as σzs in the dispersion equation; σza appears 

as σzas. In reference to the CBL, σz appears as σzj for the dispersion expression applicable to the 

direct and indirect sources and σza appears as σaj; for the penetrated source σz appears as σzp in 

the dispersion expression. 

σ y,z σ ya,za σ b , 2 2 2 = +
68

5.5.1 Dispersion from ambient turbulence

5.5.1.1 Lateral dispersion from ambient turbulence

In general terms, the ambient component of the lateral dispersion is based upon Taylor

(1921) such that:

(76)

where p = 0.5, u is the wind speed, σv is the root-mean-square lateral turbulence velocity, and TLy

is the Lagrangian integral time for the lateral turbulence. Application of eq. (76) in a preliminary 

version of AERMOD yielded poor concentration estimates in comparison to those found in the 

Prairie Grass field experiments (Barad 1958). Specifically, the lateral spread was not well 

matched. Therefore, the lateral dispersion expression was reformulated to allow for an empirical 

fit to the Prairie Grass data.

Using an approach similar to that of Venkatram et al. (1984) TLy is found to be l/σv where 

l is an appropriate length scale for lateral turbulence. Equation (76) can be written in terms of 

the non-dimensional downwind distance X and a non-dimensional height scale α as:

(77)

where XX(= σσ�vv xx⁄uu� zz) is the non-dimensional distance with u and σv given by effective 

parameters, where α = zi /l, and zi is the mixed layer height.

Based on a preliminary comparison of σya (eq. (77)) with selected stable and convective 

cases from the Prairie Grass experiment (Barad 1958) α was found equal to 78 and p equal to 

σ

σ

ya

v

x

u

x u

T

Ly

p =

+

















~

~ ~

1

2

( )

σ

σ

α

ya v

p

x

u X = +

~

~ 1
69

0.3. As such, α is treated as a fitting parameter. In later comparisons against the full Prairie 

Grass data set (Figure 16), eq. (77) tended towards the lower envelope of this widely scattered 

data (i.e., lateral dispersion estimates are on the lower end of the distribution of measurements). 

However, the preliminary values of α (= 78) and p (= 0.3) produced good agreement between 

AERMOD concentration predictions and observations (Brode 2002). Therefore, these 

preliminary values were retained in AERMOD, and eq. (77) applies for the calculation of σya for 

all plumes in both the SBL and CBL.

Figure 16. Lateral spread (Fy) as a function of non-dimensional distance (X). The data is taken 

from the Prairie Grass experiment (Barad, 1958)

The ambient component of the lateral dispersion for the penetrated source, i.e. a source 

which has been released below zi, but penetrates above, is calculated using eq. (77) with hes set 

equal to hep (the height of the penetrated source). However, for the injected source, i.e. source 

released above zi, no substitution is needed since these sources are modeled as a stable source.
70

To account for the increase in the turbulence length scale and hence the Lagrangian time 

scale with release heights greater than that at Prairie Grass, α is scaled as follows:

(78)

where zPG = 0.46 m (Prairie Grass release height), and zzmmmmmm = max(zz; zzPPPP). To insure that α

does not become unrealistically large for surface releases, z is not allowed below zPG (i.e., 0.46 

m). In the SBL, z = hes; in the CBL z = hs; for penetrated sources, z = hep. . As α becomes small 

for large release heights, σya would tend to grow linearly with downwind distance.

5.5.1.2 Vertical dispersion from ambient turbulence

For sources in the SBL (and for sources in the CBL that are emitted directly into the 

stable layer above the mixed layer), the ambient portion of the vertical dispersion (σzas) is 

composed of an elevated (σzes) and near-surface (σzgs) component. For hes < zi simple 

interpolation provides a smooth transition between the two components.

(79)

For hes ≥ zi σzas is set equal to σzes. The expressions for calculating hes are found in 

Section 0. It should be noted, for sources in the SBL, that σzas is the specific form of the ambient 

portion of the vertical dispersion (i.e., σza in eq. (75)).

In the SBL, the elevated portion of the vertical dispersion follows the form of eq. (76):

(80)

α = 



 

 78  z

z

PG

max

σ zas σ σ es

i

zgs es

i

zes

h

z

h

z = − 



 



 +





 

 1  .

σ zes σ wT ( )

Lzs

x u

x u

T = +





 



 ~ ~ ~ 1 ' 2

1 2
71

where σwT is the vertical turbulence due to the mechanical mixing (Cimorelli et al., 2004).

As with the lateral component, the Lagrangian time scale (TLzs) for the vertical 

turbulence can be written in the form (Venkatram et al. 1984)

(81)

The length scale l is an interpolation between the limiting length scales for neutral 

conditions and stable conditions

(82)

where llnn = 0.36 h and llss = 0.27 σσ�wwww⁄NN. Under very stable conditions or at large heights, l

approaches ls. When conditions are near neutral, N is very small and l approaches ln.

By combining eqs. (80), (81), and (82) we find the following expression that is used by 

AERMOD to compute σzes, the elevated portion of the vertical dispersion for the stable source:

(83)

Finally, to complete the description of eq. (79), the surface portion of vertical dispersion 

(σzgs) in the SBL, is calculated from Venkatram (1992) as

(84)

T l

Lzs

wT

= ~

σ

1 1 1

l l l n s

= + .

σ

σ

σ

σ

zes

wT

wT

es wT

t

t

h

N

=

+ +





 



 



 





~

~

. . ~

.

1

2

1

0 36 0 27

1 2

σ

π zgs

u x

u

x

L = 













 + 



 





− 2 1 0 7

1 3

*

~

/

.
72

For the direct and indirect sources in the CBL, the ambient portion of the vertical 

dispersion (σza of eq. (75)) is denoted as σzaj (j = 1, 2) to distinguish between updrafts and 

downdrafts. σzaj is composed of an elevated (σzej) and surface (σzg) portion and is given by

(85)

where the elevated portion (σzej) is obtained from Weil et al. (1997)as 

(86)

where σwj is a parameter in the bi-Gaussian pdf (eq.(53)). 

The expression αb = min(0.6 + 4Hp/zi, 10.0) designed to be 1.0 above the surface layer 

(Hp > 0.1 zi) and to otherwise match Venkatram’s (1992) result for vertical dispersion from a 

surface source in a neutral boundary layer.

For the CBL, the vertical dispersion from a source within the surface layer (Hp{x} < 0.1

zi) is parameterized by

(87)

where bc = 0.5, u* is the friction velocity, and L is the Monin-Obukhov length; above the surface 

layer (Hp > 0.1zi), σzg is assumed to equal zero. In the limit of a surface release (Hp = 0), the 

parameterization of eq. (87) follows the form suggested by Venkatram (1992) for vertical 

dispersion in the unstable surface layer; i.e., σσzz ∝ (uu∗⁄uu�)2xx2/|LL| . The parameterization is 

designed to: 1) agree with Venkatram’s result in the limit of a surface release, 2) provide good 

agreement between the modeled and observed concentrations from the Prairie Grass experiment 

(Perry et al., 2005), and 3) decrease with source height in the surface layer and ultimately vanish 

σ zaj σ zej σ zg

2 2 2 = + ,

σ α

σ

zej b

wj x

u = ~ ,

σ zg c ( ) ( ) p

i

b H

z = − u u x L 



 



 



 



 1 10 ⋅ ⋅

2 2

* ~
73

for above the surface layer. The constant bc was chosen to satisfy the second design 

requirement. In the limit of a neutral boundary layer σzg is equal to zero.

The total vertical dispersion for the penetrated source σzp (= σz in eq. (75)) is a 

combination of both ambient and buoyancy effects. The ambient portion of the vertical 

dispersion for the penetrated source contains only an elevated component σzes (= σzss) since it is 

assumed to be decoupled from the ground surface by its location above zi and therefore 

unaffected by the underlying surface. The ambient vertical dispersion for the penetrated source 

is computed as the elevated portion of a stable source (σzes of eq. (83)) with N = 0 and with no 

contribution from the surface component. The Brunt-Vaisala frequency, N, is set to zero 

because the penetrated plume passes through the well mixed layer (where N ≈0) prior to 

dispersing to receptors within the mixed layer.

5.5.2 Buoyancy induced dispersion (BID) component of σy and σz

For all plumes, the buoyancy induced dispersion (BID) is calculated following Pasquill 

(Pasquill 1976) and Weil (1988b) as

(88)

where Δh is the plume rise appropriate for each of the plume types (direct, indirect, penetrated, 

and stable plumes). The direct source plume rise is calculated from eq. (91), stable plume rise 

(Δhs ) is calculated from eq. (95) and the plume rise for the penetrated source Δhpp = h − hss

where hep is calculated from eq. (94)).

5.5.3 Treatment of building downwash

AERMOD incorporates the Plume Rise Model Enhancements (PRIME) (Schulman et al. 

2000) algorithms for estimating enhanced plume growth and restricted plume rise for plumes 

affected by building wakes (U.S. Environmental Protection Agency, 1995). PRIME partitions 

plume mass between a cavity recirculation region and a dispersion enhanced wake region based 

σ b

h = 0 4

2

. , ∆
74

upon the fraction of plume mass that is calculated to intercept the cavity boundaries. These 

boundaries are established from estimates of the locations of the lateral and vertical separation 

streamlines. Dispersion of the recirculated cavity mass is based on building geometry and is 

assumed to be uniformly mixed in the vertical. At the boundary of the cavity region, cavity 

mass is emitted into the wake region. Here, it is combined with plume mass that was not 

captured by the cavity and dispersed at an enhanced rate based on source location, release height 

and building geometry. The enhancement of turbulence within the wake decays gradually with 

distance, allowing for a smooth transition to ambient levels of turbulence in the far-field. A 

probability density function model and an eddy diffusivity model (Weil 1996) are used for 

dispersion estimates in the near-wake and far-wake regions, respectively. Plume rise, for 

sources influenced by a building, is estimated using a numerical model that includes effects 

from streamline deflection near the building, vertical wind speed shear, enhanced dilution from 

the turbulent wake and velocity deficit. In general, these building induced effects act to restrict 

the rise that the plume would have in the absence of the building.

PRIME was originally designed (Schulman et al., 2000) to enhance plume growth using 

Pasquill Gifford (PG) dispersion (Pasquill 1961; Gifford 1961). AERMOD’s estimate of plume 

growth is based on dispersion parameters derived from profiles of turbulence (see Section 0), 

not from radiation base turbulence surrogates as is done in the PG approach. A basic design 

tenet for incorporating PRIME into AERMOD was to be as faithful as possible to the PRIME 

formulation while ensuring that 1) AERMOD’s ambient dispersion was used in place of PG 

dispersion and 2) far beyond the wake region, where building influences should be insignificant, 

concentrations approach the AERMOD estimate. Therefore, within the wake, PRIME 

algorithms are use exclusively to calculate concentration with AERMOD-derived ambient 

turbulent intensities as input. To insure a smooth transition between concentrations estimated 

by PRIME, within the wake, and AERMOD estimates in the far field, concentrations beyond the 

wake are estimated as the weighted sum of the two calculations. That is, beyond the wake the 

total concentration (Ctotal) is calculated as follows:

CTotal C ( )C (89) = ime + − AERMOD γ γ Pr 1
75

where Cprime is the concentration estimated using the PRIME algorithms with AERMOD-derived 

meteorological inputs, CAERMOD is the concentration estimated using AERMOD without 

considering building wake effects, and γ the weighting parameter. The weighting parameter, γ, is 

designed such that the contribution from the PRIME calculation decreases exponentially with 

vertical, lateral and downwind distance from the wake. It is calculated as follows:

(90)

where x is the downwind distance from the upwind edge of the building to the receptor, y is the 

lateral (crosswind) distance from the building centerline to the receptor, z is the receptor height 

above ground, σxg is longitudinal dimension of the wake, σyg is the distance from the building 

centerline to lateral edge of the wake, and σzg is the height of the wake at the receptor location.

5.6 Plume rise calculations in AERMOD

5.6.1 Plume rise in the CBL

The plume rise for the direct source is given by the superposition of source momentum 

and buoyancy effects following Briggs (1984).

(91)

where Fm = (T/Ts)ws

2

rs

2 the stack momentum flux, Fb = gws rs2 (ΔT/Ts) is the stack buoyant flux, 

rs is the stack radius corrected for stack tip downwash, and β1 (= 0.6) is an entrainment 

parameter. It should be noted that up is the wind speed used for calculating plume rise. In the 

CBL up is set equal to u{hs}.

( ) ( ) ( ) γ

σ

σ

σ

σ

σ

σ = − − 















− − 















− − 













 exp exp exp

x y z xg

xg

yg

yg

zg

zg

2

2

2

2

2

2 2 2 2

∆h F x

u

F x

u

d m

p

b

p

= + ⋅ 















3 3

1 2 2 21

2

2

3

1 3

β β
76

As shown in Figure 13, the indirect plume, which is included to treat the no flux 

condition at z = zi, is modeled as a reflected version of the direct plume with an adjustment 

(Δhi) to the reflected plume height to account for the delay in vertical mixing due to plume 

lofting at the top of the boundary layer. That height adjustment is given by

(92)

where ry and rz are the lofting plume half-widths in the lateral and vertical directions, up is the 

wind speed used for plume rise, and αr = 1.4. The produce of cross-wind dimensions of the 

assumed elliptical plume is calculated from Weil et al. (1997) as

(93)

where rrh = ββ2(zz − hss), β2 = 0.4, λy = 2.3, and ae = 0.1 (dimensionless entrainment parameter). 

For a derivation and discussion of Δhi see Weil et al. (1997).

The height that the penetrated source achieves above zi is calculated as the equilibrium 

plume rise in a stratified environment and is determined by the source buoyancy flux, the stable 

stratification above zi, and the mean wind speed. In line with Weil et al. (1997), the penetrated 

source plume height, hep, is taken as the centroid of plume material above the inversion. For 

complete penetration (fp = 0) hep = hs+Δheq. However, for partial penetration ( fp > 0), hep is 

chosen as the average of the heights of the upper plume edge hs + 1.5 Δheq and zi, or

(94)

where Δheq is defined in eq.(56).

∆h F z

u r r

x

u

i

b i

r p y z p

=

















2

1 2

α ,

r r r

a w x

u y z h

e y

p

= + ⋅ 2

3 2 2 2

2 4

λ *

h h z h ep s i = eq

+

+

2

0.75∆ .
77

5.6.2 Plume rise in the SBL

Plume rise in the SBL is taken from Weil (1988b), which is modified by using an 

iterative approach which is similar to that found in Perry et al. (1989). When a plume rises in an 

atmosphere with a positive potential temperature gradient, plume buoyancy decreases because 

the ambient potential temperature increases as the plume rises; thus, plume buoyancy with 

respect to the surroundings decreases. To account for this, the plume rise equations have to be 

modified. With this modification, AERMOD computes stable plume rise, Δhs, from Weil et al. 

(1988b) as

(95)

where Nʹ = 0.7N with N given by eq. (58). N and u are evaluated initially at stack height. Once 

plume rise has been computed, subsequent plume rise estimates are made (iteratively until 

convergence) by averaging the u and N values at stack top with those at hs + Δhs/2. Equation

(95) is used for downwind distances that are less than the distance to final rise (xf). Beyond xf, 

Δhs remains constant. The distance at which the stable plume reaches its maximum rise is given 

by

(96)

Upon substituting eq. (95) for x in eq. (97) the maximum final rise of the stable plume 

Δhs{xf} reduces to:

(97)

∆h F

N u

N F

F

N x

u

N x

u s b

p

m

b p p

= 













 ⋅ 













 + − 































2 66 1 2

1 3 1 3

. sin cos ,

' ' '

x

u

N

F N

F f

p m

b

= ′

 − ′



 

 arctan .

∆h x

F

u N s f

b

p

{ } = . . 













 2 66 2

1 3
78

As with eq. (95), the velocity, up, and N in eq. (97) are evaluated initially at stack height 

and then iteratively.

When the atmosphere is close to neutral, the Brunt Vaisala frequency, N, is close to zero, 

and eq. (95) can predict an unrealistically large plume rise. Under, these circumstances, plume 

rise is limited by atmospheric turbulence. This happens when the rate of plume rise under 

neutral conditions is comparable to σw. Under these conditions, stable plume rise (eq. (97)) is 

limited by the neutral rise calculated from Weil (1985) as 

(98)

where the neutral length scale LLnn = FFbb (uuppuu∗

2 ⁄ ).

As the wind speed approaches zero, eq. (95) again predicts unrealistic values. In these 

near-calm conditions the stable plume rise (eq. (97)) is limited by the calm rise expression that 

is based on the work of Morton et al. (1956) and Briggs (1969) such that,

(99)

Finally, the stable plume rise is limited by a calculation of the unstable rise (see Section 

0).

5.7 Source characterization

AERMOD gives the user the ability to characterize a source as either a point, an area, or 

a volume. AERMOD additionally has the capability of characterizing irregularly shaped area 

sources.

Point sources are characterized exactly as in the ISC3 model (U.S. Environmental 

Protection Agency, 1995). The input to the model includes the location, elevation, emission rate, 

∆hs =12 Ln (hs +12 Ln ) 3 5 2 5 . .

∆h F

N s b = 4 1

4

3

4

.
79

stack height, stack gas temperature, stack gas exit velocity, and stack inside diameter. The 

temperature, exit velocity, and diameter are required for plume rise calculations.

Similarly, volume sources require the same input as the ISC3 model. This includes the 

location, elevation height (optional), height of release, emission rate, the initial lateral plume 

size (σy) and initial vertical plume size (σz). AERMOD differs from ISC3 in the treatment of 

volume sources only in how the initial plume size is implemented. Where ISC3 uses the virtual 

source technique to account for initial plume size, AERMOD adds the square of the initial 

plume size to the square of the ambient plume size: 

(100)

where σyo is the initial horizontal plume size, σyl is the plume size before accounting for the initial 

size, and σy is the resultant plume size after accounting for the initial size.

The area source treatment is enhanced from that available in ISC3. In addition to being 

input as squares or rectangles, area sources may be input as circles or polygons. A polygon may 

be defined by up to 20 vertices. A circle is defined by inputting its center location and radius. 

The AERMOD code uses this information to create an equivalent nearly-circular polygon of 20 

sides, with the same area as the circle.

As with ISC3, AERMOD allows for the calculation of a simple half-life decay.

5.8 Plume volume molar ration method (PVMRM)

PVMRM was first introduced in AERMOD in version 04300 as an option for modeling 

the conversion of NOx to NO2 in the presence of ozone. The implementation is based on the 

work of Hanrahan (1999) and adapted for AERMOD. Details regarding the formulation of the 

PVMRM option in AERMOD, and preliminary model evaluation results are available in U.S. 

Environmental Protection Agency (2015).

σ y σ yl σ yo

2 2 2 = +
80

5.8.1 Definition of plume volume

5.8.1.1 Total vs. relative dispersion

The PVMRM determines the conversion rate for NOx to NO2 based on a calculation of 

the NOx moles emitted into the plume, and the amount of O3 moles contained within the volume 

of the plume between the source and receptor. The dispersion algorithms in AERMOD and 

other steady-state plume models are based on the use of total dispersion coefficients, which are 

formulated to represent the time-averaged spread of the plume. A more appropriate definition of 

the volume of the plume for purposes of determining the ozone moles available for conversion 

of NOx is based on the instantaneous volume of the plume, which is represented by the use of 

relative dispersion coefficients, (Cole and Summerhays, 1979; Bange, 1991). The 

implementation of PVMRM in AERMOD is based on the use of relative dispersion coefficients 

to calculate the plume volume. Weil (1996 and 1998) has defined formulas for relative 

dispersion that are consistent with the AERMOD treatment of dispersion, and which can be 

calculated using meteorological parameters available within AERMOD.

5.8.1.2 Calculation of relative dispersion coefficients

The formula for relative dispersion combines the effects of buoyancy-induced 

turbulence, which should dominate close to the source, and ambient turbulence, which begins to

dominate further downwind. Since the travel time from the source to the receptor is important 

for defining relative dispersion, the relative dispersion coefficients are calculated based on the 

radial distance from source to receptor. Weil (1996 and 1998) assumes relative dispersion (σr)

to be isotropic, so that σrx = σry = σrz = σr. The relative dispersion (σr) due to the combined 

effects of buoyancy- induced turbulence (σrb) and ambient turbulence (σra) is parameterized as

follows:

σσrr = (σσrr

3 + σσrr

3 )1⁄3 (101)

The buoyancy-induced dispersion term, σrb, is calculated in AERMOD as
81

σσ

0.4 ΔΔh

√2

(102)

where Δh is the plume rise. Relative dispersion due to ambient turbulence, σra, is parameterized 

by

σσrr = aa1εε1/2tt3/2

1 + aa2tt/TTLL

(103)

where a1 is a constant (= 0.57), a2 = 0.62 a1, t is the plume travel time (= x/U), and TLr is a 

Lagrangian time scale for relative dispersion defined as

TTLL = aarr1

zz

σσww

(104)

where ar1 = 0.46, zi is the mixing height, and σw is the vertical turbulence parameter. The 

turbulence dissipation rate, γ, is calculated as follows, based on Weil (1996):

εε = bbσσww

2

TTLL

(105)

where b is a constant (= 0.78). The values of wind speed (U) and σw used in eqs. (103) through 

(105) are the effective values, calculated as averages across the layer from the plume centroid 

height to the receptor height (up to 2.15Φz), following the procedure used in AERMOD to 

calculate effective values. Using the effective values of σw, AERMOD calculates effective 

values of the turbulence dissipation rate, γ. 

Since the relative dispersion coefficients are source- and meteorology-dependent in 

AERMOD, the model generates a table of relative dispersion coefficients as a function of 

distance for the dominant source for each receptor and each hour in order to complete the plume 

volume calculation.
82

The original PVMRM utilized the relative dispersion coefficients described above to 

define the plume volume. These relative dispersion coefficients are applicable to unstable/

convective conditions, but are likely to overpredict the plume volume for stable conditions, 

resulting in overpredictions of NO2 concentrations. The PVMRM algorithm was modified for 

version 15181 to use the “standard” total dispersion coefficients incorporated in AERMOD to 

define the plume volume during stable conditions.

5.8.1.3 Treatment of volume and area sources

If the dominant source is a volume source, then the initial lateral and vertical dimensions 

of the volume source are included in the calculation of the relative dispersion coefficients for 

purposes of calculating the plume volume, as follows:

σσrr = (σσrr

3 + σσrr

3 + σσ0)1/3 (106)

where σ0 is the initial dispersion coefficient of the volume source calculated as �σσyy0σσzz0 based 

on the initial lateral (Φy0) and vertical (Φz0) dimensions input by the user. If a volume source is 

included among the major contributing sources it is treated the same as a point source in defining 

the combined plume volume.

For application of PVMRM to area sources, the plume volume is extended laterally if 

necessary to include the projected width of the area source or sources that are included among 

the major contributing sources. The emissions from an area source are included in the 

calculation of the NOx moles emitted into the plume if the centroid of the area source is within 

the box defined by the alongwind and crosswind extent of major contributing sources. In 

addition, if an area source is the dominant source, then the relative dispersion coefficients are 

calculated based on the radial distance from the centroid of the area source to the receptor.
83

5.8.1.4 Defining extent of plume

Since relative dispersion coefficients are used to define the plume volume, the number of 

standard deviations from the plume centerline, nz, used in the calculation of plume volume was 

increased from the value used by Hanrahan (1999) for ISCST3. The ISCST3 postprocessor 

version used a value of 1.282 for nz, corresponding to 80 percent of the area under the normal 

curve. The plume volume calculations for AERMOD are based on a value of nz = 4.0, which 

corresponds to about 99.99 percent of the area under the normal curve. The minimum value of 

the dispersion coefficient was also reduced from the 15m minimum used with ISCST3 to a 

minimum of 5m for AERMOD in order to maintain approximately the same minimum plume 

volume in AERMOD as used for ISCST3. A minimum value of 4.8m in AERMOD would 

provide the same minimum plume volume as used by ISCST3 with nz = 1.282 and a minimum 

dispersion coefficient of 15m.

The original implementation of PVMRM used four (4) times the relative dispersion 

coefficients to define the plume volume, which accounts for more than 99.99% of the plume. 

Given the fact that the PVMRM option in AERMOD assumes full and instantaneous mixing of 

the NO and O3 within the plume, using such a large percentage of the plume volume may 

introduce a bias to overpredict ambient concentrations of NO2. Beginning with version 15181, 

the PVMRM option uses 2.58 times the relative dispersion coefficients to define the plume 

volume for unstable conditions, which accounts for about 99% of the plume. For stable 

conditions, the PVMRM option uses 1.282 times the total dispersion coefficients to define the 

plume volume, consistent with the original approach proposed by Hanrahan (1999), which 

accounts for about 80 percent of the plume volume. However, since AERMOD incorporates a 

horizontal meander algorithm that increases lateral plume spread beyond that accounted for 

based on dispersion coefficients, the number of sigmas used to define the plume volume for 

stable conditions is adjusted to account for meander, i.e., 

NN NN NN = min(2.15,1,282 ∗ (SSSS /SSSS)

where SYEFF is the effective σy value that replicates the plume centerline associated with 

meander, but based on a standard Gaussian plume calculation.
84

The original implementation of PVMRM used the radial distance from source to 

receptor to calculate the plume volume and the moles of NOx contained in the plume. 

Beginning with version 15181, the downwind distance is used to calculate these values. Use of 

the downwind distance provides a more realistic estimate of NOx conversion consistent with a 

straight-line, steady-state plume model, such as AERMOD. 

5.8.1.5 Adaption for AERMOD terrain algorithm

The vertical dimension of the plume volume is based on the relative dispersion 

coefficient for the dominant source and the range in plume heights for the major contributing 

sources. Since the effective plume heights differ for the terrain following and terrain responding 

components, the vertical dimension was modified to calculate the range of plume heights 

separately for both the terrain following and terrain responding components, and then use a 

weighted value for the vertical dimension based on the terrain (plume state) weighting factor, f, 

defined in Section 0.

5.8.1.6 Treatment of penetrated plumes

For unstable conditions with partial or full plume penetration above the mixing height, 

zi, separate relative dispersion coefficients are calculated for the penetrated portion of the 

dominant plume. For cases with partial penetration for the dominant plume, AERMOD 

calculates two plume volumes, one based on relative dispersion coefficients for the direct source 

and another based on the relative dispersion coefficients for the penetrated source. Since 

AERMOD uses the same dispersion coefficients for the direct and indirect sources, separate 

values of relative dispersion coefficients for the indirect source are not needed. The effective 

plume volume used in the application of PVMRM is based on a weighted average of the direct 

and penetrated plume volumes using the plume penetration factor (PPF) for the dominant 

source. The model stores the plume centroid heights for both the direct and penetrated plumes

for all sources at each receptor, and these are used to incorporate the effect of the major 

contributing sources on the volumes for the direct and penetrated plumes.
85

5.8.2 Minimum ozone concentration for stable conditions

Conditions for implementation with ISCST3, due to the fact that surface measurements 

may be artificially low during nighttime stable conditions due to the formation of a stable 

vertical temperature gradient. Since the AERMOD model does not use Pasquill-Gifford (P-G) 

stability categories, this minimum ozone concentration was modified to use Monin-Obukhov 

length as the stability parameter. The AERMOD model first keeps track of the maximum ozone 

concentration over the previous 24 hours. If the Monin-Obukhov length is positive (i.e. stable), 

with a value of less than 50 meters (very stable), then the maximum ozone concentration over 

the previous 24 hours is used as the minimum value. If the Monin-Obukhov length is positive 

and the value is over 500 meters (nearly neutral), then no minimum ozone concentration is 

applied for that hour. If the Monin-Obukhov length is between 50 meter and 500 meters, then 

the minimum ozone concentration is determined by linear interpolation, i.e., the minimum value 

is calculated as O3MAX * (500 – L)/450, where O3MAX is the maximum ozone concentration 

over the previous 24 hours, and L is the Monin-Obukhov length in meters.

5.9 Adjustments for the urban boundary layer

Although urban surface characteristics (roughness, albedo, etc.) influence the boundary 

layer parameters at all times, the effects of the urban sublayer on the structure of the boundary 

layer is largest at night and relatively absent during the day (Oke 1998). An urban “convectivelike” boundary layer forms during nighttime hours when stable rural air flows onto a warmer 

urban surface. Following sunset, the urban surface cools at a slower rate than the rural surface 

because buildings in the urban area trap the outgoing thermal radiation and the urban subsurface 

has a larger thermal capacity. AERMOD accounts for this by enhancing the turbulence above 

that found in the rural stable boundary layer (i.e., a convective-like urban contribution to the 

total turbulence in the urban SBL). The convective contribution is a function of the convective 

velocity scale, which in turn, depends on the surface heat flux and the urban mixed layer height. 

The upward heat flux is a function of the urban-rural temperature difference.

The urban-rural temperature difference depends on a large number of factors that cannot 

easily be included in applied models such as AERMOD. For simplicity, the data presented in 
86

Oke (1973; 1982) is used to construct an empirical model. Oke presents observed urban-rural 

temperature differences for a number of Canadian cities with populations varying from about 

1000 up to 2,000,000. These data represent the maximum urban effect for each city since they 

were collected during ideal conditions of clear skies, low winds, and low humidities. An 

empirical fit to the data yields the following relationship

(107)

where ΔTmax = 12 ̊C, Po = 2,000,000 (the city population associated with the maximum 

temperature difference in Oke’s data), and P is the population of the urban area being modeled.

Since the ambient nighttime temperature of an urban area is higher than its surrounding 

rural area, an upward surface heat flux must exist in the urban area. It is assumed that this 

upward surface heat flux is related to the urban-rural temperature difference through the 

following relationship

(108)

where α is an empirical constant, ρ is the density of air, and cp is the specific heat at constant 

pressure. This expression is analogous to the bulk transfer parameterization of heat flux over a 

homogeneous surface (e.g., Businger (1973)), with α defined as the “bulk” transfer coefficient. 

We chose α to ensure that the upward heat flux is consistent with maximum measured values of 

the order of 0.1 m s-1 o

C . Because ΔTu-r has a maximum value on the order of 10 oC, and u* is 

on the order of 0.1 m s-1, α should have a maximum value on the order of 0.1. Although we 

assume that α has a maximum (city center) value of about 0.1, AERMOD uses an effective value 

of α that is averaged over the entire urban area. Assuming a linear variation of α from 0 at the 

edge of the urban area to about 0.1 at the center of the urban area results in an areal average 

equal to one-third of that at the center (since the volume of cone is one-third of that of a right 

circular cylinder of the same height). Therefore, AERMIC tested an area-averaged value of α

equal to 0.03 against the Indianapolis data. This choice for α is consistent with measured values 

∆T ∆T P u r Po − = 



 



 + 



 



 max 0.1 ln 1.0 ,

H c T u u = α ρ p ∆ u−r *,
87

of the upward heat flux in Canadian cities reported by Oke (1973; 1982). The results of the 

developmental testing indicated that this choice for α resulted in an adequate fit between 

observations and AERMOD-predicted concentrations.

The mixing height in the nighttime urban boundary layer, ziu, is based on empirical 

evidence presented in Oke (1973; 1982) that, in turn, suggests the following relationships:

(109)

where R is a measure of the city size and P is the population of the city. The first relationship is 

based on the observed growth of the internal convective boundary layer next to shorelines 

(Venkatram 1978). The second relationship implicitly assumes that population densities do not 

vary substantially from city to city.

Equation (103) leads to the following equation for the nocturnal urban boundary layer 

height due to convective effects alone:

(110)

where ziuo is the boundary layer height corresponding to Po. Based on lidar measurements taken 

in Indianapolis (1991), and estimates of ziu found by Bornstein (1968) in a study conducted in 

New York city, ziuo is set to 400 m in AERMOD.

In addition, since effects from urban heating should not cause ziu to be less than the 

mechanical mixing height, ziu is restricted from being less than zim. Therefore, the mixed layer 

height for the nighttime urban boundary layer is computed as:

(111)

ziu ≈ R and R ≈ P 1 2 1 2

,

ziuc = ziuo (P Po )

1/ 4

ziu = Max[ziuc ; zim].
88

Once the urban mixing height has been estimated, a surrogate convective velocity scale 

(appropriate for the magnitude of convective turbulence present) is computed by substituting ziu

and Hu into the definitional equation for w* (Deardorff 1970). That is, 

(112)

where w*u is the urban nighttime convective velocity scale and T is the near-surface air 

temperature.

Having estimated w*u the turbulence in the nighttime urban can be enhanced using the 

expressions found in Section 0. However, since for low level sources σwT depends primarily on 

u* (see eqs. (34) and (35)) it is not possible to directly enhance σwT for these sources using w*u. 

Therefore, an effective friction velocity (u*eff) is developed as a surrogate for w*u in the lower 

portion of urban PBL. We define u*eff as the friction velocity that is consistent with σwm = σwc at 

z = 7zo. Assuming that z = 7zo is always less 0.1ziuc, u*eff is estimated by equating σwc (eq. (35)) 

with σwm (eq. (37)) and solving for u*. Once u*eff is found, the urban friction velocity for 

nighttime conditions (u*u) is calculated as the maximum of u*eff and u* (the rural and daytime 

urban friction velocity).

Then using the enhanced velocity scales u*u and w*u, the nighttime convective portion of 

the turbulence in the urban boundary layer is computed using the expressions turbulence found 

in Section 0. That is, σwc and σwm are calculated from eqs. (35) and (37) , respectively, with w*u

used in place of the daytime convective velocity scale (w*) and u*u substituted for the rural u*. 

Furthermore, for consistency purposes, a urban nighttime Monin-Obukhov length is calculated 

using eq. (8) with substitutions u*u for u* and Hu (eq. (108)) for H.

Finally, the total nighttime turbulence in the urban boundary layer is calculated as the 

sum (in quadrature) of the convective and mechanical portions. With these enhanced levels, 

vertical dispersion due to ambient turbulence (σza) in the urban boundary layer is calculated 

from eq. (83) (the SBL formulation for σza ) with the urban PBL assumed to be neutral (i.e., N = 

w g H z

c T u u iuc

p

*

/

=















 ρ

1 3
89

0). For the lateral dispersion in the urban boundary layer, σya is calculated using the SBL 

formulation given by eq.(76). 

The potential temperature gradient in the night-time urban boundary layer is set equal to 

the upwind rural profile (Section 0) for all heights above ziu, and is assumed to be equal to a 

small positive value below ziu; i.e.,

(113)

For plumes below ziu , the effective reflection surface is set equal to the height of the 

urban boundary layer (i.e., zieff = ziu). Plumes that rise above ziu (hes > ziu) are modeled with a zieff

that is calculated from eq. (68) with zim replaced by ziu. Plume rise in the urban stable boundary 

layer is calculated from eqs. (95)-(99) with ∂θ/∂z taken from eq.(113).

Use of this value for ∂θ/∂z provides an appropriate near-neutral plume rise formulation 

that is expected within the nocturnal urban boundary layer. However, plume height in these 

conditions is not allowed to exceed 1.25 ziu.

For daytime conditions (L < 0) in urban areas, AERMOD uses the same formulations as 

in rural areas (i.e., no urban-related adjustments to boundary layer characteristics).

∂θ ∂

∂θ ∂

z for z z

z rural value for z z

iu

iu

= ≤

= >

− 10 5

.
90

6. List of symbols

Bo Bowen ratio - ratio of the sensible to latent heat fluxes (dimensionless)

CAERMOD concentration estimated using AERMOD without considering building 

wake effects (g m-3)

Cc,s{xr,yr,zr} concentration contribution from the horizontal plume state - convective 

and stable (g m-3)

Cc,s{xr,yr,zp} concentration contribution from the terrain-following plume state -

convective and stable (g m-3)

Cc{xr,yr,zr} total concentration (CBL) (g m-3

)

Cd{xr,yr,zr} concentration contribution from the direct source (CBL) (g m-3)

Cp{xr,yr,zr} concentration contribution from the penetrated source (CBL) (g m-3

)

Cr{xr,yr,zr} concentration contribution from the indirect source (CBL) (g m-3

)

Cs{xr,yr,zr} total concentration (SBL) (g m-3

)

CT{xr,yr,zr} total concentration (CBL) (g m-3)

Cch concentration from the coherent plume used in meander calculations (gm-3

)

CR concentration from the random plume used in meander calculations (g m-3)

CD neutral drag coefficient (cal g-1 °C-1)

Cprime concentration estimated using the PRIME algorithms with AERMODderived meteorological inputs (g m-3

)

cp specific heat at constant pressure (= 1004 J g-1 K-1)

Fb plume buoyancy flux (m4 s3)

Fy total horizontal distribution function - with meander (m-1)

FyC horizontal distribution function for a coherent plume (m-1)

FyR horizontal distribution function for a random plume (m-1)

FG flux of heat into the ground (W m-2)

Fm plume momentum flux (m4s2) 

Fz total vertical distribution function (m-1

)

f plume state weighting function (dimensionless)

fp fraction of plume mass contained in CBL = (1 - penetration factor) 

dimensionless)

g acceleration due to gravity (9.8 m s-2)

H sensible heat flux (W m-2

)

Hc critical dividing streamline (m)

Hp plume centroid height (m)

Hu heat flux in the nighttime boundary layer (W m-2

)

hc receptor specific terrain height scale (m)

hep penetrated source plume height above stack base (m)

hs stack height corrected for stack tip downwash (m)

Δh general symbol for distance dependent plume rise (m)
91

Δhd plume rise for the direct source (m)

Δheq equilibrium plume rise in a stable environment (m)

Δhh depth of the layer between zi and the stack top (m)

Δhp plume rise for the penetrated source (m)

Δhi plume rise for the indirect source (m)

Δhs plume rise for the stable source (m)

iz vertical turbulence intensity

k von Karman constant k = 0.4 (dimensionless)

l length scale used in determining the Lagrangian time scale (m)

ln neutral length scale - a component of l (m)

ls stable length scale - a component of l (m)

L Monin-Obukhov length (m)

N Brunt-Vaisala frequency (s-1

)

Nh Brunt-Vaisala frequency above zi (s-1

)

n cloud cover (fractional)

P population of urban area

py lateral probability density function

pz vertical probability density function

pw probability density function of the instantaneous vertical velocities

Q source emission rate (g/s)

R solar insolation (W m-2

)

Rn net radiation (W m-2

)

Ro clear sky solar insolation (W m-2)

r(φ) Albedo {solar elevation} (dimensionless)

rʹ noontime albedo (dimensionless)

rs stack radius - corrected for stack tip downwash (m)

ry lateral dimension of an elliptical plume

rz vertical dimension of an elliptical plume

S skewness factor (dimensionless)

T ambient temperature (K)

TLy lateral lagrangian time scale (sec)

TLzc vertical lagrangian time scale for the CBL (sec)

TLzs vertical lagrangian time scale for the SBL (sec)

Tr time scale used in the meander algorithm (sec)

Tref ambient temperature - at reference temperature height (K)

Ts stack gas temperature (K)

Tu urban surface temperature (K)

t time (sec)

ΔT difference between stack gas and ambient temperature (K)

ΔTu-r urban-rural temperature difference (K)
92

u wind speed (m s-1

)

ucr minimum speed for which the expression for u*, in the SBL, has a real 

valued solution (m s-1)

uo defined in eq. (14) and used in eq. (15)

up wind speed that is used for plume rise (m s-1)

uref wind speed at reference height (m s-1

)

uth wind speed instrument threshold - separate value for each data set (offsite 

& onsite) (m s-1

)

u* surface friction velocity (m s-1)

u*eff effective surface friction velocity - surrogate for w*u (m s-1

)

u*u surface friction velocity for nighttime urban conditions (m s-1)

w random vertical velocity in the CBL (m s-1)

wwjj mean vertical velocity for the updraft (j = 1) and the downdraft (j = 2) 

distributions (m-s-1

)

ws stack exit gas velocity (m-s-1)

w* convective velocity scale (m-s-1)

w*u urban nighttime convective velocity scale (m-s-1)

X non-dimensional downwind distance (dimensionless)

xr downwind distance to a receptor (m)

xf distance to final plume rise (m) - eq. (44) for the CBL; eq. (96) for the 

SBL

xm downwind distance at which plume material uniformly mixed throughout 

the boundary layer (m)

(xr,yr,zr) receptor location

(xt,yt,zt) terrain point location

zbase user specified elevation for the base of the temperature profile (i.e., 

meteorological tower)

zc total height of the plume in the CBL considering both plume rise and 

effects from convective turbulence (m)

zi mixing height (m): zi = MAX [zic; zim] in the CBL and zi = zim in the SBL

zic convective mixing height (m)

zie equilibrium height of stable boundary layer

zieff height of the reflecting surface in the SBL or in the stable layer above the 

above the CBL (m)

zim mechanical mixing height (m)

ziu urban nighttime boundary layer mixing height (m)

ziuc urban nighttime boundary layer mixing height due to convective effects 

alone (m)

zmsl height of stack base above mean sea level (m)

zo surface roughness length (m)
93

zPG release height used in the Prairie Grass experiment (m)

zp receptor “flagpole” height - height of a receptor above local terrain (m)

zr height of the receptor above local source base (m)

zref reference height for wind (m)

zTref reference height for temperature (m)

zt height of the terrain above mean sea level (m)

αα� general symbol used to represent the effective parameters in the treatment 

of the inhomogeneous boundary layer. In the text the effective values of 

the parameters u, σw, σv and TL are denoted by underscoring the character

γ parameter used to weight CAERMOD and CPrime in estimating concentrations 

that are influenced by building downwash (dimensionless)

θ potential temperature (K)

θ* temperature scale (K)

λj weighting coefficient for the updraft (j = 1) and downdraft (j = 2) 

distributions of eqs. (53),(59) and (65)

ρ density of air (kg m-3

)

σb buoyancy induced dispersion for the direct & indirect sources (m)

σh

2 total horizontal wind “energy” used in the meander algorithm (m2

)

σr

2 random “energy” component of the total horizontal wind “energy” used in 

the meander algorithm (m2)

σSB Stephen Boltzman constant (5.67x10-8 Wm-2K-4)

σu along-wind turbulence (m s-1

)

σv lateral turbulence (m s-1

)

σvc convective portion of the lateral turbulence (m s-1)

σvo surface value of the lateral turbulence (m s-1)

σvm mechanical portion of the lateral turbulence (m s-1)

σvT total lateral turbulence (m s-1

)

σw vertical turbulence (m s-1)

σwc convective portion of the vertical turbulence (m s-1)

σwm mechanical portion of the vertical turbulence (m s-1)

σwml mechanical portion of the vertical turbulence generated in the PBL (m s-1

)

σwmr mechanical portion of the vertical turbulence above the PBL (residual) (m s-1)

σwT total vertical turbulence (m s-1)

σxg longitudinal dimension of the building wake (m)

σy total lateral dispersion for the direct & indirect sources (m)

σya,zaj ambient turbulence induced dispersion for the direct & indirect sources (m)

σzas ambient dispersion for the stable source (m)

σyg distance from the building centerline to lateral edge of the building wake (m)

σyl lateral spread from combined effects of ambient turbulence and building 

downwash (m)
94

σzp total dispersion for the penetrated source (m)

σzs total dispersion for the stable source (m)

σzaj ambient vertical dispersion for the updraft & downdrafts plumes (j = 1,2), 

respectively, for both the direct & indirect sources (m)

σzej elevated portion of σzaj (m)

σzes elevated portion of σzas (m)

σzg height of the building wake at the receptor location (m)

σzj total vertical dispersion for the updrafts and downdrafts (j = 1,2 respectively), 

for both the direct and indirect sources

σzg surface portion of σzaj (m)

σzgs surface portion of σzas (m)

τ time constant controlling the temporal interpolation of zim (sec)

φ solar elevation angle

φp fraction of plume mass below Hc (dimensionless)

Ψdj total height of the direct source plume (i.e. release height + buoyancy + 

convection) (m)

Ψrj total height of the indirect source plume (m)

ψm similarity function for momentum (stability correction) - eq. (7) for the 

CBL and eq. (29) for the SBL (dimensionless)
95

Appendix A: Input / output needs and data usage

A.1 AERMET input data needs

Besides defining surface characteristics, the user provides several files of hourly 

meteorological data for processing by AERMET. At the present time AERMET is designed to 

accept data from any for the following sources: 1) standard hourly National Weather Service 

(NWS) data from the most representative site; 2) morning soundings of winds, temperature, and 

dew point from the nearest NWS upper air station; and 3) site-specific wind, temperature, 

turbulence, pressure, and radiation measurements (if available). 

The minimum measured and/or derived data needed to run the AERMOD modeling 

system are as follows:

A.1.1 Meteorology

• wind speed (u); 

• wind direction; 

• cloud cover - opaque first then total (n);

• ambient temperature (t); 

• morning sounding.

Cloud cover is also used in dry deposition calculations in the AERMOD model. 

Therefore, if cloud cover is missing and the Bulk Richardson Number Scheme is being used (see 

3.3.1) then an equivalent could cover is calculated as follows, based on van Ulden and Holtslag 

(van Ulden and Holtslag 1985):

(114)

where θ* is the temperature scale as calculated from eq. (18).

neq =  −



 



 1 0 09

0 5

0 5 θ*

. .

.
96

A.1.2 Directionally and/or Monthly Varying Surface Characteristics

• noon time albedo (rʹ );

• Bowen ratio (Bo);

• roughness length (zo) -

For AERMET, the user can specify monthly variations of three surface characteristics for 

up to 12 upwind direction sectors. These include: the albedo (r), which is the fraction of 

radiation reflected by the surface; the Bowen ratio (Bo), which is the ratio of the sensible heat 

flux to the evaporation heat flux; and the surface roughness length (zo) , which is the height 

above the ground at which the horizontal wind velocity is typically zero. 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. The information presented in the user’s guide is not be 

considered regulatory guidance. The user is encouraged to research the literature to determine 

the most appropriate values for surface characteristics, for a specific application.

A.1.3 Other

• Latitude;

• longitude;

• time zone; 

• wind speed instrument threshold for each data set ( uth ).

A.1.4 Optional

• solar radiation;

• net radiation (Rn); 

• profile of vertical turbulence (σw); 

• profile of lateral turbulence (σv)
97

A.2 Selection and use of measured winds, temperature, and turbulence in AERMET

A.2.1 Threshold Wind Speed

The user is required to define a threshold wind speed (uth) for site-specific data sets. 

Although the current version of AERMOD cannot accept a separate uth for NWS data, a separate 

uth should be selected for each on-site data set being used.

A.2.2 Reference Temperature and Height

The reference height for temperature (zTref), and thus the reference temperature, is 

selected as the lowest level of data which is available between zo & 100 m.

A.2.3 Reference Wind Speed and Height

The reference height for winds (zref), and thus the reference wind speed (uref), is selected 

as the lowest level of data which is available between 7 zo & 100m. Although the current 

version of AERMOD cannot accept a separate zref for offsite data, we believe that a separate zref

should be selected for each data set being used.

If no valid observation of the reference wind speed or direction exists between these 

limits the hour is considered missing and a message is written to the AERMET message file. 

For the wind speed to be valid its value must be greater than or equal to the threshold wind 

speed. AERMOD processes hours of invalid wind speed, e.g. calms, in the same manner as ISC 

(EPA calms policy).

All observed wind speeds in a measured profile that are less than uth are set to missing 

and are therefore not used in the construction of the wind speed profile (profiling of winds is 

accomplished in AERMOD).
98

A.2.4 Calculating the Potential Temperature Gradient above the Mixing Height from Sounding 

Data

AERMET calculates dθ/dz for the layer above zi as follows:

• If the sounding extends at least 500 m above zi the first 500 m above zi is used to 

determine dθ/dz above zi.

• If the sounding extends at least 250 m above zi (but not 500 m) then the available 

sounding above zi is used to determine dθ/dz above zi.

• AERMET limits dθ/dz above zi to a minimum of 0.005 K m-1.

• If the sounding extends less than 250 m above zi then set dθ/dz = 0.005 K m-1 (a default 

value).

A.2.5 Measured Turbulence

All measured turbulence values are passed to AERMOD if the hour is non-missing. This 

is true even for those levels where the wind speed is below uth. Based on measurements with 

research grade instruments, reasonable minimum turbulence levels in non-calm conditions for 

vertical turbulence (σw) and lateral turbulence (σv) values are set by AERMOD to 0.02 m s-1 and

0.2 m s-1, respectively. Although these lower limits are applied to the measured values of the 

turbulence the calculated profile values of σw and σv are not subjected to any lower limits. We 

do not restrict these estimated profiles because it would bias the calculation of the effective 

values of turbulence, which are averages through the layer between the receptor and the plume 

height, in determining the dispersion of the plume. However, as discussed in Section 0, these 

limits are applied to the effective values of turbulence and wind speed.

A.2.6 Data Substitution for Missing On-Site Data

If on-site data are missing for an hour, the hour is considered missing unless the user 

specifies a substitute data set. AERMET does not default to NWS (or any other offsite) data.
99

A.3 Information passed by AERMET to AERMOD

The following information is passed from AERMET to AERMOD for each hour of the 

meteorological data record.

• All observations of wind speed (u); wind direction; ambient temperature (T); lateral 

turbulence (σv); & vertical turbulence (σw) with their associated measurement heights.

• Sensible heat flux (H), friction velocity (u*), Monin Obukhov length L, zim (for all hours), 

zic & w* (for convective hours only), zo , r{φ}, & Bo, dθ/dz (above zi), uref , wind direction 

at the reference height, zref , ambient temperature at the reference height (Tref) (not used 

in AERMOD), & the reference height for temperature (zTref)

A.4 Restrictions on the growth of the PBL height

AERMET restricts the growth of zi to a reasonable maximum of 4000 m. This 

restriction applies to both calculated and measured mixing heights. Although mixing heights in 

excess of 4000 m may occur on rare occasions, in desert climates, the additional effect on 

surface concentration is most likely insignificant.

A.5 Initializing the mechanical mixing height smoothing procedure

If {t + Δt}, in eq. (26), is the first hour of the data set then no smoothing takes place. 

Furthermore, if a missing value occurs at time step t then smoothing is not performed at time 

step {t + Δt} but is restarted for subsequent hours.

A.6 Determining the mixing height when the sounding is too shallow

The left hand side of eq. (22) is determined from the morning temperature sounding and 

the right hand side from the daytime history of surface heat flux. When the temperature 

sounding, obtained from the NWS, does not reach a height which is greater than the convective 

mixing height, we must assume a profile for the potential temperature gradient in order to 

estimate zic. This is accomplished as follows:
100

• Determine dθ/dz in the top 500 m layer of the sounding. However, if part of the 500 m 

layer is within the first 100 m of the PBL, the layer should be reduced (to a minimum 

thickness of 250 m) to avoid using the portion of the sounding that is below 100 m. If the 

above conditions cannot be satisfied then zic is defined as missing.

• Extend the sounding by persisting dθ/dz up and recomputing zic.

• Provide warning messages which tell users

- the height of the actual sounding top,

- that dθ/dz has been extrapolated above the sounding zic, and

- that zic has been recomputed.

• Allow the user to reject the “fixed-up” value for zic by defining it as missing.

A.7 Input data needs for AERMAP

The following data is required input for AERMAP

• DEM formatted terrain data ( xt, yt, zt ),

• Design of receptor grid; AERMAP accepts either polar, Cartesian or discrete receptors.

A.8 Information passed by AERMAP to AERMOD

AERMAP passes the following parameters to AERMOD: xr, yr, zr, zt, & the height scale 

(hc ) for each receptor.

A.9 Wind speed and turbulence limits used in model calculations

When calculating the effective parameters, limits are placed on the such that:

(115)

These limits are also applied when selecting the turbulence for plume rise calculations.

[ ]

[ ]

σ σ

σ σ

w w

v v s

z Max z ms

z Max z u h ms

{ } { }; .

{ } { }; . { }; . .

=

=

−

−

0 02

0 05 0 2

1

1
101

Dilution of the plume is determined by the wind that corresponds to the average over the 

magnitudes of the wind vectors during a given time interval. But measurements only give the 

vector averaged wind, which can be zero, even though the dilution wind is not zero. We can 

estimate the dilution wind by assuming that the vector wind, uv, can be expressed as

(116)

where uu is the mean measured wind, and the primed quantities refer to the turbulent fluctuations. 

The assumption being made is uuvv = uu. If we assume that the measured velocity fluctuations 

correspond only to the angular variations of a constant vector, uv, we can write from eq. (116)

that

(117)

In this simple model, uv , is the dilution wind. If we take σu = σv, the dilution wind can 

be written as

(118)

This formulation assures that the dilution wind is not zero as long as either uu or σv is not 

zero. Similarly, at the time of plume rise calculations, the effective turbulence and effective 

wind speed will be recalculated using eqs. (115) and (118), where the turbulence and winds will 

be evaluated at stack top.

A.10 Using profiles for interpolating between observations 

When observations are available AERMOD uses the similarity profile functions to 

interpolate adjacent measurements. Figure 17 illustrates how AERMOD’s INTERFACE uses 

the expected shape of a meteorological profile to interpolate between observations.

u u u v v = ( + ′, ′)

u u v v u

2 2 2 2 = +σ +σ .

u h u h s s v { } = { } + . 2 2 2σ
102

Figure 17. AERMOD’s construction of a continuous meteorological profile by interpolating 

between observations.

For a gridded profile height between two observed profile heights, the observations are 

interpolated to the gridded height while maintaining the shape of the similarity profile. This is 

accomplished as follows:

1. the observations are linearly interpolated to the gridded profile height;

2. the similarity function is evaluated at the gridded profile height;

3. the similarity function is evaluated at the observed profile heights immediately 

above and below the grid height and linearly interpolated to the grid height;
103

4. the ratio of the value obtained in 2 to the value obtained in 3 is applied to the 

value obtained in 1.

For a gridded profile height above the highest observation, the procedure is modified 

slightly:

1. the observation at the highest observed profile height is extrapolated by 

persisting the value upward;

2. the similarity function is evaluated at the grid height;

3. the similarity function is evaluated at the highest height in the observed profile;

4. the ratio of the value obtained in 2 to the value obtained in 3 is applied to the 

value obtained in 1.

A similar procedure for extrapolating to heights above the observed profile is applied to 

heights below the lowest observed profile height.

A.11 Using measured mixing heights

If measured mixing heights are available, then they are treated in the following manner: 

If L>0 (SBL) the measured mixing height is defined as zie and it is treated the same as a 

calculated mechanical mixing height (smoothed as explained in Section 0). If L<0 (CBL) the 

measured mixing height is defined as zic, and zie is calculated from eq. (24), smoothed, then 

proceed as if both zic and the smoothed zim had been calculated values.

If a user has “measured” mixing heights available (and chooses to use them), AERMET 

defaults to substituting calculated mixing heights for missing measurements and a message is 

written that a substitution has occurred. If the user elects to substitute calculations for missing 

measurements, AERMET will print out a message to the message file for each hour that a 

substitution has occurred.
104

Appendix B: Model evaluation results

B.1 Introduction

This evaluation presents a benchmark of model performance based on the original field 

studies presented in Cimorelli, et al, 2005 and Perry, et al, 2005. The evaluation focused on the 

benchmarking the performance of the final (16216) versions of the AERMOD modeling system 

associated with the 2016 update to the Guideline on Air Quality Models and the current update 

(18081). The statistical analysis determines the best performing version of the model for 14 of 

the original 17 databases, including the adjust u* option formally adopted as a regulatory option 

in the final version of AERMOD (16216r). 

B.2 Database descriptions

The 14 databases used in this evaluation are briefly described in this section and 

summarized in Table 1. The stack heights, terrain complexity, urban/rural status, importance of 

downwash, inclusion of turbulence parameters and meteorological data included for the 

database are listed for each area. A more complete description of these databases can be found 

in U.S. EPA, 2003. The databases are arranged by the following hierarchy: Two categories of 

turbulence inclusion (inclusion of turbulence or no turbulence). Within each of those categories, 

databases were ordered by complexity of terrain (complex or flat), and within those two 

categories, databases were ordered by increasing height.
105

Table 1. AERMOD evaluation databases used for comparisons of AERMOD 15181 and 

AERMOD 16216. Databases in gray are also subject to the EPA’s protocol for determining best 

performing model.

Location Stack 

heights

Urban

/rural

Terrain Downwash Turbulence 

parameters

Site specific 

AERMET inputs

Martins 

Creek

59, 76, 183 

m

Rural Complex Yes 10 m σv, 

σw

10m wind, 

temperature; 90-

420 m wind (every 

30 m).

Tracy 91 m Rural Complex No σv, σw 10 and 50-400 m 

(every 25 m) 

wind, temperature

Lovett 145 m Rural Complex No σv, σw 10, 50, and 100 m 

wind, temperature

Westvaco 190 m Rural Complex No σv, σw 30, 210, 326, 366, 

and 416 m wind, 

temperature1

DAEC 1 m, 24 m, 

46 m

Rural Flat Yes σv Insolation, 10, 23.5 

and 50 m wind, 

temperature

EOCR 1, 25, 30 m Rural Flat Yes σv 4, 10, and 30 m 

wind, temperature

Alaska 39.2 m Rural Flat Yes σv, σw 33 m wind, 

temperature

Indianapolis 84 m Urban Flat No σv, σw Station pressure, 

net radiation, 10 m 

wind, temperature

Kincaid 187 m Rural Flat No σv, σw Net radiation, 

insolation, 10 , 30, 

and 50 m wind, 

temperature

AGA 9.8, 14.5, 

24.4 m

Rural Flat Yes None 10 m wind and 

temperature

Millston 3 stacks 29 

m (freon) 

48 m (SF6)

Rural Flat Yes None 10 m wind speed ; 

43.3 m wind and 

temperature

Bowline 2 stacks 

86.87 m

Rural Flat Yes None 100 m winds and 

temperature

Baldwin 3 stacks 

184.4 m 

Rural Flat Yes None2 10 and 100 m 

wind, temperature

Clifty 

Creek

3 stacks 

207.9 m

Rural Flat/Elev No None 10 m temperature; 

60 m wind

• 1. 30 m observations removed from AERMOD profile before running AERMOD.
106

B.2.1 Martin’s Creek

The Martins Creek Steam Electric Station is located in a rural area along the Delaware 

River on the Pennsylvania/New Jersey border, approximately 30 km northeast of Allentown, PA 

and 95 km north of Philadelphia, PA (Figure B-1). The area is characterized by complex terrain 

rising above the stacks. Sources included multiple tall stacks ranging from 59 to 183 m in height, 

including Martins Creek and three background sources located between 5 and 10 km from 

Martins Creek. The seven SO2 monitors were located on Scotts Mountain, which is about 2.5 - 8 

km southeast of the Martins Creek facility. On-site meteorological data covered the period from 

1 May 1992 through 19 May 1993. Hourly temperature, wind speed, wind direction, and sigmatheta (standard deviation of the horizontal wind direction) at 10 m were recorded from an 

instrumented tower located in a flat area approximately 2.5 km west of the plant. In addition, 

hourly multi-level wind measurements were taken by a sodar located approximately three km 

southwest of the Martins Creek station.
107

Figure B-1. Martin's Creek study area.
108

B.2.2 Tracy

The Tracy Power Plant is located 27 km east of Reno, Nevada in the rural Truckee River 

valley completely surrounded by mountainous terrain (Figure B-2). A field tracer study was 

conducted at the power plant in August 1984 with SF6 being released with the moderately 

buoyant plume from a 91-m stack. A total of 128 hours of data were collected over 14 

experimental periods. Stable atmospheric conditions were dominant for this study. Site-specific

meteorological data (wind, temperature, and turbulence) for Tracy were collected from an 

instrumented 150-m tower located 1.2 km east of the power plant. The wind measurements 

from the tower were extended above 150 meters using a Doppler acoustic sounder and 

temperature measurements were extended with a tethersonde.
109

Figure B-2. Tracy power plant study area.

B.2.3 Lovett

The Lovett Power Plant study consisted of a buoyant, continuous release of SO2 from a 

145-m tall stack located in a complex terrain, rural area in New York State (Figure B-3). The 

data spanned one year from December 1987 through December 1988. Data were collected from 

12 monitoring sites (ten on elevated terrain and two near stack-base elevation) that were located 

about 2 to 3 km from the plant. The monitors provided hourly-averaged concentrations. The 

important terrain features rise approximately 250 m to 330 m above stack base at about 2 to 3 

km downwind from the stack. Meteorological data include winds, turbulence, and ΔT from a 

towe instrumented at 10 m, 50 m, and 100 m. National Weather Service surface data were 

available from a station 45 km away.
110

Figure B-3. Lovett study area.
111

B.2.4 Westvaco

The Westvaco Corporation’s pulp and paper mill in rural Luke, Maryland is located in a 

complex terrain setting in the Potomac River valley (Figure B-4). A single 183-m buoyant 

source was modeled for this evaluation. There were 11 SO2 monitors surrounding the facility, 

with eight monitors well above stack top on the high terrain east and south of the mill at a 

distance of 800 - 1500 m. Hourly meteorological data (wind, temperature, and turbulence) were 

collected between December 1980 and November 1981 at three instrumented towers: the 100-m 

Beryl tower in the river valley about 400 m southwest of the facility; the 30-m Luke Hill tower 

on a ridge 900 meters north-northwest of the facility; and the 100-m Met tower located 900 m 

eastsoutheast of the facility on a ridge across the river.

Figure B-4. Westvaco study area.
112

B.2.5 DAEC

The Duane Arnold Energy Center (DAEC) is located in rural Iowa, located about 16 km 

northwest of Cedar Rapids. It is located in a river valley with some bluffs on the east side. 

Terrain varies by about 30 m across the receptor network with the eastern half of the 

semicircular receptor arcs being flat and the western half elevated. The tracer study35 involved 

SF6 releases from two rooftops (46-m and 24-m levels) and the ground (1-m level). Building 

tiers for the rooftop releases were 43 and 24 m high, respectively. The 1-m and 24-m releases 

were non-buoyant, non-momentum, while the 46-m release was close to ambient, but had about 

a 10 m/s exit velocity. The number of tracer release hours was 12, 16 and 11 from the release 

heights of 46 m, 24 m, and 1 m, respectively. There were two arcs of monitors at downwind 

distances of 300 and 1000 m (see Figure B-5). Meteorological data consisted of winds at 10, 24, 

and 50 m. The meteorological conditions were mostly convective (30 out of 39 hours), with 

fairly light wind speeds. Only one hour had a wind speed above 4 m/s (4.6) , and almost half of 

the hours were less than 2 m/s.
113

Figure B-5. DAEC study area (SF6 releases).

B.2.6 EOCR

The EOCR study involved the simultaneous release of three tracer gases (SF6, F12, and 

Freon-12B2) at three levels around the Experimental Organically Cooled Reactor (EOCR) test 
114

reactor building at the Idaho National Engineering Laboratory in Southeast Idaho. The terrain 

was flat with low-lying shrubs. The main building was 25 m high with an effective width of 25 

m. The tracer releases typically occurred simultaneously, and were conducted during 22 

separate time periods. Tracer sampler coverage was provided at eight concentric rings at 

distances of about 50, 100, 200, 400, 800, 1200, and 1600 m from the release points (see Figure 

B-6). The stability classes ranged from stable to unstable. The 10-m wind speeds for the cases 

selected ranged from 3 to 8 m/s.

Figure B-6. Terrain map featuring the entire EOCR grid with the source at the 

grid center (SF6 releases). Arcs are at distances of about 40, 80, 200, 400, 800, 

1200, and 1600 m.

B.2.7 Alaska

The Alaska North Slope tracer study (see Figure B-7) involved 44 hours of buoyant SF6 

releases from a 39-m high turbine stack. Tracer sampler coverage ranged over seven arcs from 

50 to 3,000 m downwind. Meteorological data, including wind speed, wind direction,
115

temperature, sigma-theta, and sigma-w, were available from an on-site tower at the 33-m level. 

Atmospheric stability and wind speed profiles were influenced by the smooth snow-covered 

tundra surface with negligible levels of solar radiation in the autumn months. All experiments 

(44 usable hours) were conducted during the abbreviated day light hours (0900 – 1600). Wind 

speeds taken at the 33-m level during the tests were less than 6 m/s during one and part of 

another test, between 6 and 15 m/s during four tests, and in excess of 15 m/s during three tests. 

Stability conditions were generally neutral or slightly stable.

Figure B-7. Depiction of Alaska North Slope Oil Gathering Center turbine 

stack, meteorological tower (X), and camera locations used to visualize plume 

rise.
116

B.2.8 Indianapolis

The Indianapolis study consisted of an elevated, buoyant tracer (SF6) released in a flatterrain urban to suburban area from a single 84-m stack (Figure B-8). Data are available for 

approximately a four- to five-week period with 177 monitors providing 1-hour averaged 

samples along arcs from 250 m to 12 km downwind for a total of 1,297 arc-hours. 

Meteorological data included wind speed and direction, sigma-theta on a 94-meter tower; and 

wind speed, ΔT (2m - 10m) and other supporting surface data at three other 10-m towers (Figure 

B-9). Observed plume rise and estimates of plume sigma-y are also available from the database.
117

Figure B-8. Map showing the location of the Perry-K Station (A), the Hoosier 

Dome (B), and the central Indianapolis business district (C). The downtown 

surface meteorological site is located at (D) and the "bank tower" site was on top 

of the building at (E).
118

Figure B-9. Indianapolis meteorological sites and emissions site (Perry K 

Station).
119

B.2.9 Kincaid

The Kincaid SO2 study was conducted in a flat rural area of Illinois. It involved a 

buoyant, continuous release of SO2 from a 187-m stack in rural flat terrain. The study included 

about six months of data between April 1980 and June 1981 (a total of 4,614 hours of samples). 

There were 30 SO2 monitoring stations providing 1-hour averaged samples from about 2 km to 

20 km downwind of the stack. Meteorological data included wind speed, direction, and 

temperature from a tower instrumented at 2, 10, 50, and 100 m levels, and nearby National 

Weather Service (NWS) data.

Figure B-10. Kincaid study area.
120

B.2.10 AGA

The AGA experiments39 occurred during spring and summer 1980 at gas compressor 

stations in Texas and Kansas. At each test facility, one of the gas compressor stacks was 

retrofitted to accommodate SF6 tracer gas emissions. In addition, stack height extensions were 

provided for some of the experiments (with the normal stack height close to 10 m). The stack 

height to building height ratios for the tests ranged from 0.95 to 2.52. There were a total of 63 

tracer releases over the course of the tests, and the tracer samplers were located between 50 and 

200 m away from the release point (see Figure B-11). An instrumented 10-m tower was 

operated at both experimental sites. The tracer releases were generally restricted to daytime 

hours. Stability classes range from neutral to extremely unstable, except for three hours that 

were slightly stable. Wind speeds range from 2 to 11 m/s over the 63 hours.

Figure B-11. Plan view of the locations of tracer samplers at Site 1, AGA field 

study (SF6 releases)
121

B.2.11 Millstone

The Millstone nuclear power plant is located on the Connecticut coast, near Niantic. 

The model evaluation database features 36 hours of SF6 emissions from a 48-m reactor stack and 

26 hours of Freon emissions from a 29-m turbine stack. Exit temperatures were close to 

ambient (about 295K) with exit velocities of about 10 m/s for both the reactor stack (48.3 m) 

and the three turbine stacks (29.1 m). These stacks were associated with 45-m and 28-m 

building tiers, respectively. The monitoring data consisted of three arcs at 350, 800 and 1,500 m. 

Meteorological data were available from an on-site tower at the 10-m and 43-m levels. There 

was about an even split between stable and unstable hours, with mostly on-shore winds and 

fairly high wind speeds. There were only 3 stable hours with wind speed less than 4 m/s, and 

the majority was above about 7 m/s and several above 10 m/s. Figure B-12 shows the layout of 

the study area.

Figure B-12. Millstone study area (SF6 and freon releases)
122

B.2.12 Bowline

The Bowline Point site33, located in the Hudson River valley in New York State, is 

shown in Figure B-13 (topographic map). The electric utility site included two 600-MW units, 

each with an 86.9-m stack and a dominant roof tier with a height of 65.2 m high in a rural area. 

There were four monitoring sites as shown in Figure B-13 that ranged from about 250 to 850 m 

from the stacks. Hourly emissions data was determined from load data, coal analyses, and sitespecific relationships between loads and fuel consumption. Meteorological data was obtained 

from a 100-m tower at the site. This site was also used as an independent evaluation database 

with the entire year included.

Figure B-13. Bowline Point study area (SO2 releases)

B.2.13 Baldwin

The Baldwin Power Plant is located in a rural, flat terrain setting of southwestern Illinois 

and has three identical 184-m stacks aligned approximately north-south with a horizontal 

spacing of about 100 m. There were 10 SO2 monitors that surrounded the facility, ranging in 
123

distance from two to ten km. On-site meteorological data was available during the study period 

of 1 April 1982 through 31 March 1983 and consisted of hourly averaged wind speed, wind 

direction, and temperature measurements taken at 10 m and wind speed and wind direction at 

100 m.

Figure B-14. Baldwin study area.
124

B.2.14 Clifty Creek

The Clifty Creek Power Plant is located in rural southern Indiana along the Ohio River 

with emissions from three 208-m stacks during this study (Figure B-15). The area immediately 

north of the facility is characterized by cliffs rising about 115 m above the river and intersected 

by creek valleys. Six nearby SO2 monitors (out to 16 km from the stacks) provided hourly 

averaged concentration data. Meteorological data from a nearby 60-m tower covered the twoyear period from 1 January 1975 through 31 December 1976, although only the data from 1975 

were used in this evaluation. This database was also used in a major EPA-funded evaluation of 

rural air quality dispersion models in the early 1980s.
125

Figure B-15. Clifty Creek study area.
126

B.3 Evaluation methodology

B.3.1 AERMET/AERMOD comparisons

Two versions of AERMET/AERMOD will be compared using Robust highest 

concentrations and the EPA Protocol for determining best performing model. AERMET 

15181/AERMOD 15181 will be compared against AERMET 16216/AERMOD 16216 with 

various combinations of adjusted or non-adjusted surface friction velocity (u*) and 

inclusion/exclusion of turbulence parameters (sv and sw). The modeled scenarios are:

• 15181_no_u*_with_turb: AERMET/AERMOD 15181 with no u* adjustment and 

turbulence included in the meteorological data

• 15181_no_u*_no_turb: AERMET/AERMOD 15181 with no u* adjustment and no 

turbuluence included in the meteorological data 

• 16216_no_u*_witih turb: AERMET/AERMOD 16216 with no u* adjustment and 

turbulence included in the meteorological data

• 16216_with_u*_no_turb: AERMET/AERMOD 16216 with u* adjustment and no 

turbulence included in the meteorological data.

• 16216_no_u*_no_turb:AERMET/AERMOD 16216 with no u* adjustment and no 

turbulence included in the meteorological data.

B.3.2 Evaluation procedures

B.3.2.1 Robust highest concentrations

Robust highest concentrations (RHC) were calculated for each averaging period of each 

database. The RHC statistic is calculated as:

HH = Χ(NN) + [Χ� − Χ(NN)] × ln �

3NN − 1

2 � (119)

where Χ(N) is the Nth largest value, X is the average of N-1 values, and N is the number of 
127

values exceeding the threshold value, usually 26. 

For the 1-hour RHC, the RHC is calculated based on N=26 across all modeled and 

monitored values, i.e. not paired in time or space. For the 3-hour and 24-hour the RHC is 

calculated separately for each monitor within the network for observations and modeled values. 

The highest observed RHC is then compared to the highest modeled RHC.

B.3.2.2 EPA Protocol for determining best performing model

AERMOD output among the different meteorological datasets was evaluated using the 

EPA Protocol for determining the best performing model, or Cox-Tikvart method (U.S. EPA, 

1992; Cox and Tikvart, 1990). The protocol uses a two-step process for determining the better 

performing model when comparing models. The first step is a screening test that fails to perform 

at a minimal operational level. The second test applies to those models that pass the screening 

test that uses bootstrapping to generate a probability distribution of feasible outcomes (U.S. 

EPA, 1992). This section will discuss the methodology using the evaluation cases as examples.

The first step is to perform a screening test based on fractional bias:

FF = 2 �

OO − PP

OO + PP � (120)

where FB is the fractional bias, OB is the average of the highest 25 observed concentrations and 

PR is the average of the highest 25 predicted averages. The fractional bias is also calculated for 

the standard deviation where OB and PR refer to the standard deviation of the highest 25 

observed and predicted concentrations respectively. This is done across all monitors and 

modeled receptors, unpaired in time and space for the 3-hour and 24-hour averaging periods. The 

fractional bias of the means is plotted against the fractional bias of the standard deviation. Biases 

that exceed a factor-of-two under-prediction or over-prediction are considered grounds for 

excluding a model for further evaluation (U.S. EPA, 1992).

Models that pass the screening test are subjected to a more comprehensive statistical 

comparison that involves both an operational and scientific component using the RHC (Eq. 1). 
128

For the evaluations presented here, the screening step was skipped. The operational component 

is to measure the model’s ability to estimate concentration statistics most directly used for 

regulatory purposes and the scientific component evaluates the model’s ability to perform 

accurately throughout the range of meteorological conditions and the geographic area of 

concern (U.S. EPA, 1992).

The operational component of the evaluation compares performance in terms of the 

largest network-wide RHC test statistic. The RHC is calculated separately for each monitor 

within the network for observations and modeled values. The highest observed RHC is then 

compared to the highest modeled RHC using equation 2, where RHC now replaces the means of 

the top 25 values of observed or modeled concentrations. Absolute fractional bias (the absolute 

value of fractional bias), AFB is calculated for 3 and 24-hour averages.

The scientific component of the evaluation is also based on absolute fractional bias but 

the bias is calculated using the RHC for each meteorological condition and monitor. The 

meteorological conditions are a function of atmospheric stability and wind speed. For the 

purposes of these studies, six unique conditions were defined based on two wind speed 

categories (below and above 2.0 m/s) and three stability categories: unstable, neutral, and stable. 

1 In this evaluation, only 1-hour concentrations are used and the AFB is based on RHC values 

paired in space and stability/wind speed combination. 

A composite performance measure (CPM) is calculated from the 1-hour, 3-hour, and 24-

hour AFB’s:

CC CC = 1

3

× �AA ,jj� +

2

3

× �

AA 3 − AA 24

2 � (121)

where AFBi,j is the absolute fractional bias for monitor i and meteorological condition j, AA ,jj

is the average absolute fractional bias across all monitors and meteorological conditions, AFB3 is 

 

1 In U.S. EPA (1992), the three stability categories are related to the Pasquill-Gifford categories, unstable 

being A, B, and C, neutral being D, and stable being E and F. Since AERMOD does not use the stability categories, 

the stability class was determined using Monin-Obukhov length and surface roughness using methodology from 

AERMOD subroutine LTOPG.
129

the absolute fractional bias for the 3-hour average, and AFB24 is the absolute fractional bias for 

the 24-hour average. Once CPM values have been calculated for each model, a model 

comparison measure is calculated to compare the models:

AA, = CC CCAA − CC CC (122)

where CPMA is the CPM for model A and CPMB is the CPM for model B. When more than two 

models are being compared simultaneously, the number of MCM values is equal to the total of 

the number of unique combinations of two models. For Martins Creek, Lovett, Westvaco, and 

Kincaid, there are four scenarios each, so there were six MCM comparisons for each location. 

For Bowline, Baldwin, and Clifty Creek, there are three scenarios each, resulting in three MCM 

comparisons for each location. 

In order to determine if the difference between models was statistically significant, the 

standard error was calculated. A bootstrapping technique was used to create 1000 sample years 

based on methodology outlined in U.S. EPA (1992). The original data is divided into 3-day 

blocks. Within each season, the 3-day blocks are sampled with replacement until a total season 

is created. The process is repeated until a 1000 boot-strap years are created2

. The standard 

error is calculated as the standard deviation of the bootstrap generated outcomes for the MCM.

The magnitude and sign of the MCM are indicative of relative performance of each pair 

of models. The smaller the CPM the better the overall performance of the model. This means 

that for two models, A and B, a negative difference between the CPM for A and CPM for B 

implies that model A is performing better (Model A has a smaller CPM) while a positive 

difference indicates that Model B is performing better.

Since more than two scenarios are being evaluated in these studies, simultaneous 

confidence intervals of 90 and 95 percent were calculated. These were calculated by finding the 

 

2 The bootstrapping was completed using the SAS© SURVEYSELECT procedure with resampling for 1000 

replicates.
130

90th and 95th percentiles of the distribution across all MCM values from the bootstrapping 

procedure for all model comparisons. The confidence intervals were then found by:

CCCCXX,AA, = AA, ± ccXXssAA, (123)

where CIX,A,B is the confidence interval for X percent (90 or 95th) for models A and B, MCMA,B

is as defined in Equation 4, cX is the X percentile of the MCM values from the bootstrap results 

and sA,B is the standard deviation of the bootstrap MCM results for models A and B. Note that in 

Equation 5, MCMA,B is the MCM value from the original data, not the bootstrap results.

For each pair of model comparisons, the significance of the model comparison measure 

depended on whether the confidence interval overlapped zero. If the confidence interval 

overlapped zero, then the two models were not performing at a level which was considered 

statistically different. Otherwise, if they did not overlap zero, then there was a statistically 

significant difference between the two models.

B.4. Results

B.4.1 Turbulence cases

Table 2 lists the hourly observed and modeled RHC, as well as 3-hour and 24-hour RHC 

for applicable database, for the databases that initially included turbulence. Table 3 lists the 

RHC values for those databases initially without turbulence. The modeled scenario(s) closest to 

the observed RHC are highlighted in gray for each database. 

Results in Table 2 indicate that with the exception of DAEC (24 m release) and EOCR, 

the 16216 and 18081 modeled RHC’s are identical. For DAEC and EOCR, the differences 

between 16216 and 18081 are due the correction of the downwash bug that was in 16216 and 

corrected in 18081. Results in Table 2 also indicate that for the most part for the databases with 

turbulence data, the 16216 or 18081 cases without the u* adjustment and with turbulence data 

were the better performers against observations. For a few instances, depending on the 

averaging period, the cases with the u* adjustment and no turbulence, or the cases with no u* 

adjustment and no turbulence were the better performers.
131

Table 3 indicates that for the non-turbulence databases, the use of adjusted u* increased 

modeled performance in some cases depending on the averaging period or stack height. while 

decreasing or not changing model performance in other cases, depending on averaging period or 

stack height. For the databases that had multiple averaging periods (Martins Creek, Lovett, 

Westvaco, and Kincaid), there was not a consistent better performing model across the averaging 

periods. For example, for Martins Creek, 16216_with_u*_no_turb performed better for the 1 

and 24-hour averaging period, while 16216_no_u*_with_turb performed better for the 3-hour 

period. For DAEC, which had observed concentrations for emissions from different stack 

heights, the better performing modeling appeared to be dependent on stack height. Overall, it

appears that the use of adjusted u* did not increase model performance for most of the cases and 

that the inclusion of turbulence is more important to model performance than the u* adjustment.
132

Table 2. Hourly, 3-hour, and 24-hour RHC for turbulence cases. Best performing model compared to observed RHC are highlighted 

in gray.

Database

Avg. 

period 

(hr)

RHC

Observed

AERMOD version

16216 18081

No_u*_with 

turb

With_u*_no_

turb

No_u*_

no_turb

No_u*_

with 

turb

With_u*

_no_tur

b

No_u*_

no_turb

Martins Creek 1 1216 1133 1034 1427 1133 1034 1427

3 461 497 505 655 497 505 655

24 79 141 129 156 141 129 156

Tracy 1 15 13 18 25 13 18 25

Lovett 1 426 374 538 622 374 538 622

3 187 169 239 254 169 239 254

24 52 48 63 68 48 63 68

Westvaco 1 2757 2460 1252 2091 2460 1252 2091

3 1575 1731 783 1654 1731 783 1654

24 480 524 457 615 524 457 615

DAEC (h=1m) 1 346 240 188 222 240 188 222

DAEC (h=24m) 1 253 88 70 74 84 71 75

DAEC (h=46m) 1 140 91 59 99 91 59 99

EOCR 1 3763 5797 5712 8225 5820 5752 8248

Alaska 1 6 5 8 8 5 8 8

Indianapolis 1 6 4 4 5 4 4 5

Kincaid 1 1611 1312 717 717 1312 717 717

3 618 635 470 470 635 470 470

24 113 103 167 167 103 167 167
133

B.4.2 Non-turbulence cases

Table 3 lists the RHC values for the non-turbulence databases for 16216 and 18081. In 

these databases, because of the lack of turbulence in the meteorological data, the effect of the u* 

adjustment has more impact in improving model performance. Also, the results indicate the 

changes made to AERMOD between 16216 and 18081 did not impact these findings.

Table 3. Hourly, 3-hour, and 24-hour RHC for non-turbulence cases. Best performing model 

compared to observed RHC are highlighted in gray.

Database

Avg. 

period 

(hr)

RHC

Observed

AERMOD version

16216 18081

With_u*_

no_turb

No_u*_

no_turb

With_u*

_no_turb

No_u*_

no_turb

AGA 1 296 262 281 262 281

Millston 

(Freon)

1 76 96 101 96 101

Millston 

(SF6)

1 79 33 35 33 35

Bowline 1 763 552 547 552 547

3 469 514 522 514 523

24 204 307 290 307 290

Baldwin 1 2348 3531 3531 3531 3531

3 920 1183 1184 1183 1184

24 209 230 230 230 230

Clifty Creek 1 1451 1360 1360 1360 1360

3 796 871 870 871 870

24 243 170 165 170 165

B.4.3 Statistical evaluations

While the review of RHC can indicate general model performance, the use of the EPA 

Protocol for Determining Best Model provides a statistical basis of determing the best 

performing model. Tables 4 and 5 show the composite performance measure (CPM) for the 

turbulence databases and non-turbulence databases respectively. For the databases with 
134

turbulence (Table 4), the best performing models for Martins Creek were the cases with adjusted 

u* and not turbulence but for the remaining areas, the better performing models were the 

adjusted u* and no turbulence scenarios. This means the use of adjusted u* did not increase 

model performance and the use of turbulence was important to model performance. For the nonturbulence databases (Table 5), the use of adjusted u* increased model performance for Baldwin 

and Clifty Creek, while for Bowline, the use of adjusted u* slightly decreased model 

performance. For all cases, the CPM values were identical for the 16216 and 18081 model 

versions, suggesting the changes between 16216 and 18081 had minimal to no impact on model 

performance, which was expected based on the changes made to AERMET and AERMOD and 

no changes to the adjusted u* equations.

Table 4. Composite Performance Measure (CPM) for turbulence cases. Scenarios with lowest 

CPM’s for each study location are highlighted in gray.

Scenario

Database

Martins Creek Lovett Westvaco Kincaid

16216_no_u*_with_turb 0.35 0.40 0.41 0.37

16216_with_u*_no_turb 0.31 0.52 0.60 0.56

16216_no_u*_no_turb 0.49 0.58 0.44 0.56

18081_no_u*_with_turb 0.35 0.40 0.41 0.37

18081_with_u*_no_turb 0.31 0.52 0.60 0.56

18081_no_u*_no_turb 0.49 0.58 0.44 0.56

Table 5. Composite Performance Measure (CPM) for non-turbulence databases. Scenarios 

with lowest CPM’s for each study location are highlighted in gray.

Scenario

Database

Bowline Baldwin Clifty Creek

16216_no_u*_no_turb 0.47 0.46 0.51

16216_with_u*_no_turb 0.50 0.45 0.49

18081_no_u*_no_turb 0.47 0.46 0.51

18081_with_u*_no_turb 0.50 0.45 0.49

Tables 6 through 9 show the model comparison measure (MCM) for the turbulence 

databases while Tables 10 through 12 show the MCM for the non-turbulence databases. Also 
135

shown are the 90 and 95% confidence intervals of the MCM based on the bootstrapping results. 

Confidence intervals highlighted in gray indicated statistical significance in the specific MCM 

cases. The original pairings of 16216 scenarios to other 16216 scenarios are shown for 

comparison to the analogous 18081 pairings. MCM pairings for the same u*/turbulence pairings 

between 18081 and 16216 are also shown to show if model changes made differences to results. 

For most cases, such comparisons are zero or very small. Non-zero differences in these cases are 

a result of the downwash bug fix in 18081. 

Martins Creek (Table 6): The better performing models were 16216 and 18081 with u* 

and no turbulence. Also, the MCM results indicate that the use of adjusted u* with no turbulence 

is not statistically significant when compared to no adjusted u* with turbulence for both 16216 

and 18081. There were three statistically significant MCM pairings that were statistically 

significant at the 90% confidence interval and these were the difference between no u* 

adjustment and no turbulence and the cases (no adjusted u* with turbulence or adjusted u* with 

no turbulence) for both 16216 and 18081, indicating that not using adjusted u* and not using 

turbulence noticeably decreases model performance. At the 95% confidence interval, the only 

statistically significant difference was between no adjusted u*/ no turbulence and adjusted u*/no 

turbulence for both 16216 and 18081. 

Lovett (Table 7): All cases of AERMET/AERMOD 18081 are statistically insignificant 

when compared AERMET/AERMOD 16216 at both the 90% and 95% CI, all18081 cases are 

statistically insignificant compared to each other, similar to the 16216 cases. 

Westvaco (Table 8): The use of adjusted u* decreases model performance significantly 

at both the 90% and 95% CI for both 16216 and 18081. The use of no adjusted u* and no 

turbulence also decreases model performance at a statistically significant level for both 16216 

and 18081.

Kincaid (Table 9): None of the MCM differences were statistically significant at 90% or 

95% CI. The better performers were 16216 or 18081 with no u* adjustment and inclusion of 

turbulence, but as previously stated, were not statistically different from the adjusted u* case or 

the case with no adjusted u* and no turbulence. 
136

For the non-turbulence databases (Tables 10-12), the use of adjusted u* was statistically 

insignificant compared to not using adjusted u* and the only differences between 16216 and 

18081 were due to the correction of the downwash bug in 16216.
137

Table 6. Martins Creek Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are 

significant at that percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower 

bound

Upper 

bound

Lower 

bound

Upper 

bound

16216_with_u*_no_turb - 16216_no_u*_with_turb -0.04 -0.14 0.07 -0.18 0.10

16216_no_u*_no_turb - 16216_no_u*_with_turb 0.14 0.02 0.27 -0.01 0.30

16216_no_u*_no_turb - 16216_with_u*_no_turb 0.18 0.05 0.32 0.01 0.35

18081_no_u*_no_turb - 16216_no_u*_no_turb 0 -0.12 0.12 -0.15 0.15

18081_no_u*_with_turb - 16216_no_u*_with_turb 0 -0.10 0.10 -0.13 0.13

18081_with_u*_no_turb – 16216_with_u*_no_turb 0 -0.11 0.11 -0.14 0.14

18081_with_u*_no_turb - 18081_no_u*_with_turb -0.04 -0.15 0.07 -0.18 0.10

18081_no_u*_no_turb - 18081_no_u*_with_turb 0.14 0.02 0.27 -0.01 0.3

18081_no_u*_no_turb - 18081_with_u*_no_turb 0.18 0.05 0.32 0.01 0.35
138

Table 7. Lovett Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are significant 

at that percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower 

bound Upper bound

Upper 

bound

Lower 

bound

16216_with_u*_no_turb - 16216_no_u*_with_turb 0.13 -0.06 0.31 -0.11 0.36

16216_no_u*_no_turb - 16216_no_u*_with_turb 0.18 -0.02 0.38 -0.07 0.43

16216_no_u*_no_turb - 16216_with_u*_no_turb 0.05 -0.10 0.20 -0.13 0.24

18081_no_u*_no_turb - 16216_no_u*_no_turb 0 -0.14 0.14 -0.18 0.18

18081_no_u*_with_turb - 16216_no_u*_with_turb 0 -0.11 0.11 -0.15 0.15

18081_with_u*_no_turb – 16216_with_u*_no_turb 0 -0.12 0.17 -0.15 0.15

18081_with_u*_no_turb - 18081_no_u*_with_turb 0.13 -0.4 0.29 -0.08 0.34

18081_no_u*_no_turb - 18081_no_u*_with_turb 0.18 -0.001 0.36 -0.05 0.41

18081_no_u*_no_turb - 18081_with_u*_no_turb 0.05 -0.08 0.19 -0.12 0.22
139

Table 8. Westvaco Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are significant at that 

percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower bound

Upper 

bound

Lower 

bound Upper bound

16216_with_u*_no_turb - 16216_no_u*_with_turb 0.19 0.09 0.29 0.05 0.32

16216_no_u*_no_turb - 16216_no_u*_with_turb 0.03 -0.03 0.08 -0.05 0.10

16216_no_u*_no_turb - 16216_with_u*_no_turb -0.16 -0.26 -0.06 -0.29 -0.03

18081_no_u*_no_turb - 16216_no_u*_no_turb 0 -0.05 0.05 -0.07 0.07

18081_no_u*_with_turb - 16216_no_u*_with_turb 0 -0.06 0.06 -0.08 0.08

18081_with_u*_no_turb – 16216_with_u*_no_turb 0 -0.05 0.05 -0.07 0.07

18081_with_u*_no_turb - 18081_no_u*_with_turb 0.19 0.08 0.30 0.04 0.33

18081_no_u*_no_turb - 18081_no_u*_with_turb 0.03 -0.03 0.08 -0.05 0.10

18081_no_u*_no_turb - 18081_with_u*_no_turb -0.l6 -0.26 -0.05 -0.30 -0.02
140

Table 9. Kincaid Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are significant at that 

percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower 

bound

Upper 

bound

Lower

bound Upper bound

16216_with_u*_no_turb - 16216_no_u*_with_turb 0.19 -0.14 0.52 -0.20 0.57

16216_no_u*_no_turb - 16216_no_u*_with_turb 0.19 -0.15 0.52 -0.20 0.58

16216_no_u*_no_turb - 16216_with_u*_no_turb -5.1x10-4 -0.15 0.15 -0.18 0.18

18081_no_u*_no_turb - 16216_no_u*_no_turb 0 -0.13 0.13 -0.16 0.16

18081_no_u*_with_turb - 16216_no_u*_with_turb 0 -0.34 0.34 -0.41 0.41

18081_with_u*_no_turb – 16216_with_u*_no_turb 0 -0.13 0.13 -0.16 0.16

18081_with_u*_no_turb - 18081_no_u*_with_turb 0.19 -0.09 0.47 -0.15 0.53

18081_no_u*_no_turb - 18081_no_u*_with_turb 0.19 -0.09 0.47 -0.15 0.53

18081_no_u*_no_turb - 18081_with_u*_no_turb -5.1x10-4 -0.13 0.12 -0.15 0.15
141

Table 10. Bowline Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are significant at that 

percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower 

bound

Upper 

bound Lower bound

Upper 

bound

16216_no_u*_no_turb -16216_with_u*_no_turb -0.03 -0.11 0.05 -0.13 0.07

18081_no_u*_no_turb - 16216_no_u*_no_turb 2.7x10-5 -0.08 0.07 -0.10 0.10

18081_with_u*_no_turb - 16216_with_u*_no_turb 7.4x10-5 -0.07 0.07 -0.09 0.09

18081_no_u*_no_turb - 18081_with_u*_no_turb -0.03 -0.10 0.04 -0.12 0.06
142

Table 11. Baldwin Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are significant at that 

percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower 

bound

Upper 

bound

Lower 

bound Upper bound

16216_no_u*_no_turb - 16216_with_u*_no_turb 0.002 -0.14 0.14 -0. 19 0.20

18081_no_u*_no_turb - 16216_no_u*_no_turb 0 -0.13 0.13 -0.18 0.18

18081_with_u*_no_turb - 16216_with_u*_no_turb 0 -0.12 0.12 -0.17 0.17

18081_no_u*_no_turb - 18081_with_u*_no_turb 0.002 -0.13 0.13 -0.18 0.18
143

Table 12. Clifty Creek Model Comparison Measure (MCM) results. Confidence intervals highlighted in gray are significant at that 

percent.

MCM Comparison MCM

Confidence Intervals

90% 95%

Lower 

bound

Upper 

bound

Lower 

bound Upper bound

16216_no_u*_no_turb - 16216_with_u*_no_turb 0.02 -0.06 0.09 -0.08 0.11

18081_no_u*_no_turb - 16216_no_u*_no_turb -4.2x10-8 -0.07 0.07 -0.09 0.09

18081_with_u*_no_turb - 16216_with_u*_no_turb 0 -0.07 0.07 -0.09 0.09

18081_no_u*_no_turb - 18081_with_u*_no_turb 0.02 -0.06 0.09 -0.08 0.11
144

B.5 Summary/Conclusions

Based on the results the RHC comparisons and the EPA protocol for determining best 

performing model, in situations involving turbulence, the use of turbulence without adjusting u* 

usually led to better performance than using adjusted u* without turbulence, especially in areas 

of complex terrain. In some instances, the differences between the adjusted u* cases were 

statistically worse than non-adjusted u* cases. For situations where turbulence is not in the 

meteorological data, the use of adjusted u* often resulted in little change or some increase in 

model performance. However, the databases without turbulence were in flat terrain and had talk

stacks, so model performance for non-turbulence cases with complex terrain cannot be 

determined from these results. The results of the RHC and EPA protocol also indicate that 

changes made to AERMOD 18081 had no unexpected changes from AERMOD 16216, which 

was expected. The only differences were due to the correction of the downwash bug in 16216 

and those differences were not unexpected. 
145

Appendix C: Evaluation of modified urban option

The urban option in AERMOD was modified, beginning with version 11059, to address 

potential issues associated with the transition from the nighttime urban boundary layer to the 

daytime convective boundary layer. Prior to version 11059, the enhanced dispersion associated 

with the urban nighttime heat island effect was ignored once the boundary layer turned 

convective. This could result in an unrealistic drop in the mixed layer height during early 

morning hours for urban applications, which could contribute to unrealistically high 

concentrations for low-level plumes. This effect was observed in the application of AERMOD 

for the Risk and Exposure Assessment (REA) in support of the NO2 NAAQS review in 

association with mobile source emissions (EPA, 2008). Beginning with version 11059 the 

urban option in AERMOD continues application of the urban boundary layer approach for urban 

sources until the daytime convective mixing height exceeds the urban nighttime mixing height, 

based on the user-specified population (EPA, 2018a)). This revision to AERMOD will 

generally reduce concentrations during the early morning transition to convective conditions for 

low-level urban sources, but may increase daytime concentrations for elevated urban sources. 

The modified implementation of the urban option was evaluated using data from the 

1985 Indianapolis SF6 urban field study (Perry, et al, 2005), and model-to-monitor comparisons 

at four ambient monitors for 2002 from the Atlanta NO2 REA (EPA, 2008). The Indianapolis 

study involved an elevated buoyant release and the Atlanta REA study involved mostly lowlevel mobile source emissions. Results from the Indianapolis study are presented in the form of 

Q-Q plots of ranked 1-hour modeled vs. observed concentrations, unpaired in time and space. 

Figure E-1 shows model performance for all stabilities and Figure E-2 shows model 

performance for convective conditions only. The revised urban option does not affect results 

under stable conditions. Results from the Atlanta NO2 REA are also presented in the form of QQ plots of 1-hour ranked modeled vs. observed concentrations, unpaired in time, for each of the 

four ambient NO2 monitors, shown in Figures E-3 through E-6. Both of these evaluations show 

improved model performance with the modified urban option in AERMOD.
146

Figure E-1. 1-hour Q-Q Plot for Indianapolis SF6 Study for All Stabilities

Figure E-2. 1-hour Q-Q Plot for Indianapolis SF6 Study for Convective Conditions

0.01

0.1

1

10

0.01 0.1 1 10

PREDICTED

OBSERVED

INDIANAPOLIS SF6 1-HR Q-Q PLOT (CONC) - All Stabilities

AERMOD v09292

AERMOD v11059
147

Figure E-3. 1-hour Q-Q Plot for Atlanta NO2 Study for Monitor 0002

Figure E-4. 1-hour Q-Q Plot for Atlanta NO2 Study for Monitor 3001

Atlanta NO2 Study - 1-hr QQ Plot for Monitor 0002 - Urban Transition Adjustment

0

50

100

150

200

0 50 100 150 200

Observed Conc (ppb)

Predicted Conc (ppb)

0002 - New

0002 - Old

Atlanta NO2 Study - 1-hr QQ Plot for Monitor 3001 - Urban Transition Adjustment

0

50

100

150

200

0 50 100 150 200

Observed Conc (ppb)

Predicted Conc (ppb)

3001 - New

3001 - Old
148

Figure E-5. 1-hour Q-Q Plot for Atlanta NO2 Study for Monitor 0048

Figure E-6. 1-hour Q-Q Plot for Atlanta NO2 Study for Monitor JST

Atlanta NO2 Study - 1-hr QQ Plot for Monitor 0048 - Urban Transition Adjustment

0

50

100

150

200

0 50 100 150 200

Observed Conc (ppb)

Predicted Conc (ppb)

0048 - New 

0048 - Old

Atlanta NO2 Study - 1-hr QQ Plot for Monitor JST - Urban Transition Adjustment

0

50

100

150

200

0 50 100 150 200

Observed Conc (ppb)

Predicted Conc (ppb)

JST - New

JST - Old
1
2
3

References

AECOM, 2010: AERMOD Low Wind Speed Evaluation Study Results. Prepared for the 

American Petroleum Institute and Utility Air Regulatory Group. Prepared by AECOM, 

Westford, MA. March 22, 2010.

Andre, J. C. and L. Mahrt, 1982: The nocturnal surface inversion and influence of clean-air 

radiative cooling. J.Atmos.Sci., 39, 864-878.

Baerentsen, J. H. and R. Berkowicz, 1984: Monte Carlo simulation of plume dispersion in the 

convective boundary layer. Atmos.Environ., 18, 701-712.

Bange, P., L. Jannsen, F. Nieuwstadt, H. Visser, and J. Erbrink, 1991. “Improvement of the 

modeling of daytime nitrogen oxidation in plumes by using instantaneous plume 

dispersion parameters,” Atmos. Environ., 25A (10), 2321-2328.

Barad, M. L., 1958: Project Prairie Grass, A Field Program in Diffusion. Geophysical Research 

Papers, No. 59, Vols. I and II, AFCRC-TR-58-235, Air Force Cambridge Research 

Center, 439pp.

Berkowicz, R., Olesen, J. R., and Torp, U., 1986: The Danish Gaussian air pollution model 

(OLM): Description, test and sensitivity analysis, in view of regulatory applications. Air 

Pollution Modeling and Its Application. De Wispelaire, V. C., Schiermeier, F. A., and 

Grillani, N. V, Plemum, 453-481pp.

Bornstein, R. D., 1968: Observations of urban heat island effects in New York City. 

J.Appl.Meteor., 7, 575-582.

Bowne, N. E., R. J. Londergan, D. R. Murray, and H. S. Borenstein, 1983: Overview, Results, 

and Conclusions for the EPRI Plume Model Validation and Development Project: Plains 

Site. EPRI Report EA-3074, Project 1616-1, Electric Power Research Institute, Palo Alto, 

CA, 234 pp. 1983.
4

Brett, A. C. and S. E. Tuller, 1991: Autocorrelation of hourly wind speed observations. 

J.Appl.Meteor., 30, 823-833.

Briggs, G. A., 1969: Plume rise. USAEC Critical Review Series, TID-25075, NTIS, 81pp.

Briggs, G. A., 1971: Some recent analyses of plume rise observations. Proceedings of the Second 

International Clean Air Congress. Englund, H. M. and Berry, W. T., Academic Press, 

1029-1032.

Briggs, G. A., 1973: Diffusion Estimation for Small Emissions. 1973 Annual Report, ATDL106, Air Resources Atmospheric Turbulence and Diffusion Laboratory, Environmental 

Res. Lab., NOAA, Oak Ridge, TN. 

Briggs, G. A., 1975: Plume rise predictions. Lectures on Air Pollution and Environmental Impact 

Analysis. Haugen, D. A., American Meteorological Society, 59-111pp.

Briggs, G. A., 1984: Plume rise and buoyancy effects. Atmospheric Science and Power 

Production. Randerson, D., U.S. Dept. of Energy, 327-366pp.

Briggs, G. A., 1988: Analysis of diffusion field experiments. Lectures on Air Pollution 

Modeling. Venkatram, A. and Wyngaard, J. C., American Meteorological Society, 

63-117.

Briggs, G. A., 1993: Plume dispersion in the convective boundary layer. Part II: Analysis of 

CONDORS field experiment data. J.Appl.Meteor., 32, 1388-1425.

Brode, R. W., 2002: Evaluation of the AERMOD Dispersion Model. 

Brost, R. A., J. C. Wyngaard, and D. H. Lenschow, 1982: Marine stratocumulus layers: Part II: 

Turbulence budgets. J.Atmos.Sci., 39, 818-836.

Businger, J. A., 1973: Turbulent transfer in the atmospheric surface layer. Workshop on 

Micrometeorology. Haugen, D. A., American Meteorological Society. 
5

Businger, J. A., J. C. Wyngaard, Y. Izumi, and E. F. Bradley, 1971: Flux-profile relationships in 

the atmospheric surface layer. J.Atmos.Sci., 28, 181-189.

Carruthers, D. J., and Coauthors, 1992: UK atmospheric dispersion modelling system. Air 

Pollution Modeling and Its Application. Plenum Press, New York.

Carson, D. J., 1973: The development of a dry inversion-capped convectively unstable boundary 

layer. Quart.J.Roy.Meteor.Soc., 99, 450-467.

Caughey, S. J. and S. G. Palmer, 1979: Some aspects of turbulence structure through the depth of 

the convective boundary layer. Quart.J.Roy.Meteor.Soc., 105, 811-827.

Cimorelli, A. J., S. G. Perry, R. F. Lee, R. J. Paine, A. Venkatram, J. C. Weil, and R. B. Wilson, 

1996: Current Progress in the AERMIC Model Development Program. Preprints, 89th

Annual Meeting Air and Waste Management Association, Air and Waste Management 

Association, Pittsburgh, PA.

Cimorelli, A. J., S. G. Perry, A. Venkatram, J. C. Weil, R. J. Paine, R. B. Wilson, R. F. Lee, 

W. D. Peters, R. W. Brode, and J. O. Paumier, 2002: AERMOD: Description of Model 

Formulation (Version 02222). EPA 454/R-02-002d. U. S. Environmental Protection 

Agency, Research Triangle Park, NC.

Cimorelli, A. J., S. G. Perry, A. Venkatram, J. C. Weil, R. J. Paine, R. B. Wilson, R. F. Lee, W. 

D. Peters, and R. W. Brode, 2005: AERMOD: A dispersion model for industrial source 

applications Part I: General model formulation and boundary layer characterization. 

J.Appl.Meteor. 44, 682-693

Clarke, R. H., A. J. Dyer, R. R. Brook, D. G. Reid, and A. J. Troop, 1971: The Wangara 

experiment: boundary layer data. Technical Report No. 19, Division of Meteorological 

Physics CSIRO, Australia. 

Cole, H.S. and J.E. Summerhays, 1979. “A review of techniques available for estimating shortterm NO2 concentrations,” J. Air Pollut. Control Assoc., 29(8), 812-817.

Collier, L. R. and J. G. Lockwood, 1975: Reply to comment. Quart.J.Roy.Meteor.Soc., 101, 
6

390-392.

Cox, W.M. and J.A. Tikvart, 1990: A statistical Procedure for Determining the Best Performing 

Air Quality Simulation Model. Atmos. Environ., 24A (9): 2387-2395.

Deardorff, J. W., 1970: Convective velocity and temperature scales for the unstable boundary 

layer for Rayleigh convection. J.Atmos.Sci., 27, 1211-1213.

Deardorff, J. W., 1972: Numerical investigation of neutral and unstable planetary boundary 

layers. J.Atmos.Sci., 29, 91-115.

Deardorff, J. W., 1979: Prediction of convective mixed-layer entrainment for realistic capping 

inversion structure. J.Atmos.Sci., 36, 424-436.

Deardorff, J. W., 1980: Progress in Understanding Entrainment at the Top of a Mixed Layer. 

Preprints, Workshop on the Planetary Boundary Layer, American Meteorological 

Society, Boston, MA. 

DiCristofaro, D. C. et al., 1985: EPA Complex Terrain Model Development: Fifth Milestone 

Report - 1985. EPA-600/3-85-069, U. S. Environmental Protection Agency, Research 

Triangle Park, North Carolina.

Dyer, A. J., 1974: A review of flux-profile relationships. Bound.Layer Meteor., 7, 363-372.

Garratt, F. R., 1992: The Atmospheric Boundary Layer. Cambridge University Press, New York, 

New York, 334pp.

Gifford, F. A., 1961: Uses of routine meteorological observations for estimating atmospheric 

dispersion. Nuclear Safety, 2, 47-51.

Hanna, S. R., 1983: Lateral turbulence intensity and plume meandering during stable conditions. 

J.Appl.Meteor., 22, 1424-1430.
7

Hanna, S. R. and J. S. Chang, 1991: Modification of the Hybrid Plume Dispersion Model 

(HPDM) for urban conditions and its evaluation using the Indianapolis data set, Volume 

III: Analysis of urban boundary layer data. EPRI Project No. RP-02736-1, Electric Power 

Research Institute, Palo Alto, CA.

Hanna, S. R. and J. S. Chang, 1993: Hybrid Plume Dispersion Model (HPDM), improvements 

and testing at three field sites. Atmos.Environ., 27A, 1491-1508.

Hanna, S. R. and R. J. Paine, 1989: Hybrid Plume Dispersion Model (HPDM) development and 

evaluation. J.Appl.Meteor., 28, 206-224.

Hanna, S. R., J. C. Weil, and R. J. Paine, 1986: Plume Model Development and Evaluation -

Hybrid Approach. EPRI Contract No. RP-1616-27, Electric Power Research Institute, 

Palo Alto, CA.

Hanrahan, P.L., 1999a. “The plume volume molar ratio method for determining NO2/NOx ratios 

in modeling. Part I: Methodology,” J. Air & Waste Manage. Assoc., 49, 1324-1331.

Haugen, D. A. (Editor), 1959: Project Prairie Grass, A field program in diffusion. Geophysical 

Research Paper, No. 59, Vol. III. Report AFCRC-TR-58-235, Air Force Cambridge 

Research Center, 439 pp.

Hayes, S. R. and G. E. Moore, 1986: Air quality model performance: a comparative analysis of 

15 model evaluation studies. Atmos.Environ., 20, 1897-1911.

Hicks, B. B., 1985: Behavior of turbulent statistics in the convective boundary layer. 

J.Appl.Meteor., 24, 607-614.

Holtslag, A. A. M., 1984: Estimates of diabatic wind speed profiles from near-surface weather 

observations. Bound.Layer Meteor., 29, 225-250.

Holtslag, A. A. M. and A. P. van Ulden, 1983: A simple scheme for daytime estimates for the 

surface fluxes from routine weather data. J.Climate Appl.Meteor., 22, 517-529.
8

Irwin, J. S., J. O. Paumier, and R. W. Brode, 1988: Meteorological Processor for Regulatory 

Models (MPRM) User's Guide. EPA-600/3-88-043, U.S. Environmental Protection 

Agency, RTP, NC. 

Izumi, Y., 1971: Kansas 1968 Field Program Data Report. No. 379, AFCRL-72-0041, Air Force 

Cambridge Research Laboratory, Bedford, MA, 79pp.

Kaimal, J. C., J. C. Wyngaard, D. A. Haugen, O. R. Cote', Y. Izumi, S. J. Caughey, and C. J. 

Readings, 1976: Turbulence structure in the convective boundary layer. J.Atmos.Sci., 33, 

2152-2169.

Kasten, F. and G. Czeplak, 1980: Solar and terrestial radiation dependent on the amount and type 

of cloud. Solar Energy, 24, 177-189.

Lamb, R. G., 1982: Diffusion in the convective boundary layer. Atmospheric Turbulence and Air 

Pollution Modelling. Nieuwstadt, F. T. M. and van Dop, H., Reidel, 159-229pp.

Lee, R. F., R. J. Paine, S. G. Perry, A. J. Cimorelli, J. C. Weil, A. Venkatram, and R. B. Wilson, 

1998: Developmental Evaluation of the AERMOD Dispersion Model. Preprints, 10th 

Joint AMS/AWMA Conference on Application of Air Pollution Meteorology, American 

Meteorological Society, Boston.

Lee, R. F., S. G. Perry, A. J. Cimorelli, R. J. Paine, A. Venkatram, J. C. Weil, and R. B. Wilson, 

1995: AERMOD - the developmental evaluation. Preprints, Tewnty-First NATO/CCMS 

International Technical Meeting on Air Pollution Modeling and Its Application, 

Baltimore, MD, U.S.A.

Liu, M. K., and G. E. Moore: 1984: Diagnostic validation of plume models at a plains site. EPRI 

Report No. EA-3077, Research Project 1616-9, Electric Power Research Institute, Palo 

Alto, CA. (1984)

Luhar, A.K., and K. N. Rayner, 2009: “Methods to Estimate Surface Fluxes of Momentum and 

Heat from Routine Weather Observations for Dispersion Applications under Stable 

Stratification”, Boundary-Layer Meteorology, 132, 437–454.
9

Misra, P. K., 1982: Dispersion of non-buoyant particles inside a convective boundary layer. 

Atmos.Environ., 16, 239-243.

Morton, B. R., G. I. Taylor, and J. S. Turner, 1956: Turbulent gravitational convection from 

maintained and instantaneous sources. Proc.Roy.Soc.London, A234, 1-23.

Nieuwstadt, F. T. M. and H. van Dop, 1982: Atmospheric Turbulence and Air Pollution 

Modelling. Reidel, 358pp.

NOAA, 1974: Technical Memorandum ERL ARL-52, 1974. “Diffusion under Low Wind Speed, 

Inversion Conditions.” Sagendorf, J. F., C. Dickson. Air Resources Laboratory, Idaho 

Falls, Idaho. 

NOAA, 1976: Technical Memorandum ERL ARL-61, 1976. “Diffusion under Low Wind Speed 

Conditions near Oak Ridge, Tennessee.” Wilson, R. B., G. Start, C. Dickson, N. Ricks. 

Air Resources Laboratory, Idaho Falls, Idaho.

Oke, T. R., 1973: City size and the urban heat island. Atmos.Environ., 7, 769-779.

Oke, T. R., 1978: Boundary Layer Climates. John Wiley and Sons, New York, New York, 

372pp.

Oke, T. R., 1982: The energetic basis of the urban heat island. Quart.J.Roy.Meteor.Soc., 108, 

1-24.

Oke, T. R., 1998: An algorithmic scheme to estimate hourly heat island magnitude. Preprints,

2nd Urban Environment Symposium, American Meteorological Society, Boston, MA,

80-83.

Paine, R. J. and B. A. Egan, 1987: User's guide to the Rough Terrain Diffusion Model (RTDM) -

Rev. 3.20. ERT Document PD-535-585, ENSR, Acton, MA, 260pp.
10

Paine, R. J. and S. B. Kendall, 1993: Comparison of observed profiles of winds, temperature, and 

turbulence with theoretical results. Preprints, Joint conference of the American 

Meteorological Society and Air & Waste Management Association Specialty Conference: 

The Role of Meteorology in Managing the Environment in the 90s, Scottsdale, AZ. 

Publication VIP-29, Air & Waste Management Association, Pittsburgh, PA.

Panofsky, H. A. and J. A. Dutton, 1984: Atmospheric Turbulence: Models and Methods for 

Engineering Applications. John Wiley and Sons, New York, 417pp.

Panofsky, H. A., H. Tennekes, D. H. Lenschow, and J. C. Wyngaard, 1977: The characteristics 

of turbulent velocity components in the surface layer under convective conditions. 

Bound.Layer Meteor., 11, 355-361.

Pasquill, F., 1961: The estimation of the dispersion of windborne material. Meteorol.Mag., 90, 

33-49.

Pasquill, F., 1976: Atmospheric dispersion parameters in Gaussian plume modeling - Part III: 

possible requirements for change in the Turner's Workbook values. EPA-600/4-76-030B, 

U. S. Environmental Protection Agency, Research Triangle Park, NC. 

Pasquill, F. and F. R. Smith, 1983: Atmospheric Diffusion. John Wiley and Sons Inc., New York, 

440pp.

Paumier, J. O., S. G. Perry, and D. J. Burns, 1992: CTDMPLUS: A dispersion model for sources 

near complex topography. Part II: Performance characteristics. J. Appl. Meteor., 31, 646–

660.

Perry, S. G., 1992: CTDMPLUS: A dispersion model for sources in complex topography. Part I: 

Technical formulations. J.Appl.Meteor., 31, 633-645.

Perry, S. G., D. J. Burns, R. J. Adams, R. J. Paine, M. G. Dennis, M. T. Mills, D. G. Strimaitis, 

R. J. Yamartino, and E. M. Insley, 1989: User's Guide to the Complex Terrain Dispersion 

Model Plus Algorithms for Unstable Situations (CTDMPLUS) Volume 1: Model 

Description and User Instructions. EPA/600/8-89/041, U.S. Environmental Protection 

Agency, RTP, NC, 196pp.
11

Perry, S. G., A. J. Cimorelli, R. F. Lee, R. J. Paine, A. Venkatram, J. C. Weil, and R. B. Wilson, 

1994: AERMOD: A dispersion model for industrial source applications. Preprints, 87th

Annual Meeting Air and Waste Management Association, Air and Waste Management 

Association, Pittsburgh, PA.

Perry, S. G., A. J. Cimorelli, R. J. Paine, R. W. Brode, J. C. Weil, A. Venkatram, R. B. Wilson, 

R. F. Lee, and W. D. Peters, 2005: AERMOD: A dispersion model for industrial source 

applications Part II: Model performance against seventeen field-study databases. 

J.Appl.Meteor. 44, 694-708.

Qian, W., and A. Venkatram, 2011: "Performance of Steady-State Dispersion Models UnderLow

Wind-Speed Conditions", Boundary Layer Meteorology, 138, 475-491.

Readings, C. J., D. A. Haugen, and J. C. Kaimal, 1974: The 1973 Minnesota atmospheric 

boundary layer experiment. Weather, 29, 309-312.

Schulman, L. L., D. G. Strimaitis, and J. S. Scire, 2000: Development and evaluation of the 

PRIME plume rise and building downwash model. Journal of the Air & Waste 

Management Association, 50, 378-390.

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

Smith, M. E., 1984: Review of the attributes and performance of 10 rural diffusion models.

Bull.Amer.Meteor.Soc., 65, 554-558.
12

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, 1985: The structure of the strongly stratified flow over hills: 

Dividing streamline concept. J.Fluid Mech., 152, 249-288.

Stull, R. B., 1983: A heat flux history length scale for the nocturnal boundary layer. Tulles, 35A, 

219-230.

Sykes, R. I., D. S. Henn, and S. F. Parker, 1996: SCIPUFF - A generalized hazard dispersion 

model. Preprints, 9th Joint Conference on Applications of Air Pollution Meteorology with 

AWMA, Amer. Meteor. Soc., American Meteorological Society, Boston, 184-188.

Taylor, G. I., 1921: Diffusion by continuous movements. Proc.London Math.Soc., Ser. 2(20), 

196-211.

Turner, D. B., T. Chico, and J. Catalano, 1986: TUPOS - A Multiple Source Gaussian Dispersion 

Algorithm Using On-site Turbulence Data. EPA/600/8-86/010, U.S. Environmental 

Protection Agency, RTP, NC, 39pp.

U.S. Environmental Protection Agency, 1992: Protocol for Determining the Best Performing 

Model. EPA-454/R-92-025, U.S. Environmental Protection Agency, RTP, NC. 

U.S. Environmental Protection Agency, 1995a: User's Guide for the Industrial Source Complex 

(ISC3) Dispersion Models (revised) Volume I - User Instructions. EPA-454/b-95-003a, 

U.S. Environmental Protection Agency, Research Triangle Park, NC.

U.S. Environmental Protection Agency, 1995b: Modeling Fugitive Dust Impacts from Surface 

Coal Mining Operations – Phase III. U.S. Environmental Protection Agency, Office of 

Air Quality Planning and Standards, EPA-454/R-96-002. December 1995.

U.S. Environmental Protection Agency, 2003: AERMOD: Latest Features and Evaluation 

Results. EPA-454/R-03-003, U.S. Environmental Protection Agency, RTP, NC. 
13

U.S. Enviromental Protection Agency, 2002: Compendium of Reports from the Peer Review 

Process for AERMOD. U.S. Environmental Protection Agency, RTP, NC. 

http://www.epa.gov/scram001/7thconf/aermod/dockrpt.pdf.

U.S. Environmental Protection Agency, 2015a: Technical Support Document (TSD) for NO2-

related AERMOD modifications. EPA document EPA-454/B-15-004. July 2015.

U.S. Environmental Protection Agency, 2015b: Addendum: User’s Guide for the AMS/EPA 

Regulatory Model – AERMOD. EPA-454/B-03-001. June 2015.

USGS, 1994: The 1994 plan for the National Spatial Data Infrastructure - Building the 

foundation of an information based society. Federal Geographic Data Committee Report. 

U.S. Geological Survey, Reston, VA. https://www.fgdc.gov/policyandplanning/NSDI

Strategy 1994.pdf.

van Ulden, A. P. and A. A. M. Holtslag, 1983: The stability of the atmospheric surface layer 

during nightime. Preprints, Sixth Symposium on Turbulence and Diffusion, American 

Meteorological Society, Boston, 257-260.

van Ulden, A. P. and A. A. M. Holtslag, 1985: Estimation of atmospheric boundary layer 

parameters for diffusion applications. J.Climate Appl.Meteor., 24, 1196-1207.

Venkatram, A., 1978: Estimating the convective velocity scale for diffusion applications. 

Bound.Layer Meteor., 15, 447-452.

Venkatram, A., 1980: Estimating the Monin-Obukhov length in the stable boundary layer for 

dispersion calculations. Bound.Layer Meteor., 19, 481-485.

Venkatram, A., 1982: A semi-empirical method to compute concentration associated with 

surface releases in the stable boundary layer. Atmos.Environ., 16, 245-248.

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

Venkatram, A., 1988: Dispersion in the stable boundary layer. Lectures on Air Pollution 

Modeling. Venkatram, A. and Wyngaard, J. C., American Meteorological Society, 229-

265pp.

Venkatram, A., 1992: Vertical dispersion of ground-level releases in the surface boundary layer. 

Atmos.Environ., 26A, 947-949.

Venkatram, A., R. Brode, A. Cimorelli, R. Lee, R. Paine, S. Perry, W. Peters, J. Weil, and R. 

Wilson, 2001: A complex terrain dispersion model for regulatory applications. 

Atmospheric Environment, 35, 4211-4221.

Venkatram, A., D. G. Strimaitis, and D. Dicristofaro, 1984: A semiemperical model to estimate 

vertical dispersion of elevated releases in the stable boundary layer. Atmos.Environ., 18, 

923-928.

Venkatram, A. and J. C. Wyngaard, 1988: Lectures on Air Pollution Modeling. American 

Meteorological Society, Boston, 390pp.

Weil, J. C., 1985: Updating applied diffusion models. J.Climate Appl.Meteor., 24(11), 1111-

1130.

Weil, J. C., 1988a: Dispersion in the convective boundary layer. Lectures on Air Pollution 

Modeling. Venkatram, A. and Wyngaard, J. C., American Meteorological Society, 167-

227pp.

Weil, J. C., 1988b: Plume rise. Lectures in Air Pollution Modeling. Venkatram, A. and 

Wyngaard, J. C., American Meteorological Society, 119-162pp.

Weil, J. C., 1992: Updating the ISC model through AERMIC. Preprints, 85th Annual Meeting of 

Air and Waste Management Association, Air and Waste Management Association, 

Pittsburgh, PA. 
15

Weil, J. C., 1996: A new dispersion model for stack sources in building wakes. Preprints, Ninth 

Joint Conference on Applications of Air Pollution Meteorology with the Air & Waste 

Management Association, American Meteorological Society, Boston, MA, 333-337.

Weil, J.C., 1998. The SERDP Open Burn/Open Detonation Dispersion Model (SOBODM), 

Volume IIa – Technical Description, and Volume IIb – Meteorological Inputs (DRAFT), 

CIRES, University of Colorado, Boulder, CO.

Weil, J. C. and R. P. Brower, 1983: Estimating convective boundary layer parameters for 

diffusion applications. PPSP-MD-48, Maryland Power Plant Siting Program, Maryland 

Department of Natural Resources, Baltimore, MD, 45pp.

Weil, J. C. and R. P. Brower, 1984: An updated Gaussian plume model for tall stacks. J.Air 

Poll.Control Assoc., 34, 818-827.

Weil, J.C., B. Templeman, R. Banta, R. Weber, and W. Mitchell, 1996. “Dispersion model 

development for open burn/open detonation sources,” Preprints 9th Joint Conference on 

Applications of Air Pollution Meteorology, American Meteorological Society, Boston, MA, 

610-616.

Weil, J. C., L. A. Corio, and R. P. Brower, 1997: A PDF dispersion model for buoyant plumes in 

the convective boundary layer. J.Appl.Meteor., 36, 982-1003.

Willis, G. E. and J. W. Deardorff, 1981: A laboratory study of dispersion in the middle of the 

convectively mixed layer. Atmos.Environ., 15, 109-117.

Wyngaard, J. C., 1988: Structure of the PBL. Lectures on Air Pollution Modeling. 

Venkatram, A., and Wyngaard, J. C., eds., American Meteorological Society, 9-57pp.

Zilitinkevich, S. S., 1972: On the determination of the height of the Ekman boundary layer. 

Bound.Layer Meteor., 3, 141-145.
United States

Environmental Protection

Agency

Office of Air Quality Planning and Standards

Air Quality Assessment Division

Research Triangle Park, NC

Publication No. EPA-454/R-18-003

April, 2018