AERMOD Model Formulation
                                       EPA-454/B-26-003
                                               July 2026




    AERMOD Model Formulation




U.S. Environmental Protection Agency
   Office of State Air Partnerships
  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 document has since
undergone multiple revisions, including a title change when the evaluation results were added,
and the document was retitled AERMOD Model Formulation and Evaluation. With the release
of AERMOD version 22112, the evaluation results were moved to a separate document,
AERMOD Model Evaluation, and the title of this document also changed to AERMOD Model
Formulation. This version of the document reflects additions and changes to the AERMOD
Modeling System released as version 26135.


         This report has been reviewed by the Office of State Air Partnerships, 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
Table of Contents


Section                                                                                                                                         Page
Notice ............................................................................................................................................. iii

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

Figures........................................................................................................................................... vii

1. Introduction ................................................................................................................................. 8

     1.1 Background ......................................................................................................................... 8

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

     1.3 Model development process ............................................................................................. 12

     1.4 Purpose of the document ................................................................................................... 14

2. Model overview ........................................................................................................................ 15

3. Meteorological preprocessor (AERMET)................................................................................. 20

     3.1 Energy balance in the PBL ............................................................................................... 21

           3.1.1 Net radiation............................................................................................................. 21

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

     3.2 Derived parameters in the CBL ........................................................................................ 23

           3.2.1 Friction velocity (u*) and Monin Obukhov length (L) in the CBL .......................... 23

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

     3.3 Derived parameters in the SBL ......................................................................................... 26

           3.3.1 Original friction velocity (u*) in the SBL ................................................................ 26

           3.3.2 Adjusted friction velocity (u*) in the SBL ............................................................... 29

           3.3.3 Sensible heat flux (H) in the SBL ............................................................................ 30

           3.3.4 Monin-Obukhov length (L) in the SBL ................................................................... 31

     3.4 Mixing height .................................................................................................................... 31

           3.4.1 Convective mixing height (zic)................................................................................. 31
                                                                                iv
         3.4.2 Mechanical mixing height (zim)................................................................................ 32

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

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

    4.1 General profiling equations............................................................................................... 37

         4.1.1 Wind speed profiling................................................................................................ 37

         4.1.2 Wind direction profiles ............................................................................................ 39

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

         4.1.4 Potential temperature profiling ................................................................................ 42

         4.1.5 Vertical turbulence calculated.................................................................................. 43

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

    4.2 Vertical inhomogeneity in the boundary layer as treated by the interface ....................... 52

5. The AMS/EPA Regulatory Model: AERMOD ........................................................................ 58

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

    5.2 Concentration predictions in the CBL .............................................................................. 66

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

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

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

    5.3 Concentrations in the SBL ................................................................................................ 79

    5.4 Treatment of lateral plume meander ................................................................................. 80

    5.5 Estimation of dispersion coefficients ................................................................................ 83

         5.5.1 Dispersion from ambient turbulence ........................................................................ 85

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

         5.5.3 Treatment of building downwash ............................................................................ 92

    5.6 Plume rise calculations in AERMOD ............................................................................... 93

         5.6.1 Plume rise in the CBL .............................................................................................. 93


                                                                      v
           5.6.2 Plume rise in the SBL .............................................................................................. 95

     5.7 Source characterization ..................................................................................................... 97

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

           5.8.1 Definition of plume volume ..................................................................................... 99

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

     5.9 Generic Reaction Set Method (GRSM) .......................................................................... 105

     5.10 Adjustments for the urban boundary layer.................................................................... 107

6. List of symbols........................................................................................................................ 113

Appendix A. Input / output needs and data usage ...................................................................... 119

Appendix B. Description of ALPHA Options in AERMOD 26135........................................... 129

References ................................................................................................................................... 164




                                                                             vi
Figures


Figure                                                                                                                                        Page

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

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

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

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

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

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

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

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

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

Figure 10. AERMOD two state approach. The total concentration predicted by AERMOD is the
weighted sum of the two extreme possible plume states. ............................................................. 61

Figure 11. Treatment of Terrain in AERMOD. Construction of the weighting factor used in
calculating total concentration. ..................................................................................................... 64

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

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

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

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. ............................................................ 73

Figure 16. Lateral spread (Fy) as a function of non-dimensional distance (X). The data is taken
from the Prairie Grass experiment (Barad, 1958). ........................................................................ 87

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


                                                                              vii
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 1970s 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-1980s (e.g., see Briggs
(1988) and Venkatram (1988)).


       During the mid-1980s, 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

                                                      8
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-1980s, new applied dispersion models based on this technology had been
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 1980s, a substantial scientific base on the PBL and new dispersion
approaches existed for revamping regulatory dispersion models. 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,


                                                    9
other problems such as plume interaction with terrain, surface releases, building downwash and
urban dispersion were recognized as needing attention.


       The model developed by AERMIC is aimed at short-range dispersion from stationary
industrial sources, the same scenario handled by its predecessor, the EPA Industrial Source
Complex Model, ISC3 (U.S. EPA, 1995). This work clearly has benefitted from the model
development activities of the 1980s, especially in the parameterization of mean winds and PBL
turbulence, dispersion in the CBL, and the treatment of plume/terrain interactions. Techniques
used in this model for PBL parameterizations and CBL dispersion are like 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 was on the regulatory models designed for estimating near-field
impacts from a variety of industrial source types. EPA’s regulatory platform for near-field
modeling, with few exceptions, has remained fundamentally unchanged. 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) ensuring 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.
                                                   10
         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 pre-
processor, AERMAP (U.S. EPA, 2024c), both characterizes the terrain and generates receptor
grids for the dispersion model (AERMOD).


         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 (e.g., airport), on-site (site-specific), or forecast data from a
prognostic weather model, 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. AERMAP can process multiple
formats of gridded elevation data 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 determined 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 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


                                                      11
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 AERMOD, a regulatory replacement for the ISC3 model. The process
followed was: 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 what was then EPA’s Office of Air
Quality Planning and Standards (OAQPS) for consideration as a regulatory model.


       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 to identify areas needing improvement. The
developmental evaluation provided a basis for selecting formulation options.


       This developmental evaluation was conducted using five databases. Three consisted of
event-based tracer releases, while the other two each contain up to a full year of continuous SO2
measurements. These databases 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. EPA, 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 that was repeated
after revisions based on the developmental evaluation and peer review.




                                                     12
       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 databases 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 databases (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. EPA, 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. EPA, 1992).


       On April 21, 2000, EPA proposed 1 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 May 19, 2000, EPA announced 2 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 Plume Rise
Model Enhancements (PRIME) algorithms for building downwash (Schulman et al., 2000), 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




       1
           40 CFR Part 51 pages 21506-21546
       2
           Federal Register on May 19, 2000 (Volume 65, Number 98)
                                                    13
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 details that are not included in the published journal articles (Cimorelli et al. 2004; Perry
et al. 2003).




                                                    14
2. Model overview


       This section provides a general overview of the most important features of AERMOD.
AERMOD serves as a complete replacement for ISC3. Thus, the AERMOD model described
here is applicable to rural and urban areas, flat and complex terrain, limited overwater
applications, surface and elevated releases, and multiple source types (including, point, area,
volume, buoyant line, and mobile sources). Every effort has been made to avoid model
formulation discontinuities wherein large changes in calculated concentrations result from small
changes in input parameters.


       Since it was originally developed and promulgated as a preferred EPA model, the
AERMOD Modeling System has been updated with added capabilities including improved
characterization of certain emission sources with new source types and the ability to process
offshore meteorology.


       AERMET was originally designed to process land-based meteorology (i.e.,
meteorological data collected over land surfaces) to estimate PBL parameters influenced by
land-air interactions where due to the diurnal cycle of surface heating and cooling, the
atmosphere is generally convective during daytime hours and stable during the nighttime.
Beginning with version 23132, AERMET was updated to include the Coupled Ocean
Atmosphere Response Experiment (COARE) (Fairall, et al., 2003) to appropriately characterize
the marine boundary layer using data collected offshore. This addition to AERMET is the first
in a series of planned enhancements to the AERMOD Modeling system to expand AERMOD’s
capability to model offshore emission sources. For information on the formulation and usage of
the COARE algorithm in AERMET, the reader is directed to Fairall, et al. (2003). The
formulation of AERMET as described in this document continues to focus on the treatment of
land-based meteorology.




                                                    15
       The BUOYLINE source type for characterizing buoyant line sources was added
beginning with version 15181. The BUOYLINE source type serves as an implementation of the
Buoyant Line and Point (BLP) source dispersion model originally developed to estimate
impacts from aluminum smelter potline emissions. Integration and testing of the BLP model and
source type with AERMOD is described in EPA (2016). While a few functional enhancements
have been made to the BUOYLINE source type beyond the capability of the BLP model, the
formulation of the BUOYLINE source type has not been updated from the formulation
described in the Buoyant Line and Point Source Dispersion Model User’s Guide (BLP User’s
Guide) (Schulman and Scire, 1980). Users should refer to the BLP User’s Guide for a
description of the formulation of the BUOYLINE source type in AERMOD.


       Similarly, the RLINE source type was added beginning with version 19191 to better
characterize mobile source emissions (i.e., roadways). The RLINE source type is an adaptation
of the Research Line (R-LINE) model (Snyder and Heist, 2013). For a description of the
formulation of the RLINE source type in AERMOD, refer to Snyder and Heist (2013) and EPA
(2023e).


       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 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.


                                                     16
        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 stable
stratified conditions.


        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. AERMOD can operate using data that are 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. 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 to calculate the convective mixing height throughout the day. Surface
characteristics (surface roughness, Bowen ratio, and albedo) are also needed to construct
similarity profiles of the relevant PBL parameters.


        Unlike prior 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.
                                                      17
       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.




                   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
                                                  18
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 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
Elevation Model (DEM) data. AERMAP can also be used to create receptor grids. The
elevation for each specified receptor can be 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 26135 version of AERMOD, AERMET, and AERMAP. In addition, all the symbols used
for the many parameters and variables that are referred to in this document are defined, with
their appropriate units, in Section 6 titled “List of Symbols.”




                                                     19
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 4.


       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 like 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 can estimate 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 necessary for
AERMOD.


       It should be noted as mentioned in section 1.2 of the Introduction to this document, the
formulation descriptions that follow focus on the treatment of land-based observed data and the
characterization of the planetary boundary layer over land. For information on the formulation


                                                     20
of the COARE algorithm and the treatment of offshore meteorology to characterize the marine
boundary layer, refer to Fairall, et al. (2003).


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
calculate the sensible heat flux, H. We begin with the following simple characterization of the
energy balance in the PBL:


                                         𝐻𝐻 + 𝜆𝜆 𝐸𝐸 + 𝐺𝐺 = 𝑅𝑅𝑛𝑛                             (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


                                                     0.9 𝑅𝑅𝑛𝑛                               (2)
                                          𝐻𝐻 =                   .
                                                   �1 + 1�𝐵𝐵 �
                                                              𝑜𝑜


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



                                                          21
                                                           6
                           𝑅𝑅𝑛𝑛 = (1 − 𝑟𝑟{𝜑𝜑}) 𝑅𝑅 + 𝑐𝑐1 𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟             4
                                                                 − 𝜎𝜎𝑆𝑆𝑆𝑆 𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟 − 𝑐𝑐2 𝑛𝑛     (3)
                                                                                            ,
                                                  1 + 𝑐𝑐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


        {𝜑𝜑} = 𝑟𝑟 ′ + (1 − 𝑟𝑟 ′ ) exp(𝑎𝑎𝑎𝑎 + 𝑏𝑏),


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


                                           𝑅𝑅 = 𝑅𝑅𝑜𝑜 (1 − 0.75𝑛𝑛3.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 𝜑𝜑 = ((𝜑𝜑{𝑡𝑡_𝑝𝑝 } + 𝜑𝜑{𝑡𝑡})) ⁄ (2)) is the solar elevation (tp and t are the
previous and present hours, respectively). 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



                                                               22
                                             1    −𝑐𝑐1 𝑇𝑇 6 + 𝜎𝜎𝑆𝑆𝑆𝑆 𝑇𝑇 4 − 𝑐𝑐2 𝑛𝑛          (5)
                       sin(𝜑𝜑𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 ) =      �                                   + 30�.
                                            990 (1 − 𝑟𝑟{𝜑𝜑})(1 − 0.75𝑛𝑛3.4 )


       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

                                                                      1⁄3.4
                                                        𝑅𝑅 ⁄𝑅𝑅𝑜𝑜�
                                       𝑛𝑛𝑒𝑒𝑒𝑒 = �1 −             0.75�


3.2 Derived parameters in the CBL


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


       In the CBL, AERMET computes the surface friction velocity, u*, and the Monin-

Obukhov 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, like 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,




                                                             23
                                                            𝑘𝑘 𝑢𝑢𝑟𝑟𝑟𝑟𝑟𝑟
                        𝑢𝑢∗ =                                                                    ,   (6)
                                 ln �𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟 ⁄𝑧𝑧0 − 𝛹𝛹𝑚𝑚 �𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟 ⁄(𝐿𝐿} + 𝛹𝛹𝑚𝑚 {𝑧𝑧0 ⁄(𝐿𝐿})��


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:


                              𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟            1+µ         1 + µ2                 𝜋𝜋
                     𝛹𝛹𝑚𝑚 �            � = 2 ln �     � + ln �       � − 2 tan−1 µ +
                                 𝐿𝐿                2            2                    2
                                                                                                     (7)
                           𝑧𝑧0        1 + µ0          1 + µ20                  𝜋𝜋
                     𝛹𝛹𝑚𝑚 � � = 2 ln �       � + ln �         � − 2 tan−1 µ0 +
                            𝐿𝐿          2               2                      2


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)):


                                                       𝜌𝜌 𝑐𝑐𝑝𝑝 𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟 𝑢𝑢∗3
                                                𝐿𝐿 = −                                               (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 equations (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
                                                                 24
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 where 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 7zo (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, to estimate turbulence in the CBL, an
estimate of w* is needed. AERMET calculates the convective velocity scale from its definition
as:

                                                               1⁄3
                                               𝑔𝑔 𝐻𝐻 𝑧𝑧𝑖𝑖𝑖𝑖
                                       𝑤𝑤∗ = �                 �     ,                 (9)
                                              𝜌𝜌 𝑐𝑐𝑝𝑝 𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟


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



                                                          25
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 Original friction velocity (u ) in the SBL
                                  *




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


                                             𝜃𝜃∗ = −𝐻𝐻⁄ρ cp 𝑢𝑢∗ ,                        (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


                                                       𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟 2
                                               𝐿𝐿 =             𝑢𝑢
                                                      𝑘𝑘 𝑔𝑔 𝜃𝜃∗ ∗                        (11)


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


                                             𝑢𝑢∗       𝑧𝑧  𝛽𝛽𝑚𝑚 𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟
                                      𝑢𝑢 =       �ln � � +               �,              (12)
                                              𝑘𝑘      𝑧𝑧0      𝐿𝐿

                                                            26
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


                                             𝑢𝑢   1    𝛽𝛽𝑚𝑚 𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟 𝑔𝑔 𝜃𝜃∗
                                                =    +                      .                            (13)
                                            𝑢𝑢∗ 𝐶𝐶𝐷𝐷      𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟 𝑢𝑢∗2


         Multiplying eq. (13) by CD u*2 and rearranging yields a quadratic of the form


                                           𝑢𝑢∗2 − 𝐶𝐶𝐷𝐷 𝑢𝑢 𝑢𝑢∗ + 𝐶𝐶𝐷𝐷 𝑢𝑢02 = 0                            (14)


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

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

                                                                                          1
                                                     ⎡                          2 2⎤

                                       𝐶𝐶𝐷𝐷 𝑢𝑢𝑟𝑟𝑟𝑟𝑟𝑟 ⎢              2 𝑢𝑢0          ⎥
                                 𝑢𝑢∗ =               ⎢ 1 + ⎛1 − � 1            � ⎞ ⎥                     (15)
                                            2
                                                     ⎢
                                                           ⎝     𝐶𝐶𝐷𝐷2 𝑢𝑢𝑟𝑟𝑟𝑟𝑟𝑟 ⎠ ⎥
                                                     ⎣                             ⎦


         Equation (15) produces real-valued solutions only when the wind speed is greater than or
                                                                          1
equal to the critical value 𝑢𝑢𝑐𝑐𝑐𝑐 = �4 𝛽𝛽𝑚𝑚 𝑧𝑧𝑟𝑟𝑟𝑟𝑟𝑟 𝑔𝑔 𝜃𝜃∗ ⁄𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟 𝐶𝐶𝐷𝐷 �2 . For the wind speed less than the
critical value, u* and θ* are parameterized using the following linear expression:


                                                                𝑢𝑢
                                  𝑢𝑢∗ = 𝑢𝑢∗ {𝑢𝑢 = 𝑢𝑢𝑐𝑐𝑐𝑐 } �          �   𝑓𝑓𝑓𝑓𝑓𝑓 𝑢𝑢 < 𝑢𝑢𝑐𝑐𝑐𝑐
                                                               𝑢𝑢𝑐𝑐𝑐𝑐
                                              𝑢𝑢
                                  𝜃𝜃∗ = 𝜃𝜃∗ � �                           𝑓𝑓𝑓𝑓𝑓𝑓 𝑢𝑢 < 𝑢𝑢𝑐𝑐𝑐𝑐
                                             𝑢𝑢𝑐𝑐𝑐𝑐


         These expressions approximate the u* versus θ* dependence found by van Ulden and
Holtslag (1983).

                                                                    27
       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 𝑛𝑛2 ) ,                         (16)


where n is the fractional cloud cover.


       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,


                                                   𝜃𝜃∗ ln 𝑧𝑧       𝑧𝑧
                                  𝜃𝜃{𝑧𝑧} − 𝜃𝜃0 =      �      + 𝛽𝛽𝑚𝑚 � ,                  (17)
                                                    𝑘𝑘 𝑧𝑧0         𝐿𝐿


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


                                                  𝑘𝑘(𝜃𝜃2 − 𝜃𝜃1 )
                                  𝜃𝜃∗ =                               .
                                            ln 𝑧𝑧          (𝑧𝑧 − 𝑧𝑧 )                    (18)
                                          �� 𝑧𝑧 2 � + 𝛽𝛽𝑚𝑚 2 𝐿𝐿 1 �
                                               1



       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:

                                                        28
                                                            ku
                         u* =                                                     .
                                 z          z 4.25       0.5    βm                   (19)
                                 ln  + 7 ln  +       −        +    −  .
                                                                      1648 
                                 zo        L ( z L ) ( z L )2 2          

3.3.2 Adjusted friction velocity (u*) in the SBL


          Beginning with the 2017 updated regulatory version of the model (version 16216), a
second approach for calculating the u* value based on Qian and Venkatram (2011) was added.
Under this formulation, u* is computed as a function of the drag coefficient for neutral conditions
(CDN),
                                                        𝑘𝑘 2
                                      𝐶𝐶𝐷𝐷𝐷𝐷 =
                                                      𝑧𝑧 − 𝑑𝑑 2
                                                 �ln � r 𝑧𝑧 ℎ ��
                                                            0



where zr is the measurement height and dh the zero-plane displacement (see Eq. 20 in Qian and
Venkatram (2011)). Using CDN, u* can be calculated based on eq. 26 from Qian and Venkatram
(2011),
                                                            𝑟𝑟 2
                                        𝐶𝐶DN 𝑈𝑈  1 + exp �−      �
                                                             2
                                  𝑢𝑢∗ =         �                 �,
                                           2                 2
                                                  1 − exp �− 𝑟𝑟 �

with
                                                     2𝑢𝑢0
                                            𝑟𝑟 =                1
                                                   𝑈𝑈(𝐶𝐶DN )4
and
                                                                    1
                                          𝛽𝛽g (𝑧𝑧r − 𝑑𝑑h − 𝑧𝑧0 )𝑇𝑇0 2
                                 𝑢𝑢0 = ��                          � �
                                                    𝑇𝑇0


          The minimum limit of the u* value here can be computed as
                                                        29
                        1
                        2
               𝑈𝑈𝐶𝐶𝐶𝐶 𝐶𝐶𝐷𝐷𝐷𝐷
𝑢𝑢∗,𝑚𝑚𝑚𝑚𝑚𝑚 =                   th
                    4
                                                                           2𝑢𝑢0
                                                               𝑈𝑈𝐶𝐶𝐶𝐶 =         1    .
                                                                             4
                                                                           𝐶𝐶𝐷𝐷𝐷𝐷



          The final u* value is then the maximum of u*,min and the u*. Using this value of u*, θ* is
computed based on eq. (18).


          As with the original method, 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; however, when using the adjusted u*
formulation, the computation is based on Luhar and Rayner (2009). When the zref (reference
height for the wind speed) is greater than 0.7*L (the Monin-Obukhov Length), then u* is
calculated based on the initial estimate of L derived from the two-level temperature
measurements,
                                                                    𝑘𝑘𝑘𝑘
                                    𝑢𝑢∗ =                                       𝛽𝛽
                                            𝛼𝛼�𝜁𝜁𝛽𝛽 (1 + 𝛾𝛾𝜁𝜁 (1−𝛽𝛽) ) − 𝜁𝜁0 �1 + 𝛾𝛾𝜁𝜁0 (1−𝛽𝛽) ��

where α = 4, β = 0.5and γ = 0.3 and
                                                                           𝑧𝑧
                                                                  𝜁𝜁 = 𝐿𝐿0.

Using this value of u*, θ* is computed based on eq. (18).

3.3.3 Sensible heat flux (H) in the SBL


          Having computed u* and θ* for stable conditions, AERMET calculates the surface heat
flux from eq. (10) as


                                                       𝐻𝐻 = −𝜌𝜌 𝑐𝑐𝑝𝑝 𝑢𝑢∗ 𝜃𝜃∗ .                      (20)



                                                                     30
        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


                                     [𝜃𝜃∗ 𝑢𝑢∗ ]𝑚𝑚𝑚𝑚𝑚𝑚 = 0.05𝑚𝑚𝑠𝑠 −1 𝐾𝐾 .                   (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.4 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.


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

                                                         31
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

                                            𝑧𝑧𝑖𝑖𝑖𝑖                          𝑡𝑡                (22)
                                                                            𝐻𝐻{𝑡𝑡 ′ }
                            𝑧𝑧𝑖𝑖𝑖𝑖 𝜃𝜃{𝑧𝑧𝑖𝑖𝑖𝑖 } − � 𝜃𝜃{𝑧𝑧}𝑑𝑑𝑑𝑑 = (1 + 2𝐴𝐴) �           𝑑𝑑𝑑𝑑,
                                                                             𝜌𝜌 𝑐𝑐𝑝𝑝
                                            0                              0



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 ensure 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


                                                           𝑢𝑢 𝐿𝐿 0.5                          (23)
                                             𝑧𝑧𝑖𝑖𝑖𝑖 = 0.4 � ∗ �𝑓𝑓� ,


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


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

                                                              32
                                                                       3
                                                   𝑧𝑧𝑖𝑖𝑖𝑖 = 2400 𝑢𝑢∗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)


                                               𝑑𝑑𝑧𝑧𝑖𝑖𝑖𝑖 (𝑧𝑧𝑖𝑖𝑖𝑖 − 𝑧𝑧𝑖𝑖𝑖𝑖 )
                                                       =                   .                                  (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


                  𝑧𝑧𝑖𝑖𝑖𝑖 {𝑡𝑡 + 𝛥𝛥𝛥𝛥} = 𝑧𝑧𝑖𝑖𝑖𝑖 {𝑡𝑡}�𝑒𝑒 −𝛥𝛥𝛥𝛥⁄𝜏𝜏 � + 𝑧𝑧𝑖𝑖𝑖𝑖 {𝑡𝑡 + 𝛥𝛥𝛥𝛥}�1 − �𝑒𝑒 −𝛥𝛥𝛥𝛥⁄𝜏𝜏 �� .   (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


                                                                  33
scale (and thus relaxation time) will be short if the equilibrium mixing height grows rapidly but
will be long if it decreases rapidly.


        Although equations (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 nighttime to daytime 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:


                        zi = MAX[zic ; zim]   for L < 0 (CBL)
                                                                                       (27)
                        zi = zim              for L > 0 (SBL)


        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). 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



                                                    34
AERMET control file for versions of AERMET prior to 22112 and the Stage 2 control file
beginning with version 22112.


       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.




                                                    35
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 A.1). 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 sub-sections in Section 4.1 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).




                                                         36
4.1 General profiling equations


4.1.1 Wind speed profiling


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


                                 z 
                    u = u{7 zo }                                    for   z < 7 zo
                                 7 zo 
                         u*   z         z       zo                                   (28)
                    u=       ln  − Ψ m   + Ψ m                for 7 zo ≤ z ≤ zi
                         k   zo         L       L 
                    u = u{zi }                                        for z > zi

       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 Ψms are evaluated using eq. (7) with zref replaced by z, and during
stable conditions they are calculated from van Ulden & Holtslag (1985) as


                              z                     z
                         Ψ m   = − 17 1 − exp − 0.29  
                              L                     L 
                                                                                          (29)
                              zo                    zo  
                         Ψ m   = − 17 1 − exp − 0.29   .
                              L                      L 

       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.

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




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



4.1.2 Wind direction profiles


       For both the CBL and 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.




                                                 39
4.1.3 Profiles of the potential temperature gradient


       Above the relatively shallow super adiabatic 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
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 ensure 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,




                                                       40
                      ∂θ   θ      (2) 
                         = * 1 + 5                     for z ≤ 2m
                      ∂z k (2)     L 
                                                                                       (30)
                      ∂θ θ*      z
                        =    1+ 5                  for 2m < z ≤ 100m.
                      ∂z k z     L


       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


                               ∂θ ∂θ {z mx }      ( z − z mx ) 
                                  =          exp −                                   (31)
                               ∂z    ∂z           0.44 ziθ 


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 below 100 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 equations (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) θ* =
                                                    41
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.




                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



                                                      42
                                                               g z msl
                                      θ {zTref } = Tref +              ,                    (32)
                                                                cP


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 + ∆z} = θ {z} +              ∆z                     (33)
                                                                ∂z z


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

                                                                    2
        In the CBL, the vertical velocity variance or turbulence (𝜎𝜎𝑤𝑤𝑤𝑤 ) 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:


                                            2
                                          𝜎𝜎𝑤𝑤𝑤𝑤     2
                                                 = 𝜎𝜎𝑤𝑤𝑤𝑤     2
                                                          + 𝜎𝜎𝑤𝑤𝑤𝑤 .                        (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 (σ 2wc) of the total variance is calculated as:




                                                          43
                                          23
                        2        z
                    σ   wc = .
                            16               ⋅ w*2                        . z ic
                                                                   for z ≤ 01
                                z ic 
                      2
                    σ wc = 0.35w*2                                      . z ic < z ≤ z ic
                                                                   for 01                   (35)

                      2                         6( z − z ic ) 
                    σ wc = 0.35w*2 exp −                          for z > z ic
                                                    z ic 


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 σ 2wc 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.




                                                            44
              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,


                          2      2       2
                        σ wm = σ wml + σ wmr .                                         (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 σ 2wm with z
near z = zi. The expression for σwml following the form of Brost et al. (1982) is


                                                    45
                                                      12                                 (37)
                                              z
                               σ wml = .
                                      13u* 1 −           for z < zi
                                            zi 
                               σ wml = 0.0                 for z ≥ zi


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.




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




                                   47
       In the SBL the vertical turbulence contains only a mechanical portion which is given by
eq. (36). The use of the same σ 2wm 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.




                      Figure 7. Profile of vertical turbulence in the SBL

4.1.6 Lateral turbulence calculated by the interface

                                                   2
       In the CBL the total lateral turbulence, 𝜎𝜎𝑣𝑣𝑣𝑣 , is computed as a combination of a
mechanical (σvm) and convective (σvc) portions such that

                                                       48
                                           2      2      2
                                         σ vT = σ vc + σ vm .                           (38)


       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 surface 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


                                                σ v 0 2 = Cu*2 ,                        (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, σ2vm is assumed to
vary linearly as



                       σ   2
                                  2
                               σ vm {zim } − σ vo2  z + σ 2        for z ≤ zim
                           vm =                           vo
                                     zim                                            (40)
                         2
                       σ vm     2
                            = σ vm {zim }                          for z > zim ,


                                                      49
         2 {𝑧𝑧 }          2            2 −2         2
where 𝜎𝜎𝑣𝑣𝑣𝑣  𝑚𝑚 = min[𝜎𝜎𝑣𝑣𝑣𝑣 ; 0.25 𝑚𝑚 𝑠𝑠 ] and σvo , the surface value of the lateral turbulence, is

equal to 3.6 u*2. This linear variation of σ 2vm 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).


        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 s-1
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.




                                                      50
   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)
                                         σ vc2 = 0.35w*2


       This constant value of 𝜎𝜎𝑣𝑣2 ⁄𝑤𝑤∗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

                                                       51
1979).


         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 layer 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 𝑢𝑢�).


         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

                                                       52
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 9. 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.




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



        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 𝜎𝜎�𝑤𝑤𝑤𝑤 {𝑥𝑥𝑟𝑟 } and 𝑢𝑢�{𝑥𝑥𝑟𝑟 } are

calculated. Once the averaging layer for a given plume and receptor is established the effective
values, 𝛼𝛼�, are computed as simple averages:


                                                               ht
                                                      1
                                            α~ =          ∫ α {z} dz                             (42)
                                                   ( t b) h
                                                    h − h       b
                                                                           ,

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



                H p {xr , yr },
                                                           if H p {xr , yr } ≤ zr    
          hb =                                                                        
               
                        {[                                    ] }
                                      215σ z {x sr } , zr , if H p {xr , yr } > zr 
                 MAX H p {xr , yr } − .
                                                                                       
                                                                                                    (43)
                        {[                                    ] }
                                     215σ z {x sr } , zr , if H p {xr , yr } ≤ zr 
                MIN H p {xr , yr } + .
               
          ht = 
                                                                                       
                                                                                       
               
                H p { x r , y r },                        if  H p { x r , y r } > z r 
                                                                                       


        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 𝛼𝛼� = 𝛼𝛼{𝑧𝑧𝑖𝑖 }.


        For plumes in stable conditions and for the penetrated source in the CBL, Hp is always
set equal to the plume centerline height (Δℎ𝑠𝑠 + ℎ𝑠𝑠 ) 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 5.6.2.


        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.


        Prior to plume stabilization, i.e., x < xf (distance to plume stabilization),


                                            𝐻𝐻𝑝𝑝 = ℎ𝑠𝑠 + 𝛥𝛥ℎ𝑑𝑑,𝑝𝑝 ,


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).

                                                           55
        The distance to plume stabilization, xf , is determined following Briggs (Briggs 1975;
Briggs 1971) as


                                      x f = 49 Fb5 8                  for Fb < 55
                                                                                               (44)
                                      x f = 119 Fb2 5                 for Fb ≥ 55
                                                                                           ,


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. EPA, 1995)
expression



                                                       (              )
                                                                          2
                                                  8rs ws + 3u p
                                           xf =                                                (45)
                                                                              ws u p
                                                                                       ,


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.


        The distance xm is taken to be the product of the average mixed layer wind speed and the
                     𝑧𝑧
mixing time scale, 𝜎𝜎 𝑖𝑖 . That is,
                      𝑤𝑤𝑤𝑤




                                                           u zi
                                               xm =               ,                            (46)
                                                           σ wT

                                                             56
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:



                                                z                 (
                                                               x − xf
                               H p = H {x f } +  i − H {x f } ⋅
                                                                            )           (47)
                                                2                (
                                                               xm − x f     ).
        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.




                                                     57
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 urban-
rural temperature difference as suggested by Oke (1978; 1982).
                                                      58
       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 a combination of two limiting cases: a horizontal plume (terrain impacting) and a terrain-
following (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
                                                      59
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.


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
removing 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


                    CT {xr , yr , zr } = f ⋅ Cc , s {xr , yr , zr } + (1 − f ) Cc , s {xr , yr , z p }       (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),
                                                     60
Cc,s{zr,yr,zr} is the contribution from terrain-following state, f is the plume state weighting
function, { x r , y r , z 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.




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
                                                   61
receptor specific terrain height scale (hc) from AERMAP, Hc is calculated from the same
algorithms found in CTDMPLUS as:


                                                   hc
                                1 / 2 ⋅ u {H c } = ∫ N 2 ( hc − z ) dz.
                                        2
                                                                                               (49)
                                                   Hc



                                                                    1�
                                                          𝑔𝑔 𝜕𝜕𝜕𝜕     2
where u {Hc} is the wind speed at height Hc, and 𝑁𝑁 = �𝜃𝜃 𝜕𝜕𝜕𝜕 �          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.


       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 contrasts with 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
                                                     62
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:


                                                Hc
                                                 Cs {xr , yr , zr } dz
                                       ϕ p = ∫∞0                                          (50)
                                             ∫0 Cs {xr , yr , zr } dz ,

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.
                                                      63
       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.




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:


                             C{x , y , z} = (Q u~ ) Py { y; x} Pz {z; x},               (51)


where Q is the source emission rate, 𝑢𝑢� is the effective wind speed, and py and pz are probability
density functions (pdf) which describe the lateral and vertical concentration distributions,
                                                    64
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
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
                                                     65
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., 𝑢𝑢�).


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
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
non-Gaussian 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 affected 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 ensure 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


                                                      66
component, the lateral velocity pdf is approximately Gaussian (Lamb 1982), and this pdf and
the resulting concentration distribution, Fy, are assumed to be Gaussian.




     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
                                                     67
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
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 re-entrained by and disperses in the growing CBL.




                                                      68
                  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
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 no-
flux 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.



                                                      69
       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


                                                           wx
                                        zc = hs + ∆h +        ,                           (52)
                                                           u

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




              pw =
                    λ1
                              w− w 2
                          exp −
                                       (
                                      1 
                                          +
                                             λ2 )       w− w 2
                                                    exp −
                                                                2 
                                                                    ,
                                                                           (        )        (53)
                   2πσ w1       2σ w1 
                                    2
                                            2πσ w 2       2σ w 2 
                                                              2
                                                               


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
                                                            3                  3
turbulent velocity), the vertical velocity skewness 𝑆𝑆 = 𝑤𝑤 �𝜎𝜎𝑤𝑤3 (where 𝑤𝑤 is the third moment of
w), and a parameter 𝑅𝑅 = 𝜎𝜎𝑤𝑤1 ⁄𝑤𝑤1 = − 𝜎𝜎𝑤𝑤2 ⁄𝑤𝑤2 = 2. An expanded discussion of the PDF
parameters is given in Weil et al. (1997).




                                                      70
       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 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 ( w 1 ) or downdraft ( w 2 )
velocity. The mean height of each trajectory, 𝑧𝑧𝑐𝑐1 or 𝑧𝑧𝑐𝑐2 , can be found by averaging eq. (53).
The 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).




                                                      71
    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 bi-
Gaussian 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).




                                                     72
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
                                                      73
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
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 re-entrainment into the CBL. Figure 13 illustrates the 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


                  Cc {xr , yr , zr } = Cd {xr , yr , zr } + Cr {xr , yr , zr } + C p {xr , yr , zr },    (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


                     fp = 0                          if ∆hh < 0.5∆heq
                     fp =1                           if ∆hh > .
                                                             15∆heq
                                                                                                  (55)
                          ∆hh
                     fp =      − 0.5                 if 0.5∆heq < ∆hh > .
                                                                       15∆heq
                          ∆heq


where Δhh = zi - hs, and Δheq is the equilibrium plume rise in a stable environment. The Δheq has
the form Berkowicz et al. (1986)



                                                         74
                                                   (
                                         ∆heq = 2.63 Ps + (2 3)          ) ∆h ,
                                                                        3 13
                                                                                 h                     (56)



where 𝑃𝑃𝑠𝑠 = 𝐹𝐹𝑏𝑏 ⁄𝑢𝑢𝑁𝑁ℎ2 Δℎℎ3 is the penetration parameter, and the stack buoyancy flux (Fb), and
Brunt-Vaisala frequency (Nh) are given respectively by


                                                                    ∆T
                                                   Fb = gws rs2                                        (57)
                                                                    Ts ,


        and


                                                                        12
                                               g ∂θ          
                                         Nh =                 .                                      (58)
                                               { i }
                                                θ z    ∂z     
                                                          z>z      i



        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:



                          Qf p         2 ∞ λ              ( z − Ψ − 2mz )2         ( z + Ψ + 2mz )2  
       Cd {xr , yr , z} =       F  ⋅ ∑    ∑
                                                  j 
                                                      exp  −
                                                                   dj     i 
                                                                               + exp −       dj     i 
                                                                                                           ,   (59)
                          2π u~
                                 y
                                     j =1 m = 0σ zj                  2                        2     
                                                     
                                                                   2 σ zj                   2 σ zj     


        where
                                                             75
                                                                    wjx
                                      Ψdj = hs + ∆hd +                    ,                       (60)
                                                                     u

                                               1            −𝑦𝑦 2
u is the wind speed at stack top, 𝐹𝐹𝑦𝑦 �=              exp � 2𝜎𝜎2 �� the lateral distribution function with
                                            √2𝜋𝜋𝜎𝜎𝑦𝑦           𝑦𝑦

meander (discussed in Section 5.4), 𝑤𝑤𝑗𝑗 = 𝑎𝑎𝑗𝑗 𝑤𝑤∗ (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 terrain-
following states, respectively. Here, Ψdj and σzj are the effective source height and vertical
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 building-
induced turbulence are calculated as discussed in Sections 5.5.1.1 and 5.5.1.2 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 5.5.1.2) and λj, the weighting coefficient
for each distribution in eq. (53), is calculated from Weil et al. (1997) as


                                                 w2         a2
                                        λ1 = +         =+
                                               w2 − w1    a2 − a1
                                                                                                  (61)
                                                  w1         a1
                                        λ2 = −         =−
                                               w2 − w1    a2 − a1


where



                                        σ~wT  α S
                                                                               12
                                               1        4                         
                                  a1 =        + α 2S 2 +                          
                                       w*  2 2        β                         
                                                                                    
                                                                                                  (62)
                                        σ~
                                                             12
                                                 1  2 2 4  
                                        wT  α S
                                  a2 =           − α S + 
                                       w*  2 2         β  
                                                                 .

                                                             76
         Recall that σ~wT is the total effective vertical turbulence and is calculated from eq. (34).

The parameters appearing in eq. (62) are given by


                           w3
                               = .
                                0125                    for H p {x} ≥ 01
                                                                       . zi
                           w*3
                                                                                              (63)
                           w3        H p {x}
                               = 125
                                  .                     for H p {x} < 01
                                                                       . zi
                           w*3         zi


where,


                                    1 + R2
                               α=
                                   1 + 3R 2
                               β = 1 + R2                                                     (64)
                                      w 3 w*3
                               S=               ≡ Skewness factor ,
                                     ~
                                    ( wT          )
                                              3
                                     σ w      *



and R is assumed to be 2.0 (Weil et al., 1997). Likewise, the term 𝑤𝑤𝑗𝑗 𝑥𝑥⁄𝑢𝑢 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.




                                                        77
5.2.2 Indirect Source contribution to concentration calculations in the CBL


        The concentration due to the indirect source is calculated from:



                                               λ j   ( z + Ψrj − 2mzi )        ( z − Ψ + 2mz )2  
                                                                          2
                         Qf p        2   ∞
                                                                              + exp −
                                                                                            rj      i 
       Cr {xr , yr , z} =    ~ ⋅ Fy ⋅ ∑ ∑ exp −                                                         , (65)
                          2π u        j =1 m =1σ zj              2                           2    
                                                     
                                                               2 σ zj                     2  σ zj    


where Ψrj = Ψdj - Δhi, and z is either zr (for the horizontal plume state) or zp (for the terrain-
following 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.


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:



                        Q (1 − f p )       ∞      ( z − h + 2mz )2 
                                                                  ieff 
                                                                                ( z + h + 2mz )2  
                                                                          + exp −
                                                           ep                           ep      ieff 
    C p {xr , yr , z} =      ~       Fy ⋅ ∑ exp −                                                         (66)
                         2π u σ zp       m = −∞             σ 2                        σ 2       
                                                 
                                                             2  zp                      2  zp       


where zieff is the height of the upper reflecting surface in a stable layer (see Section 5.3) 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 5.5.1.2.


        The penetrated plume height, hep, is taken as the height of the plume centroid above the
mixed layer and is calculated from eq. (94).



                                                            78
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



                        Q               ∞      ( z − h − 2mz )2            ( z + h + 2mz )2  
                                ⋅ Fy ⋅ ∑ exp −                       + exp −                     ,
                                                        es       ieff                 es       ieff
   Cs {xr , yr , z} =    ~                                                                                 (67)
                      2π u σ zs                          2σ 2                        2σ 2       
                                      m = −∞
                                                            zs                          zs      


where zieff is the effective mechanical mixed layer height, σzs is the total vertical dispersion in the
SBL (see discussion in Section 5.5), and hes is the plume height (i.e., stack height plus the plume
rise - see Section 5.6.2).


        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
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
                                                          79
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:


                                            [
                             zieff = MAX ( hes + .
                                                215σ zs {hes }; zim ) .   ]                 (68)


where σzs in eq. (68) is determined from equations found in Section 5.5.1.2 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.


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

                                                       80
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:


                                              1      − y2 
                        F        =              exp 2                               (69)
                            yC            2πσ y      2σ y  ,


where σy is the lateral dispersion parameter (see Section 5.5). 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:

                                     1
                        𝐹𝐹𝑦𝑦𝑦𝑦 = 2𝜋𝜋𝑥𝑥 ,                                                (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.


       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,




                                                          81
                                        (              )         (
                           Cc , s = CCh 1 − σ r2 σ h2 + CR σ r2 σ h2         ),                   (71)


where σh2 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 σr2/σh2 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 𝑢𝑢, and random components σu
and σv. Thus, a measure of the total horizontal wind “energy” (given that the along wind and
crosswind fluctuations are assumed equal i.e., σu = σv), can be represented as


                           σ h 2 = 2σ~v 2 + u 2 ,                                                 (72)


                          1
where 𝑢𝑢 = (𝑢𝑢�2 − 2𝜎𝜎�𝑣𝑣2 )2 . The random energy component is initially 2𝜎𝜎�𝑣𝑣2 and becomes equal to
σh2 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

                                            2              −𝑥𝑥𝑟𝑟
                          𝜎𝜎𝑟𝑟2 = 2𝜎𝜎�𝑣𝑣2 + 𝑢𝑢 �1 − exp �       �𝑢𝑢�𝑇𝑇 ��,                        (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)
                                                                                                           2
shows that at small travel times, 𝜎𝜎𝑟𝑟2 = 2𝜎𝜎�𝑣𝑣2 , while at large times (or distances) 𝜎𝜎𝑟𝑟2 = 2𝜎𝜎�𝑣𝑣2 + 𝑢𝑢 ,
which is the total horizontal kinetic energy (σh2) of the fluid. Therefore, the relative
contributions of the coherent and random horizontal distribution functions (eq. (71)) are based on
                                                            82
the fraction of random energy contained in the system (i.e., 𝜎𝜎𝑟𝑟2 ⁄𝜎𝜎ℎ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 σr2 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 (σr2/σh2|P and σr2/σh2|NP) into a single effective factor (σr2/σh2|CBL). That is,


                         σ r2                       2
                                           = f p ⋅σ r
                                                                                  2
                                                                 + (1 − f p ) ⋅ σ r                 (74)
                                σ h2 CBL                σ h2 P                        σ h2 NP
                                                                                                ,


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 (σr2/σh2)
replaced by (σr2/σh2|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 5.5.3) 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 4.2), which are denoted by an overscript tilde, e.g., 𝜎𝜎�𝑤𝑤𝑤𝑤 .


                                                                 83
        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 and σ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):


                                        σ 2y , z = σ 2ya , za + σ b2 ,                     (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



                                                          84
direct and indirect sources and σza appears as σaj; for the penetrated source σz appears as σzp in
the dispersion expression.


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:


                                          σ~ x
                         σ        =         v
                             ya                      p
                                             x u~                                      (76)
                                      u~ 1 +      
                                          2T 
                                               Ly 
                                                     ,


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:


                                                          σ~v x
                                            σ ya =                                       (77)
                                                     u~ (1 + αX )
                                                                  p




                                                           85
where 𝑋𝑋(= 𝜎𝜎�𝑣𝑣 𝑥𝑥⁄𝑢𝑢�𝑧𝑧𝑖𝑖 ) 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
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.




                                                      86
 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.


       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 z PG z                                                (78)
                                         max  ,

                                                     87
where zPG = 0.46 m (Prairie Grass release height), and 𝑧𝑧𝑚𝑚𝑚𝑚𝑚𝑚 = max(𝑧𝑧; 𝑧𝑧𝑃𝑃𝑃𝑃 ).o 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.


                                                  h              h 
                                σ zas = 1 − es  σ zgs +  es  σ zes .                   (79)
                                                  zi              zi 


        For hes ≥ zi σzas is set equal to σzes. The expressions for calculating hes are found in
Section 5.6.2. 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):



                                                                   x u~ 
                                                                          12
                                                              
                                σ zes   = σ~  wT   ( x u~ )   1 +       '                (80)
                                                               2TLzs 


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)
                                                              88
                                  l
                          TLzs = ~                                                            (81)
                                  σ wT .


          The length scale l is an interpolation between the limiting length scales for neutral
conditions and stable conditions


                                                   1 1 1
                                                    = + .                                     (82)
                                                   l ln ls


where 𝑙𝑙𝑛𝑛 = 0.36 ℎ𝑒𝑒𝑒𝑒 and 𝑙𝑙𝑠𝑠 = 0.27 𝜎𝜎�𝑤𝑤𝑤𝑤 ⁄𝑁𝑁 . 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 equations (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:


                                                        σ~wT t
                           σ zes =                                               .
                                            σ~
                                                                            12
                                            wT t      1       N                           (83)
                                     1 +                  +         
                                                2  0.36hes 0.27σ~wT  


          Finally, to complete the description of eq. (79), the surface portion of vertical dispersion
(σzgs) in the SBL, is calculated from Venkatram (1992) as



                                  2  u* x  
                                                          −1/ 3
                                                       x
                          σ zgs =              1 + 0.7                                      (84)
                                  π  u~           L
                                                                .


          For the direct and indirect sources in the CBL, the ambient portion of the vertical
                                                          89
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


                                            2       2       2
                                          σ zaj = σ zej + σ zg ,                                   (85)


where the elevated portion (σzej) is obtained from Weil et al. (1997) as


                                                      σ x
                                           σ zej = α b wj~ ,                                       (86)
                                                       u

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


                                                 Hp
                                    
                          σ zg = bc 1 − 10 
                                           
                                                         
                                                          
                                                      zi                   (
                                                             ⋅ (u* u~ ) ⋅ x 2 L
                                                                       2
                                                                                     )             (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
                                                                �)2 𝑥𝑥 2
                                                          (𝑢𝑢∗ ⁄𝑢𝑢
dispersion in the unstable surface layer; i.e., 𝜎𝜎𝑧𝑧 ∝          |𝐿𝐿|
                                                                           . The parameterization is designed

to: 1) agree with Venkatram’s result in the limit of a surface release, 2) provide good agreement

                                                           90
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 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


                                                  0.4∆h
                                           σb =         ,                                    (88)
                                                     2


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 Δℎ𝑝𝑝 = ℎ𝑒𝑒𝑒𝑒 − ℎ𝑠𝑠
where hep is calculated from eq. (94)).




                                                     91
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. EPA, 1995). PRIME partitions plume mass between a cavity
recirculation region and a dispersion enhanced wake region based 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 4),
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
                                                     92
algorithms are used exclusively to calculate concentration with AERMOD-derived ambient
turbulent intensities as input. To ensure 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 = γ CPr ime + (1 − γ ) C AERMOD                   (89)


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:



                            (
                       − x −σ 2
               γ = exp
                                xg    )  − y −σ 2
                                    exp          (
                                                   yg         )
                                                            − z −σ 2
                                                       exp          zg (
                                                                                 )       (90)
                             2
                           2σ xg              2
                                             2σ yg               2
                                                                2σ zg    
                                                                    


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).
                                                         93
                                                                          13
                                         3Fm x   3 Fb x 2 
                                 ∆hd =  2 2 + 2 ⋅ 3                                    (91)
                                         β1 u p 2β1 u p 


where Fm = (T/Ts)ws2rs2 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}.


        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


                                                              12
                                             2 Fb zi             x
                                     ∆hi =                        ,                    (92)
                                              α u r  r
                                             r p y z             up


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 product of crosswind dimensions of the
assumed elliptical plume is calculated from Weil et al. (1997) as


                                                     a    λ3y 2 w*2 x 2
                                        ry rz = rh2 + e            ⋅                       (93)
                                                          4            u2p
                                                                               ,


where 𝑟𝑟ℎ = 𝛽𝛽2 (𝑧𝑧𝑖𝑖 − ℎ𝑠𝑠 ), β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).



                                                      94
       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


                                             hs + zi
                                     hep =           + 0.75∆heq .                               (94)
                                                2

where Δheq is defined in eq. (56).


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 been
modified. With this modification, AERMOD computes stable plume rise, Δhs, from Weil et al.
(1988b) as


                                     1 3                                              1 3
                           F               N 'F        N 'x              N 'x 
              ∆hs = 2.66  2b           ⋅     m
                                                    sin       + 1 − cos            ,     (95)
                           N up              F
                                             b          u p               u p  


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
(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


                                               up         − F N ′
                                        xf =      arctan  m  .                          (96)
                                               N′         Fb 


        Upon substituting xf eq. (96) for x in eq. (95) the maximum final rise of the stable plume
Δhs{xf} reduces to:


                                                                     13
                                                           F 
                                       ∆hs {x f } = 2.66  b 2  .                      (97)
                                                           up N 


        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, final stable plume rise (eq. (97))
is limited by the neutral rise calculated from Weil (1985) as


                                12 L3n 5 ( hs + .
                                               12 Ln )
                                                            25
                           ∆hs = .                                                        (98)
                                                                 ,


where the neutral length scale 𝐿𝐿𝑛𝑛 = 𝐹𝐹𝑏𝑏 ⁄�𝑢𝑢𝑝𝑝 𝑢𝑢∗2 �.




                                                            96
          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,


                                                          1
                                                     4 Fb 4
                                             ∆hs =        3
                                                                  .                        (99)
                                                      N       4




          Finally, the stable plume rise is limited by a calculation of the unstable rise (see Section
5.6.1).


5.7 Source characterization


          AERMOD gives the user the ability to characterize a source as either a point, a volume,
an area, a roadway, or a buoyant line. Additionally, AERMOD has the capability of
characterizing irregularly shaped area sources.


          Point sources are characterized exactly as in the ISC3 model (U.S. EPA, 1995). The
input to the model includes the location, elevation, emission rate, 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:



                                                       97
                          σ 2y = σ 2yl + σ 2yo                                           (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.


        As was described in section 2.0, Model Overview, the buoyant line source algorithm
from the Buoyant Line and Point Source (BLP) model (Schulman and Scire, 1980) was
incorporated into AERMOD beginning with version 15181. The input to the model includes
emission rate and average release height for individual lines, and average building length,
average building height, average building width, average line source width, average building
separation, and average buoyancy for the single source. Users should refer to the BLP User’s
Guide for a description of the formulation of the BUOYLINE source type in AERMOD.


        The RLINE source type was added beginning with version 19191 to better characterize
mobile source emissions (i.e., roadways). The RLINE source type is an adaptation of the
Research Line (R-LINE) model (Snyder and Heist, 2013). The formulation of the RLINE source
type was updated in version 23132 to harmonize the RLINE source type with other source types
(POINT, AERA, and VOLUME) in AERMOD which also resulted in an update to the
dispersion coefficients. For a complete description of the formulation of the RLINE source type
in AERMOD, refer to Snyder and Heist (2013) and EPA (2023e).



                                                       98
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.
EPA (2015).


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 using 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

                                                     99
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:

                                   3
                        𝜎𝜎𝑟𝑟 = (𝜎𝜎𝑟𝑟𝑟𝑟      3 )1⁄3
                                       + 𝜎𝜎𝑟𝑟𝑟𝑟    .                                     (101)


       The buoyancy-induced dispersion term, σrb, is calculated in AERMOD as

                        𝜎𝜎𝑟𝑟𝑟𝑟 =
                                   0.4 𝛥𝛥ℎ
                                              ,                                          (102)
                                        √2



where Δh is the plume rise. Relative dispersion due to ambient turbulence, σra, is parameterized
by

                                            1 3
                                   𝑎𝑎1 𝜀𝜀 2 𝑡𝑡 2                                         (103)
                        𝜎𝜎𝑟𝑟𝑟𝑟 =      𝑎𝑎 𝑡𝑡        ,
                                    1+ 2
                                      𝑇𝑇𝐿𝐿𝐿𝐿




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

                                             𝑧𝑧
                        𝑇𝑇𝐿𝐿𝐿𝐿 = 𝑎𝑎𝑟𝑟1 𝜎𝜎 𝑖𝑖 ,                                           (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):


                        𝜀𝜀 =
                                   2
                               𝑏𝑏𝜎𝜎𝑤𝑤
                                        ,                                                (105)
                               𝑇𝑇𝐿𝐿𝐿𝐿




                                                       100
where b is a constant (= 0.78). The values of wind speed (U) and σw used in eq. (103) through
eq. (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.


       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:



                                                       1
                                   3
                        𝜎𝜎𝑟𝑟 = (𝜎𝜎𝑟𝑟𝑟𝑟      3
                                       + 𝜎𝜎𝑟𝑟𝑟𝑟 + 𝜎𝜎0 )3 ,                               (106)




where σ0 is the initial dispersion coefficient of the volume source calculated as �𝜎𝜎𝑦𝑦0 𝜎𝜎𝑧𝑧0 based

                                                             101
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 along-wind 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.


5.8.1.4 Defining extent of plume


       The original implementation of PVMRM in AERMOD used four (4) times the relative
dispersion coefficients to define the plume volume, which accounts for more than 99.99% of the
plume. The original implementation also used a minimum value 5.0 m for the plume radius to
account for near-source effects (e.g., downwash), which is consistent with the initial plume
volume assumed by Hanrahan (1999). 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 16216, 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, and a minimum plume radius of 15.0 m. 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, and minimum plume radius of 9.0 m. Since AERMOD incorporates a
horizontal meander algorithm that increases lateral plume spread beyond that accounted for


                                                    102
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.,


                                        𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆
       𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = min �2.15,1,282 ∗ � 𝑆𝑆𝑆𝑆 ��,


where NSUBZ is the number of sigmas used to define the plume volume and SYEFF is the
effective σy value that replicates the plume centerline associated with meander but based on a
standard Gaussian plume calculation. An upper bound of 2.15 is applied, which corresponds with
the point at which the concentration is 10 percent of the centerline value and is used to define the
extent of the plume in other contexts (e.g., initial dispersion coefficients for volume sources).
The minimum value of the dispersion coefficient is also adjusted based on the ratio
1.282/NSUBZ, where NSUBZ is the adjusted value.


       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 16216, 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 5.1.

                                                     103
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.


5.8.2 Minimum ozone concentration for stable conditions


       Hanrahan (1999a) applied a minimum ozone concentration of 40 ppb for stable
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, with an upper limit of 40 ppb. If the
Monin-Obukhov length is positive and the value is over 250 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
                                                     104
interpolation, i.e., the minimum value is calculated as MIN(40ppb, O3MAX24) * (250 – L)/200,
where O3MAX24 is the maximum ozone concentration over the previous 24 hours, and L is the
Monin-Obukhov length in meters.


       This approach is primarily applicable for tall stacks, where ozone concentrations aloft
may be higher than surface measurements, as surface ozone will be titrated overnight by surface
emissions of NO but can maintain daytime-high concentrations above the mixed layer (U.S.
EPA, 2006). The option can be turned off via the NOMINO3 keyword, which will use the
hourly measured ozone for NO conversion rather than the value computed as described above.
The use of the NOMINO3 option may generally be applicable for surface releases and shorter
stacks, where the plume heights are expected to be within the surface layer during typical
nighttime conditions. As with all NO2 conversion options, this option should be used in
consultation with the appropriate reviewing authority.


5.9 Generic Reaction Set Method (GRSM)


       The Generic Reaction Set Method (GRSM) was introduced into AERMOD as an
ALPHA option in version 21112 and upgraded to a BETA option in version 22112 with
additional updates in version 23132. Carruthers et al., 2017, document the development of a
technique that accounts for the equilibrium between NO, NO2, and ozone in the atmosphere
called the Atmospheric Dispersion Model Method, ADMSM. The ADMSM uses similar
calculations for plume entrainment as PVMRM but adds a “reaction rate” based on solar
radiation and travel time from source to receptor. The reaction rate is based on the generic
reaction set (GRS) chemistry scheme, which is a semi-empirical photochemical model
developed originally by CSIRO in Australia (Azzi and Johnson, 1992; Venkatram et al., 1994)
for multiple step conversions between NO, NO2, and O3 which limit NO conversion based on
the available ozone.




                                                   105
         Reaction rates are determined by three equations (see equations 4-8 in Carruthers et al.,
2017):


                d[NO2]/dt = k[NO][O3] – J[NO2] (Equation B-1)

                d[NO]/dt = J[NO2] – k[NO][O3]          (Equation B-2)

                d[O3]/dt = J[NO2] – k[NO][O3]          (Equation B-3)


         Where the reaction rate coefficient k is defined as:

                k = 4.405*exp(-1370/Tamb)              (Equation B-4)


         And the photo-dissociation rate J are defined as:


                J = 8*10-4*exp(-10/Q) + 7.4*10-6*Q              (Equation B-5)


         Where Q is solar radiation in W/m2.


         A fifth order Runge-Kutta scheme with adaptive time-stepping is used to solve the
equations, to give ‘post-chemistry’ concentrations at each receptor. The post-chemistry NO and
NO2 values at each receptor are partitioned between the sources, based on the partitioning of the
pre-chemistry NOX, to give post-chemistry NO and NO2 instantaneous plume concentrations for
each receptor and source combination. The post-chemistry ensemble plume concentrations are
then calculated from these post-chemistry instantaneous plume concentration values, again
using the ratio of the plume cross sectional areas.


         The GRSM method was further revised in AERMOD version 23132 with minor updates
and re-evaluated as described in Stocker et al. (2023) and U.S. EPA (2024b).




                                                      106
5.10 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 are largest at night and relatively absent during the day (Oke, 1998). An urban
“convective-like” 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
Oke (1973; 1982) is used to construct an empirical model. Oke presents observed urban-rural
temperature differences for several 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 humidity. An empirical fit
to the data yields the following relationship:



                              ∆Tu − r = ∆Tmax 01              10,
                                                 . ln  P  + .                     (107)
                                                      Po      


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.


                                                    107
       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:


                                       H u = α ρ c p ∆Tu − r u* ,                        (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 oC . 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 area 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 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:


                                            ziu ≈ R1 2 and R ≈ P1 2 ,                    (109)



                                                     108
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 (109) leads to the following equation for the nocturnal urban boundary layer
height due to convective effects alone:


                        ziuc = ziuo ( P Po )
                                             1/ 4
                                                                                       (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:


                                        ziu = Max [ ziuc ; zim ].                      (111)

       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,


                                                 1/ 3
                                g H u ziuc 
                        w*u =                                                      (112)
                                ρ cp T                ,


where w*u is the urban nighttime convective velocity scale and T is the near-surface air

                                                            109
temperature.


       Having estimated w*u the turbulence in the nighttime urban can be enhanced using the
expressions found in Section 4.1.5. However, since for low level sources σwT depends primarily
on u* (see eq. (34) and eq. (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 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 4.1.5. That is, σwc and σwm are calculated from eq. (35) and eq. (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, an urban nighttime Monin-Obukhov length is
calculated using eq. (8) with substitutions u*u for u* and Hu (eq. (108)) for H. This results in a
convective urban Monin-Obukhov length. For dispersion calculations, during stable hours, the
absolute value of the urban Monin-Obukhov length is used, while daytime transition hours will
use the original convective urban Monin-Obukhov length. To maintain consistency, during
stable hours, the most neutral Monin-Obukhov length is selected between the urban nighttime
Monin-Obukhov length and the rural Monin-Obukhov length, resulting in a stable-to-neutral
atmosphere. Similarly, during daytime transition hours, the most convective Monin-Obukhov
length value is selected between the urban Monin-Obukhov length and the rural Monin-
Obukhov length, resulting in a neutral-to-convective atmosphere.


       Finally, the total nighttime turbulence in the urban boundary layer is calculated as the
sum (in quadrature) of the convective (urban) and mechanical (rural) portions. With these
                                                     110
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 = 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 4.1.3) for all heights above ziu, and is assumed to be equal to a
small positive value below ziu; i.e.,


                         ∂θ ∂z = 10 − 5                         for z ≤ ziu
                                                                                           (113)
                         ∂θ ∂z = rural value                     for z > ziu .


        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 eq. (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. When AERMOD was originally
promulgated, plume height in these conditions was not allowed to exceed 1.25 ziu. This created
potential issues for tall stacks with buoyant plumes in urban areas, as the limit of 1.25 ziu,
artificially limited the plume rise of these types of sources. With 15181, a change was made to
the urban formulation to no longer limit the plume rise but treat such sources as pseudo
penetrated plume sources to allow for the plumes to rise above the urban boundary layer when
the stack height plus initial plume rise was greater than or equal to the urban mixing height. In
version 22112, a further modification was made to calculate plume penetration parameters when
the stack height plus plume rise is greater than or equal to the urban mixing height and the stack
height is less than the urban mixing height. This change was made to avoid divide by zero in

                                                      111
the plume penetration calculations. For more information about AERMOD’s penetrated plume
calculations, see Section 5.2.3.


       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). However,
beginning with AERMOD version 11059, the urban boundary layer options are still
implemented for early daytime hours when the rural convective mixing height is less than the
urban mixing height calculated by Equations 110 and 111. This change to continue the urban
boundary layer options for early daytime hours was made because for early daytime hours, the
lack of urban adjustments could result in an unrealistic drop in mixing heights for urban sources
when the urban boundary layer height is still larger than the rural convective mixing height.


       An analysis of the Indianapolis, IN urban SF6 study used as part of the AERMOD
evaluation shows that the continued application of the urban adjustments for early morning or
early evening hours increased model performance compared to observations. Observed Robust
Highest Concentrations (RHC) for Indianapolis were 5.5 µg/m3, while RHC values without the
transition procedure were 4.3 µg/m3 and RHC values with the transition procedure were 4.5
µg/m3. Limiting the analysis to convective hours, which are the hours most affected by the
transition procedure, showed observed RHC values of 5.5 µg/m3, 4.2 µg/m3 with no transition
procedure, and 4.7 µg/m3 with the transition procedure.


       Urban boundary layer adjustments were implemented in version 19191 for the RLINE
and RLINEXT source types which require the BETA and ALPHA flag, respectively. Based on
algorithms for the AERMOD area source, modification was introduced so when a RLINE or
RLINEXT source is specified as an urban source, the u* urban nighttime Monin-Obukhov
length, ziu, and zo parameters are recalculated for morning transition and stable hours. As initial
integration of the RLINE and RLINEXT source types in version 19191 is based on algorithms
of version 1.2 of the Research LINE source (R-LINE) model (Snyder et al, 2013), the R-LINE
methodology for calculation of w*, is included. The methodology incorporates a check to set
                                                  112
any negative w* values to 0.0, associated with stable meteorological conditions, unlike other
AERMOD source types that accept negative w* values. The resultant w* values are then used in
the calculation of σv for RLINE and RLINEXT sources only.


       Urban boundary layer adjustments were also implemented with the BUOYLINE source
type as an ALPHA option in version 19191. The ALPHA requirement was removed for the
BUOYLINE source type in version 22112 and the PG Stability class assignment was corrected
for the BOUYLINE source when the urban option is used. PG stability class is set to 4 for stable
hours when the urban boundary layer adjustments are applied.


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 AERMOD-
                      derived 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)
                                                   113
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)
Δ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)
                                    114
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)
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)
𝑤𝑤𝑗𝑗    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)
                                    115
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)
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 equations (53),(59) and (65)
                                         116
ρ         density of air (kg m-3)
σb        buoyancy induced dispersion for the direct & indirect sources (m)
σh2       total horizontal wind “energy” used in the meander algorithm (m2)
σr2       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)
σ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)
                                       117
φ     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)




                                  118
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:


      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):


                                                 0.5
                                1 − θ* 0.09 
                         neq =                                                      (114)
                                     0.5             ,


                                                           119
where θ* is the temperature scale as calculated from eq. (18).


      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.
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. Values for albedo and Bowen ratio should be based on a 10 km x 10
km area centered on the meteorological tower. Values for surface roughness should be based on
a 1 km radius around the meteorological tower and can be differentiated for up 12 wind sectors.
Surface characteristics should be determined with the aid of the AERSURFACE program (U.S.
EPA, 2026b). Refer to the AERSURFACE and AERMET User’s Guides (U.S. EPA, 2026a) for
additional information on determining appropriate surface characteristic values.


      Other

           •     Latitude;
           •     Longitude;
           •     Time zone;
           •     Wind speed instrument threshold for each data set ( uth ).

      Optional

           •     Solar radiation;
           •     Net radiation (Rn);
           •     Profile of vertical turbulence (σw);

                                                        120
           •   Profile of lateral turbulence (σv)

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

      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.


      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.


      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, therefore, not used in the construction of the wind speed profile (profiling of winds is
accomplished in AERMOD).



                                                    121
       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).


       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 A.9, these
limits are applied to the effective values of turbulence and wind speed.


       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.


                                                    122
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 to estimate zic.
This is accomplished as follows:
                                                       123
    •   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
                o the height of the actual sounding top,
                o that dθ/dz has been extrapolated above the sounding zic, and
                o 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 applied such that:



                                       [                 ]
                           σ w {z} = Max σ w {z}; 0.02 m s −1
                                                                                           (115)
                           σ {z} = Max [σ {z}; 0.05 u{h }; 0.2 m s ].
                            v               v               s
                                                                         −1



                                                     124
These limits are also applied when selecting the turbulence for plume rise calculations.


       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


                        uv = ( u + u′, v ′ ) ,                                         (116)


where 𝑢𝑢 is the mean measured wind, and the primed quantities refer to the turbulent fluctuations.
The assumption being made is 𝑢𝑢𝑣𝑣 = 𝑢𝑢. 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


                                        uv2 = u 2 + σ v2 + σ u2 .                      (117)


       In this simple model, uv , is the dilution wind. If we take σu = σv, the dilution wind can
be written as


                                     u{hs } = u{hs }2 + 2 σ v2 .                       (118)


       This formulation assures that the dilution wind is not zero as long as either 𝑢𝑢 or σv is not
zero. Similarly, at the time of plume rise calculations, the effective turbulence and effective
wind speed will be recalculated using equations (115) and (118), where the turbulence and
winds will be evaluated at stack top.




                                                      125
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.




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:


                                                   126
            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;

            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. 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
                                                     127
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.




                                                     128
Appendix B. Description of ALPHA Options in AERMOD 26135


       This appendix to the AERMOD Model Formulation document provides descriptions of
the ALPHA options included in version 26135 of the AERMOD modeling system. Because
these model options are not part of the promulgated regulatory formulation of AERMOD (U.S.
EPA, 2024a), they are not included in the main text of the AERMOD Model Formulation
document. It should be noted that there are no BETA options implemented within version
26135 of the AERMOD Modeling System. Refer to the AERMOD User’s Guide (U.S. EPA,
2026c) for details on the required control file syntax to experiment with the ALPHA options
described in this appendix.


B.1 Characterization of ALPHA and BETA Options

       BETA options are more mature with respect to where they are along the path of
development than ALPHA options. The scientific basis of BETA options is well established,
implementation in the model is considered complete, and performance evaluations demonstrate
their validity and justification for inclusion as a regulatory update to AERMOD. BETA options
have generally met the criteria for justification as an alternative model (with appropriate
approval from the EPA Regional Office with concurrence from the EPA Model Clearinghouse)
and proposal as a regulatory formulation update to AERMOD.


       In contrast to BETA options, ALPHA options are less mature in the development
process than BETA options and can range anywhere from an idea or concept that EPA has
identified to be explored as a priority development topic but is still in the initial research and
experimental phases, to an option that is far along the development path but has not yet met all
of the criteria necessary for justification as an alternative model and proposal as a regulatory
update to the model formulation. In some cases, a new model option may be implemented as an
ALPHA option that is not considered a change to the scientific formulation but needs to be more
fully tested by the broader user community before it is implemented for use. ALPHA options
that represent changes in the scientific formulation of the model may eventually reach maturity
                                                    129
in development and to be proposed and promulgated as scientific changes to the formulation of
AERMOD while others may never develop that level of scientific maturity and eventually
removed from the model. Each of the ALPHA options described below in version 26135
represent potential updates to the scientific formulation of AERMOD and must go through the
regulatory update process which includes rulemaking to update the Guideline before they are
considered part of the scientific formulation of the regulatory version of AERMOD.


B.2 ALPHA Options in AERMOD 26135

       In contrast to BETA options, ALPHA options are less mature in the development
process than BETA options and can range anywhere from an idea or concept that EPA has
identified to be explored as a priority development topic but is still in the initial research and
experimental phases, to an option that is far along the development path but has not yet met all
of the criteria necessary for justification as an alternative model and proposal as a regulatory
update to the model formulation. In some cases, a new model option may be implemented as an
ALPHA option that is not considered a change to the scientific formulation but needs to be more
fully tested by the broader user community before it is implemented for use. ALPHA options
that represent changes in the scientific formulation of the model may eventually reach maturity
in development and to be proposed and promulgated as scientific changes to the formulation of
AERMOD while others may never develop that level of scientific maturity and eventually
removed from the model. Each of the ALPHA options described below in version 23132
represent potential updates to the scientific formulation of AERMOD and must go through the
regulatory update process which includes rulemaking to update the Guideline before they are
considered part of the scientific formulation of the regulatory version of AERMOD.


      RLINE Extended (RLINEXT)

       The RLINEXT source type refers to an extended version of the RLINE source type and
was implemented in version 19191, concurrently with the RLINE source type, but as an
ALPHA option. As an extension of the RLINE source type, RLINEXT accepts additional input

                                                     130
parameters to represent depressed roadways (RDEPRESS), solid roadside barriers
(RBARRIER) where each segment is a separate barrier, continuous solid barriers (SBARRIER)
where multiple segments define a single solid barrier, barrier edge effect (BAREDGE), and
perform wind profile calculations without a displacement height (RLINEFDH). Each of these
are discussed separately below.


      Solid Roadway Barrier (RBARRIER, SBARRIER)

       The RBARRIER option simulates the effects of solid roadside barriers on either or both
sides of the roadway and can only be used concurrently with an RLINEXT source type. The
development of the barrier algorithms is ongoing and are based on the work of Schulte et al,
2014; Ahangar, 2017; and Venkatram et al., 2021. Algorithm development has made use of
measurements taken during several different experiments including studies conducted in EPA’s
wind tunnel (Heist et al., 2009) and field studies including the Idaho Falls experiment (Finn et
al., 2010), the Raleigh near road study in 2006 (Baldauf et al., 2008), and the Phoenix, Arizona
field study (Baldauf et al., 2016). Barrier input parameters are specified with the RBARRIER
keyword on the SO pathway.


       As an alternative to using the RBARRIER keyword, beginning with AERMOD version
26135 solid roadside barriers can also be input using barrier end points with the SBARRIER
keyword. Barriers can be broken into segments (optional) with varying lengths, heights, and
elevations. SBARRIER provides AERMOD with the needed inputs to calculate the barrier
Htwall and DCL as used in RBARRIER. The input barrier segment groups (Barid) will be
treated as a continuous barrier.


      Barrier Edge Effects

       Beginning with AERMOD version 26135, the option to treat the transition at barrier
endpoints can be enabled with the BAREDGE keyword as a model option (MODELOPT) on
the CO pathway. Enabling this option assigns an edge barrier influence zone based on the

                                                   131
location of the receptor relative to the barrier. BAREDGE increases concentrations at the edges
of the solid barriers where the plume wraps around the end of the barrier (Venkatram et al.,
2021).


      Vegetative Barriers

         Beginning with AERMOD version 26135, a vegetative barrier, based on Hashad et al.,
2024, can be input using the VBARRIER keyword. Dispersion from vegetative barriers is a
function of vegetation parameters: height, width, leaf area index, leaf area density. The
VBARRIER only influences dispersion when it is located downwind of the source for a given
wind direction. The VBARRIER keyword cannot be used in combination with RBARRIER or
SBARRIER for the same source ID (Srcid).


      Depressed Roadways (RDEPRESS)

         The RDEPRESS option accounts for dispersion associated with depressed roadways and
can only be used concurrently with an RLINEXT source type. The development and refinement
of the RDEPRESS option is ongoing by EPA’s ORD, utilizing wind tunnel studies described in
Heist et al, 2009. Road depression input parameters are specified with the RDEPRESS keyword
on the SO pathway.


      RLINE/RLINEXT Displacement Height (RLINEFDH)

         The RLINEFDH option modifies the vertical wind profile calculations for the RLINE
and RLINEXT source types without using a displacement height. This was added as an
ALPHA option for the RLINE and RLINEXT source types to make the wind profile
calculations more consistent with other AERMOD source types which do not consider a
displacement height in the wind profile calculations. This option is enabled by including the
secondary keyword RLINEFDH on the CO pathway as a model option specified with the
MODELOPT primary keyword.


                                                   132
      Building Downwash Options

        Analyses have shown AERMOD to both overpredict and underpredict ground-level
concentrations in the building wake, depending on the building dimensions; stack height; stack
location; and the orientation of the building relative to the wind direction. Overprediction and
underprediction have been demonstrated in analyses of single, one-tiered rectangular buildings.
Some examples in which AERMOD has been shown to be deficient with regard to building
downwash include elongated buildings, buildings that are angled rather than perpendicular to
the wind, and buildings with stacks located near a building corner (Perry et al., 2016; Petersen et
al., 2017).


        The building downwash algorithms in AERMOD are based on solid, square and
rectangular ground-based buildings. Porous, streamlined, and lattice-type structures that are
common at many sites have been shown to have a different influence on flow and dispersion
than solid buildings. Currently, these types of structures can only be modeled in AERMOD as
solid buildings which are not representative.


        EPA’s Air Quality Modeling Group (AQMG) has collaborated on multiple research
initiatives to improve AERMOD’s performance with respect to point emission sources that are
subject to building downwash. These include initiatives lead by EPA’s Office of Research and
Development (ORD), the PRIME2 Advisory Subcommittee (PRIME2) within the Atmospheric
Modeling and Meteorology (APM) Committee of the Air and Waste Management Association
(A&WMA), the Sibley School of Mechanical and Aerospace Engineering at Cornell University,
and the Bureau of Ocean Energy Management (BOEM).


        ORD has performed wind tunnel experiments and embedded large eddy simulations
(LES) to better understand how to parameterize buildings that are elongated and angled relative
to the wind flow and the parameterization of the plume in the cavity and far wake regions. The
ORD studies are concentrated on single rectangular buildings, specifically investigating changes
in plume parameters at discrete downwind distances from the building and source, longitudinal
                                                133
and lateral plume profiles, the lateral plume shift on the lee side of rotated buildings which
resulted in recommended modifications to the current PRIME algorithm in AERMOD.


       A&WMA involved the collaboration of technical experts, industry groups, and
representatives from the regulatory agencies with the purpose of (1) providing a technical
review forum to improve the PRIME building downwash algorithms in AERMOD; and (2)
establishing a mechanism to review, approve, and implement new science into the model.
A&WMA’s research has included the reanalysis of existing wind tunnel data, as well as the
completion of new wind tunnel experiments to investigate the decay of the building wake above
the top of the building, appropriate height at which approach turbulence and wind speed are
calculated, the reduction of wake effects for streamlined structures, and the effect of approach
roughness on the wake. Their analyses have led to recommendations for new turbulence
enhancement and velocity deficit equations that address these aspects, as well as an update to
the degree of entrainment into the plume and the application of plume rise at the receptor when
the plume is affected by downwash. Their analyses have led to recommendations for new
turbulence enhancement and velocity deficit equations that address these aspects, as well as an
update to the degree of entrainment into the plume and the application of plume rise at the
receptor when the plume is affected by downwash.


       Researchers at Cornell University have focused on the lateral plume shift on the lee side
of rotated buildings which have prompted the addition of a new source type in AERMOD to
further the research to better parameterize this phenomenon caused by building downwash.


       BOEM sponsored the integration of the platform downwash algorithm in the Offshore
Coastal Dispersion (OCD) model into AERMOD to further the research needed to enable
AERMOD to model to simulate dispersion from downwash effects caused by raised lattice and
porous type structures such as offshore platforms.




                                                     134
       The ALPHA options related to building downwash effects currently in AERMOD that
have resulted from these different research initiatives are described in the sections that follow.


B.2.7.1 ORD_DWNW


       Three ALPHA options have been implemented into AERMOD as a result of the research
conducted by ORD and are enabled in the model control file with the use of the primary
keyword ORD_DWNW on the CO pathway. These options attempt to address the following:
(1) resolve a mismatch in plume width at the transition between the cavity and the far wake; (2)
use an effective wind speed for the primary plume (instead of stack height wind speed); and (3)
adjust the maximum ambient turbulence levels (Monbureau et al., 2018).


ORDCAV - Plume Spread Matching at Boundary Between Near-wake and Far-wake


       The current implementation of PRIME in AERMOD creates a cavity plume and a re-
emitted plume to simulate two distinct regions with a weighted distribution of mass between the
two plumes. The cavity and re-emitted plumes initially have the same vertical dispersion on the
leeward side of the building. The re-emitted plume is then allowed to grow with downwind
distance while the vertical dispersion of the cavity plume remains unchanged throughout the
cavity. This creates a discontinuity of the two plumes at the near-wake boundary that results in a
reduction in ground level concentrations near cavity boundary. ORD has proposed eliminating
this discontinuity such that the vertical dispersion of both plumes is equal at the cavity
boundary. This change would result in an increase in the re-emitted plume concentration. Figure
B-1 illustrates the mismatch at the cavity boundary (LR) in the vertical dispersion (𝜎𝜎z) of the re-
emitted plume (solid line) which grows from the lee wall through the cavity and the cavity
plume (dotted line) which is constant through the cavity. This option was implemented as an
ALPHA option in AERMOD beginning with version 19191 and is specified in the AERMOD
control file with the combination of the ORD_DWNW and ORDCAV keywords.


                                                    135
           Figure B-1.  Discontinuity in vertical dispersion (σz) between re-emitted plume (solid line)
           and the cavity plume (dotted line) at the cavity boundary (LR). Inset: the vertical
           concentration profile within the cavity (left panel); the profile after plume is emitted from
           the cavity (solid line, original; dotted line, revised) (Monbureau et al., 2018).

ORDUEFF - Effective Wind Speed


       The PRIME model currently uses the wind speed at stack height for the primary plume
which likely under-estimates ground-level concentrations in the lee of the building. This is
inconsistent with other source types in AERMOD where an effective wind speed (Ueff) is used.
Ueff is defined as the averaged profiled wind speed between the receptor height and the plume
centerline allowing the wind speed of the plume to change with a changing environment, such
as the variability in surface roughness and stability. ORD has proposed the use of Ueff in PRIME
which would result in a decrease in the wind speed toward the ground, i.e., the wind speed
within the near wake, thus increasing cavity concentrations (Monbureau et al., 2018). This
option was implemented as an ALPHA option in AERMOD beginning with version 19191 and
is specified in the AERMOD control file with the combination of the ORD_DWNW and
ORDUEFF keywords.


ORDUTURB - Maximum Turbulence


       Vertical and lateral dispersion coefficients are based on the formulations of Weil (1996).
When PRIME was implemented in AERMOD, the maximum value of the ambient turbulence
intensity in the wake was reduced from Weil's published value of 0.07 to 0.06. ORD has
                                                 136
proposed increasing this maximum turbulence value to Weil’s original value of 0.07. Increased
turbulence brings the plume down toward the ground more quickly and increases the dispersion.
This will result in reduced concentrations for the primary source and a shorter downwind
distance to the maximum ground-level concentrations (Monbureau et al., 2018). This option was
implemented as an ALPHA option in AERMOD beginning with version 19191 and is specified
in the AERMOD control file with the combination of the ORD_DWNW and ORDUTURB
keywords.


B.2.7.2 AWMADWNW


       As a result of the research performed by A&WMA, five ALPHA options have been
implemented, two of which are variations of a single update. The updates are related to the
effective wind speed of the primary plume, similar to ORD’s proposal, new formulations for
turbulence enhancement and velocity deficit which incorporate building shape (streamline
versus rectangular) and approach turbulence intensity (surface roughness), which will
subsequently affect the calculations of plume rise and the dispersion coefficients in PRIME
(Petersen and Guerra, 2018), and the amount of entrainment into the plume affected by
downwash.


AWMAUEFF - Effective Wind Speed


       Similar to ORD, A&WMA also proposed to redefine how the wind speed is determined
for the primary plume, which is currently defined at the stack height. A&WMA has proposed
using the wind speed at the height of the plume centerline (Petersen and Guerra, 2018). This
option was implemented as an ALPHA option in AERMOD beginning with version 19191 and
is specified in the AERMOD control file with the combination of the AWMADWNW and
AWMAUEFF keywords.


AWMAUTURB/ AWMAUTURBHX - Turbulence Enhancement and Velocity Deficit

                                                  137
       In their analyses of wind tunnel data, A&WMA has demonstrated that building wake
effects decay rapidly back to ambient levels above the top of the building versus the current
implementation in AERMOD in which wake effects can extend up to three building heights. In
addition, the PRIME algorithm currently sets the lateral turbulence in the cavity equal to vertical
turbulence, whereas the wind tunnel shows that lateral turbulence is less than the vertical
turbulence. Finally, the current PRIME algorithm does not account for the effects of approach
roughness on the wake. A&WMA has demonstrated that wake effects decrease as the approach
roughness increases (Petersen et al., 2017).


       Further, AERMOD is not able to properly characterize streamlined structures such as
storage tanks and cooling towers. A wind tunnel study conducted by Petersen (2014) shows that
dispersion is reduced in the wake of streamlined structures and that AERMOD may
underpredict ground-level concentrations when these structures are modeled as rectangular
buildings.


       A&WMA has developed and proposed new formulations for turbulence enhancement
and velocity deficit that address several of the limitations of PRIME previously mentioned,
including rate of decay for building wake effects above the building, lateral turbulence
enhancement in the wake, approach turbulence and wind speed, and wake effects for
streamlined structures (Petersen and Guerra, 2018). Beginning with version 19191, two
additional ALPHA options proposed by A&WMA were incorporated into AERMOD. The first
enables A&WMA’s new formulation and is specified with the keyword AWMAUTURB, along
with the keyword AWMADWNW, as mentioned above. The AWMAUTURB option uses the
minimum of the final momentum plume rise or a representative PRIME plume rise height for all
calculations. It also uses the final momentum plume rise height used to compute effective wind
speed (UEFF), effective 𝜎𝜎w (SWEFF), effective 𝜎𝜎v (SVEFF), effective potential temperature
gradient (TGEFF), and initial turbulence intensities (ambiy and ambiz) and computes mean
wind speed, 𝜎𝜎w, and 𝜎𝜎v at 30 meters (U30, SW30 and SV30, respectively). The second option,


                                                   138
STREAMLINE, when used in combination with the AWMAUTURB keyword, modifies the
A&WMA equations for streamlined buildings.


       Beginning with version 21112, the A&WMA’s formulation was extended to use final
momentum plume rise height to initially compute effective wind speed (UEFF), effective 𝜎𝜎w
(SWEFF), effective 𝜎𝜎v (SVEFF), effective potential temperature gradient (TGEFF), and initial
turbulence intensities (ambiy and ambiz) and then uses the PRIME computed plume rise at each
downwind distance, and as with AWMAUTURB, computes mean wind speed, 𝜎𝜎w, and 𝜎𝜎v at 30
meters (U30, SW30 and SV30, respectively). This newer ALPHA option is enabled by using the
keyword AWMAUTURBHX. The STREAMLINE option is available for use with
AWMAUTURB and AWMAUTURBHX.


AWMAENTRAIN - Entrainment Constants


       The AERMOD PRIME algorithm defines two entrainment constants, alpha (A) and beta
(B), to solve equations for conservation of mass, energy, and momentum. These are not to be
confused with the ALPHA and BETA option designations in AERMOD. The programmed
constant values for alpha and beta constants used in the PRIME algorithm in AERMOD are 0.11
and 0.6, respectively. As discussed previously, AERMOD has been shown to overpredict in
some cases. It has been proposed by A&WMA that this overprediction is due in part to an
underprediction of plume rise, and the PRIME plume rise algorithm has not been tested against
field or wind tunnel observations for building wake applications. A&WMA also suggests that
the entrainment constants in the PRIME algorithm differ from those determined by other
researchers. Based on wind tunnel observations and variation in the values of the alpha and beta
entrainment constants (results to be published), A&WMA found that best performance was
based on a change in the beta entrainment constant value from 0.6 to 0.35. Beginning with
version 21112, a new ALPHA option was included in AERMOD based on A&WMA’s
proposed change to the beta entrainment constant of a value of 0.35. The AERMOD


                                                  139
programmed value of 0.6 can be overridden with the value 0.35 by specifying the keyword
AWMAENTRAIN in combination with AWMADWNW.


B.2.7.3 PLATFORM


       The platform downwash algorithm from OCD has been implemented in AERMOD as an
ALPHA option to account for the treatment of building downwash effects from offshore
platforms with the current AERMOD single boundary layer parameterization limitation. The
platform downwash parameterization can be applied to the POINT, POINTHOR, and
POINTCAP source types in AERMOD using the primary keyword PLATFORM in the SO
pathway.


       The OCD model was originally developed by the Mineral Management Service (MMS),
now BOEM, in the 1980s for estimating plume transport and dispersion in marine and coastal
environments (Hanna et al., 1985). OCD, like AERMOD, is a steady-state Gaussian dispersion
model that simulates plume transport and dispersion based on boundary layer meteorological
parameters. Like many of AERMOD’s predecessors, boundary layer parameterization is based
on the older P-G stability classes.


       The original platform downwash algorithms were developed by Petersen (1986) and then
slightly modified and incorporated in the OCD model (Hanna and Dicristofaro, 1988). These
algorithms use the platform parameters to enhance the initial dispersion of the emission plume
in a range of distances based on the distance from the emission source and the total height of the
buildings and platform relative to the height of the emission release.


       In Petersen’s original formulation, additional initial dispersion added due to the presence
of a platform in the lateral (y) and vertical (z) dimension is defined when 2.2 < 𝑥𝑥⁄𝐻𝐻𝑏𝑏 < 12.6:

                                   ′
                                𝜎𝜎𝑦𝑦0 = [1.9 + 48.2(2𝑥𝑥 ⁄𝑊𝑊 )−1.4 ]0.5

                                                     140
                                    ′
                                 𝜎𝜎𝑧𝑧0 = [3.0 + 40.2(𝑥𝑥⁄𝐻𝐻𝑏𝑏 )−1.4 ]0.5

       where 𝑊𝑊 is the platform width, 𝐻𝐻𝑏𝑏 is the total height of the platform which is defined as
the height above the sea surface of the top of highest influential building, and 𝑥𝑥 is the distance
downwind from the source. When 𝑥𝑥⁄𝐻𝐻𝑏𝑏 is outside this range, use 𝑥𝑥⁄𝐻𝐻𝑏𝑏 = 2.2 or 𝑥𝑥⁄𝐻𝐻𝑏𝑏 = 12.6
to determine the initial dispersion.


       Modifications to Petersen’s original platform downwash formulation were needed to
make them applicable in all conditions for which the OCD model may be applied. The additional
initial dispersion added due to the presence of a platform in the lateral (y) and vertical (z)
dimension is defined when 2.2 < 𝑥𝑥⁄𝐻𝐻𝑏𝑏 < 12.6:


                               ′
                            𝜎𝜎𝑦𝑦0 = 0.071𝑥𝑥[1.9 + 48.2(2𝑥𝑥 ⁄𝑊𝑊 )−1.4 − 1]0.5
                               ′
                            𝜎𝜎𝑧𝑧0 = 0.11𝑥𝑥 0.81 [3.0 + 40.2(𝑥𝑥⁄𝐻𝐻𝑏𝑏 )−1.4 − 1]0.5



The additional initial dispersion due to the platform was added in quadrature to the lateral and
vertical dispersion calculated without the platform.


       AERMOD considers multiple plumes to derive a total solution for pollutant impacts
under certain conditions. In a convective environment, AERMOD will define as few as one
plume when the release is injected above the mixed layer or as many as three distinct plumes
when the release is below the mixed layer which include the direct plume, the indirect plume,
and the penetrated plume. In addition, the direct and indirect plumes each have distinct updraft
and downdraft components. Separate lateral dispersion coefficients (σy) are computed for each
of the three plumes. Separate vertical dispersion coefficients (σz) are computed for each of the
updraft and downdraft components of the direct and indirect plumes and the penetrated plume
when defined. In a stable environment, a single plume is also defined.



                                                      141
       In AERMOD σy and σz are computed as the root mean square of multiple components,
commonly referred to addition in quadrature, including ambient turbulence and Buoyancy
Induced Dispersion (BID). For sources subject to platform downwash, lateral and vertical
dispersion components due to downwash have been incorporated into the appropriate root mean
square equations. The downwash component was added to the solution for σy for a plume in a
stable environment as well as the direct and penetrated plumes in an unstable environment. The
downwash component was added to the solution for σz for a plume in a stable environment, the
updraft and downdraft components of the direct plume in a convective environment, and a
penetrated plume in a convective environment when defined.


       The integration of the OCD platform downwash algorithm into AERMOD requires only
three additional inputs beyond those required for POINT, POINTHOR, and POINTCAP source
types to define the platform. These inputs generally mimic those required by the OCD model
with some modification to their definitions, for adaption to AERMOD. These inputs, as well as
stack height, as assumed and interpreted by AERMOD, are described below:


       Stack height: Height of the release point above the sea surface. NOTE: This is input
       with the POINT, POINTCAP, or POINTHOR keyword on the SO card.

       Building height (𝐇𝐇𝐛𝐛 ): Height above the sea surface of the top of the tallest structure that
       could influence downwash. NOTE: This is input with the PLATFORM keyword on the
       SO card.

       Platform width (𝐖𝐖): Assuming a rectangular platform, the shorter of the two sides when
       comparing the lateral distance across the left most and right most buildings or structures
       on the platform that can contribute to downwash. NOTE: This is input with the
       PLATFORM keyword on the SO card.

       Platform base height (Zelp): The height of the bottom of the platform above the sea
       surface. NOTE: This is input with the PLATFORM keyword on the SO card.


                                                    142
         There are several differences between the OCD and AERMOD models with respect to
default model behaviors. AERMOD, by default and part of the regulatory formulation of
AERMOD, applies stack-tip downwash, gradual plume rise, and buoyancy induced dispersion
while OCD while none of these are part of the regulatory formulation of OCD. Thus, by default,
AERMOD will apply stack-tip downwash, gradual plume rise, and buoyancy induced dispersion
for a platform source consistent with all other POINT, POINTHOR, and POINTCAP source
types.


B.2.7.4 SWPOINT


         Based on studies performed by researchers at Cornell University, a new source type,
SWPOINT, was added beginning with AERMOD version 21112 to further study the
phenomenon researchers have coined as “sidewash” which is a lateral shift of the building wake
cavity that forms on the lee side a building that occurs when the wind is oblique to one of the
longer sides of an elongated building. Yang et al. (2020) investigated the flow structures and
concentration fields under oblique wind conditions. A key finding was that the flow is not only
entrained downward (downwash) but also directed by the oblique wind along the building
leeward surface (sidewash), which creates a sidewash-downwash (S-D) vortex. This S-D vortex
causes the plume to shift in the lateral direction.


         As a research tool, the SWPOINT source has several limitations that need to be
highlighted for the user who uses the SWPOINT source type to investigate and study the
sidewash effect. These limitations include the following:


         •    Produces downwind concentrations from short stacks assumed to be centered along
              the leeward side of elongated buildings. Refer to Figure X below for an illustration
              of the assumed stack location with respect to the building and wind flow
              orientation.
         •    Model concentrations are limited to the building wake cavity. Impacts at receptors
              outside or in the transition zone from the building wake cavity are not calculated.

                                                      143
       •       Terrain impacts are not considered.
       •       Plume rise due to mechanical and thermal buoyancy is not considered.
       •       Building representation is assumed to be rectangular.
       •       PRIME downwash is not applied.
       •       SWPOINT sources have not been configured for use with EVENT processing,
               NOX-to-NO2 conversion methods, or the MAXDCONT source culpability
               processing.
       •       Stack top wind speeds below approximately 2 m/s have been shown to cause
               anomalously high concentrations. This parameter is quite sensitive and currently
               can shift from reasonable predicted concentrations to values several orders of
               magnitude with the adjustment of a few tenths of a meter per second.




        Figure B – 2. Fixed Stack location with respect to Building and Wind Flow Orientation



     Low Wind Conditions (LOW_WIND)

       To further explore the performance of AERMOD during low wind conditions and
investigate possible improvements to AERMOD performance, the primary keyword
LOW_WIND is available on the CO pathway to override default values for several variables and
default behavior. The EPA has identified two areas of model improvements to address these low
wind issues:


   1. Plume meander: adjustments to the treatment of plume meander within AERMOD and
      the modifications to the meander components, and


                                                     144
   2. Meteorology and turbulence: adjustments to the minimum turbulence values, specifically
      the minimum sigma-v value and the sigma-w value, which are components of the
      turbulence, as well as the minimum wind speed, which is closely tied to the minimum
      turbulence values.

Using the LOW_WIND primary keyword, users can override the following:


   •   Adjustments to met

          o Minimum sigma-v - standard deviation of horizontal velocity, a turbulence
            parameter that impacts the horizontal size of the plume. Estimates of sigma-v in
            low wind conditions are often too small (approaching zero), requiring a minimum
            value to be set. During these low wind conditions, plume volumes are inherently
            small, generally resulting in higher concentrations. Increasing the minimum
            sigma-v will result in lowering the maximum concentrations estimated for surface
            releases.

          o Minimum wind speed – limits the calculated wind speed computed in the wind
            profile.

                    The model formulation has the minimum wind speed correlated to the
                     minimum sigma-v.

                    AERMET uses the 2*sigma-v to reset low winds, this is applied more
                     broadly in AERMOD as a check against effective wind speeds (i.e., the
                     AERMET restriction is applied at the meteorology measurement height,
                     but applied at different heights on the profile in AERMOD)
                                                            2
                                              𝑢𝑢         = �2 ∗ 𝜎𝜎𝑣𝑣2
                                                   min

          o Minimum sigma-w - standard deviation of vertical velocity, a turbulence
            parameter that impacts the vertical size of the plume. Estimates of sigma-w in low
            wind conditions are often too small (approaching zero), requiring a minimum
            value to be set. During these low wind conditions, plume volumes are inherently
            small, generally resulting in higher concentrations. Increasing the minimum
            sigma-w will result in lowering the maximum concentrations estimated for
            surface releases.

                                                   145
   •   Adjustments to dispersion

           o Plume meander/Lower and upper limit of FRAN – The Fraction of the Random
             plume, or FRAN, is the weighting factor used to average the direct or coherent
             plume and the random or pancake plume. The default range of FRAN in
             AERMOD is 0.0 - 1.0. These user-defined adjustments were added to off-set
             other adjustments to the meteorology result from the low wind options. Lowering
             the value of FRAN in general will increase concentrations.
                                                              𝜎𝜎𝑟𝑟2            𝜎𝜎𝑟𝑟2
                                       𝐶𝐶𝑐𝑐,𝑠𝑠 = 𝐶𝐶𝐶𝐶ℎ �1 −         � + 𝐶𝐶𝑅𝑅 �       �
                                                              𝜎𝜎ℎ2             𝜎𝜎ℎ2
           o Big T – The random plume time scale, of Big T, is a scaling parameter that
             indicates the time lapse for which mean wind information at the source is no
             longer correlated with the location of plume material at a downwind receptor.
             Changing this term will increase the FRAN for any particular situation, but the
             impact is also distance-dependent, so the computed value will be larger at greater
             distances (i.e., the random portion of the plume will be larger at greater
             distances).

                      RLINE does not account for this distance dependence and in general, the
                       meander computation is slightly different in RLINE.
                                                                 −𝑥𝑥𝑟𝑟
                                𝜎𝜎𝑟𝑟2 = 2𝜎𝜎�𝑣𝑣2 + 𝑢𝑢�2 �1 − exp �         ��
                                                                  𝑢𝑢�𝑇𝑇𝑟𝑟
           o Momentum Balance FRAN Method - Alternate momentum balance (PBAL)
             approach to determine plume meander which overrides the default energy balance
             approach.

B.2.8.1 More on FRAN Momentum Balance Method


       AERMOD accounts for plume meander by interpolating between two concentration
limits: the coherent plume limit (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 (assumes an equal probability of any wind direction).


                                                        146
       Once these 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,


                                               𝜎𝜎2            𝜎𝜎2
                       𝐶𝐶𝑐𝑐,𝑠𝑠 = 𝐶𝐶𝐶𝐶ℎ �1 − 𝜎𝜎𝑟𝑟2 � + 𝐶𝐶𝑅𝑅 �𝜎𝜎𝑟𝑟2 �         (Equation B-6)
                                                ℎ              ℎ


       where σh2 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 σr2/σh2 is an indicator of the
importance of the random component and can be used to weight the two concentrations as done
in Equation B-8.


       The horizontal wind is composed of a mean component ū, and random components σu
and σv. Thus, a measure of the total horizontal wind “energy” (given that the along-wind and
crosswind fluctuations are assumed equal i.e., σu = σv), can be represented as:


                                 𝜎𝜎ℎ2 = 2𝜎𝜎�𝑣𝑣2 + 𝑢𝑢�2                (Equation B-7)

       where ū = (𝑢𝑢̃ 2 − 2𝜎𝜎̃𝑣𝑣2) 1/2, where u is the mean component of the wind. The random energy
component is initially 2𝜎𝜎̃𝑣𝑣2 and becomes equal to σh2 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:

                                                            −𝑥𝑥
                       𝜎𝜎𝑟𝑟2 = 2𝜎𝜎�𝑣𝑣2 + 𝑢𝑢�2 �1 − exp � 𝑢𝑢�𝑇𝑇𝑟𝑟 ��         (Equation B-8)
                                                               𝑟𝑟



where Tr is a time scale (= 24 hours) at which mean wind information at the source is no longer
correlated with the location of plume material at a downwind receptor.

                                                              147
       The current blending between the coherent and random plume in AERMOD is based on
the relative amount of energy indicated by the mean horizontal wind versus the energy indicated
by the standard deviation in the horizontal winds. Recall that the formula for kinetic energy
(KE) is:

                                      1
                               𝐾𝐾𝐾𝐾 = 2 𝑚𝑚 ∗ 𝑣𝑣 2              (Equation B-9)



where m is mass and v is velocity.


       Normally, the wind speed is a vector quantity. However, similarity theory, which is used
to compute the AERMOD wind profile, treats the wind speed and wind direction separately,
such that the wind speed is considered a scalar base-quantity. That is, a composite wind speed
and wind direction are used rather than explicit values of the horizontal wind components (i.e., a
north-south component and an east-west component). Accounting for wind speeds in the manner
could result in an average wind speed of zero simply because the winds were from opposite
directions in each part of the hour. The energy equation takes the dot-product of the wind
speeds, resulting a scalar quantity as well. The result of using an energy-based balancing of the
two wind speed values is that differences between the two are emphasized, particularly when
one value is greater than one and the other is less than one. Under lower wind speeds, the
horizontal fluctuations are typically larger (low winds are often characterized by frequent and
wind changes in the wind direction, even though the wind speeds may be low). As wind speeds
increase, the horizontal fluctuations typically decrease. While the horizontal fluctuations are all
generally less than 1 m/s (in fact, horizontal fluctuations are typically less than 0.5 m/s), the
wind speed is typically greater than 1 m/s. When the winds increase, the use of an energy
balance, which is a quadratic function of the wind speeds, will rapidly decreases the fraction of
the random plume.




                                                     148
       An alternative to using an energy balance would be to use a momentum balance. Recall
that the equation for momentum (P) is:


                                    𝑃𝑃 = 𝑚𝑚 ∗ 𝑣𝑣    (Equation B-10)

The use of a momentum also differs from the energy approach in several ways:



   •   The momentum will vary linearly rather than a quadratic function of the wind speeds. As
       a result, divergence in the values will be less emphasized than in the current formulation.
       The result will be that fraction of the random plume will generally be larger at higher
       wind speeds.


   •   In contrast to the energy balance, a momentum balance would be a vector-based quantity,
       which is in some ways is more consistent with the concept of transport of a pollutant via
       winds, which will have a direction dependency (i.e., a vector dependency).


The momentum-based approach is applied as the square root of the random plume fraction.




                                                   149
      PSD Credit (PSDCREDIT)

       Due to the ozone limiting effects of the PVMRM option, the predicted concentrations of
NO2 are not linearly proportional to the emission rate. Therefore, the approach of modeling
NO2 increment consumption with PSD credits through the use of a negative emission rate for
credit sources cannot be used with the PVMRM option.


       Increment is the maximum allowable increase in concentration of a pollutant above a
baseline concentration for an area defined under the Prevention of Significant Deterioration
(PSD) regulations. The PSD baseline area can be an entire State or a subregion of a State such
as a county or group of counties. Increment consumption is the additional air quality impact
above a baseline concentration.


       The baseline concentration is the ambient concentration of the pollutant that existed in
the area at the time of the submittal of the first complete permit application by any source in that
area subject to PSD regulations. A baseline source is any source that existed prior to that first
application and the baseline date is the date of the PSD application. This baseline date is
referred to as the minor source baseline date in PSD regulations. By definition, baseline sources
do not consume increment. However, any baseline source that retires from service after the
baseline date expands the increment available to new sources. Therefore, a PSD modeling
analysis performed for a new source may need to account for this increment expansion. Such an
analysis may therefore involve identification of three groups of sources: 1) increment-
consuming sources; 2) retired (increment-expanding) baseline sources; and 3) existing, non-
retired, baseline sources.


       Calculating increment consumption under the PSDCREDIT option in AERMOD is not a
simple arithmetic exercise involving the three groups of sources defined above. Since the
amount of ozone available in the atmosphere limits the conversion of NO to NO2, interactions
of plumes from the existing and retired baseline sources with those from the increment
consuming sources must be considered as part of the calculation of net increment consumption.
                                                 150
Without the PSDCREDIT option, properly accounting for the potential interaction of plumes
among the different source categories would require post-processing of results from multiple
model runs. Internal processing algorithms have been incorporated in AERMOD under the
PSDCREDIT option to account for the apportioning of the three groups of sources to properly
calculate increment consumption from a single model run.


       Defining the following three source groupings for the discussion that follows:


               A = increment-consuming sources;
               B = non-retired baseline sources; and
               C = retired baseline, increment-expanding sources.
       The calculation of the amount of increment consumption by the A sources cannot simply
be estimated by modeling the A sources alone because of the possible interaction of those
plumes with the plumes from B sources. The PVMRM algorithm is designed to account for
such plume interactions and calculate the total NO to NO2 conversion in the combined plumes
based on the amount of ozone available. Therefore, the total increment consumption by the A
sources is given by the difference between (1) the total future impact of increment consuming
sources and non-retired baseline sources (A+B) and (2) the total current impact (B), which can
be expressed as (A+B) – (B). Here (A+B) represents the value that would be compared against
the National Ambient Air Quality Standard (NAAQS) for NO2 during PSD review of the A
sources.


       In a case where some of the baseline sources have been retired from service (C sources),
the PSD regulations allow the consideration of increment expansion when assessing compliance
with the PSD increment. However, the amount of increment expansion cannot be estimated by
simply modeling the C sources alone because of the possible interaction of those plumes with
the plumes from B sources. Therefore, the total increment expansion, i.e., PSD credit, is
calculated as the difference between (1) the total impact prior to the retirement of C sources, i.e.


                                                    151
(B+C), and (2) the total impact from existing (non-retired) baseline sources (B), which can be
expressed as (B+C) – (B).


       Finally, the net increment consumption is given by the difference between total
increment consumption and the total increment expansion, or


                       [(A+B) – (B)] – [(B+C) – (B)]          (Equation B-11)


       Note that in the absence of any increment expansion, the net increment consumption is
equal to the total increment consumption [(A+B) – (B)], as described above.


       These expressions of net increment consumption and expansion cannot be interpreted as
algebraic equations. Instead, the terms within parentheses represent the results of separate
model runs that account for the combined effects of NOx conversion chemistry on specific
groups of sources. The expression shown in Equation B-11 above represents four model
simulations: (A+B), (B), (B+C), and (B) again. In this case, the two (B) terms do cancel each
other, and we are left with:


                                      [(A+B)] – [(B+C)]       (Equation B-12)


       The expression presented in Equation B-12 summarizes how the net increment
consumption calculation is performed under the PSDCREDIT option. Under this option,
AERMOD first models the A and B groups together, then models the B and C groups together,
and finally computes the difference to obtain the desired result, i.e., the value to compare to the
PSD increment standard. For AERMOD to perform the special processing associated with this
option, the user must define which sources belong to each of the groupings defined above.




                                                    152
       Deposition Options

       The AERMOD model includes algorithms for both dry and wet deposition of both
particulate and gaseous emissions. The deposition algorithms incorporated into AERMOD are
based on the draft Argonne National Laboratory (ANL) report (Wesely et al., 2002), with
modifications based on peer review. Treatment of wet deposition was revised from Wesely et
al. (2002) based on recommendations by peer review panel members (Walcek et al., 2001). A
full technical description of the deposition algorithms implemented in AERMOD is provided in
an EPA report specific to these algorithms (U.S. EPA, 2003b).


       Based on the guidance provided for application of the AERMOD model in the Guideline
(U.S. EPA, 2024a), and the history of the deposition algorithms in the AERMOD and ISC
models, the particle deposition algorithms with a user-specified particle size distribution
(referred to below as “Method 1”) can be applied simultaneously with the regulatory DFAULT
keyword. Method 1 is comparable to the particle deposition algorithm in the ISCST3 model
(EPA, 1995a). The gas deposition algorithms and the “Method 2” option for particle deposition
based on the ANL draft report (Wesely, et al, 2002) are non-regulatory ALPHA options in
AERMOD, and beginning with version 19191, the model will issue a fatal error message and
abort processing if the ALPHA keyword is not specified with the gas deposition or Method 2
particle deposition options.




                                                    153
                                     Table B - 1. Summary of Deposition Options
Pollutant Type Model Output Type               Required Keywords        Allowed under DFAULT?
Gaseous            CONC w/dry depletion        CO GASDEPVD                           No1
                   DDEP                         or
                                               CO GDSEASON,
                                               CO GDLANUSE, and
                                               SO GASDEPOS
Gaseous            CONC w/wet depletion        SO GASDEPOS                           No1
                   WDEP
Gaseous            CONC w/dry & wet            CO GDSEASON,                          No1
                     depletion                 CO GDLANUSE, and
                   DEPOS                       SO GASDEPOS
Particulate        CONC w/dry and/or wet       SO PARTDIAM,                          Yes2
(“Method 1”)         depletion                 SO PARTDENS, and
                   DEPOS                       SO MASSFRAX
                   DDEP
                   WDEP
Particulate        CONC w/dry and/or wet       SO METHOD_2                           No1
(“Method 2”)         depletion
                   DEPOS
                   DDEP
                   WDEP
1
    The ALPHA option must be included.
2
 While “Method 1” is allowed under the regulatory “DFAULT” option within AERMOD, the use of
“Method 1” for particulate emissions in regulatory modeling applications should follow the guidance
provided in Section 7.2.1.3 of the Guideline (U.S. EPA, 2024a).



          The user should be aware that one or more of the following meteorological parameters
are needed for deposition: precipitation code, precipitation rate, relative humidity, surface
pressure, and cloud cover.




                                                          154
       Time Travel Reaction Method (TTRM/TTRM2)

       Beginning with version 21112, the Travel Time Reaction Method (TTRM) was added as
an ALPHA option for NO2 conversion that considers the distance and the travel time from the
emission source to each receptor. In general, much of the conversion of NO to NO2 occurs within
the first minute of travel which limits the effectiveness of this method to the near field receptors.
PVMRM, OLM, and ARM2 assume the reaction involving NO and available O3 to form NO2
occurs instantaneously. Although this chemical reaction is relatively rapid, it is not actually
instantaneous and depends on the transport time to the downwind receptor of interest. TTRM,
was initially implemented as a stand-alone ALPHA option which can determine the initial
fraction of NO to NO2 conversion in the travel time of each source emissions to each receptor.
The conversion is capped at an upper limit, which is typically reached after a few tens of seconds
of plume travel. Beyond the distance the fraction reaches the upper limit of the equilibrium
fraction (generally 0.9), TTRM is no longer effective, and another method is needed for
receptors beyond that distance.
       Beginning with version 22112, TTRM was integrated more fully into AERMOD to be
used simultaneously with PVMRM, OLM, or ARM2. This integration was also added as the
ALPHA option, TTRM2, separate from TTRM which was retained as a stand-alone option.
When TTRM2 is specified along with PVMRM, OLM, or ARM2, TTRM will be implemented
for near field receptors where the fraction of conversion has not reached the upper limit, and the
other specified method will be used for all other receptors.
       Though the NOX chemistry regime can involve a larger number of reactions, it is fairly
well understood and can be parameterized to include a wide range of the relevant reactions.
However, in the near-field (i.e., time scales on the order of minutes), the dominant reactions
during the daytime are between NO, NO2, and O3 include:


               O3 + NO  NO2 + O2 (k1)                (Reaction 1)

               O2 + NO2 + hv  NO + O3 (k2)           (Reaction 2)


                                                    155
        where k1 and k2 are the reaction rate constants for each reaction and hv is the incident
solar radiation.3
        These reactions (i.e., Reactions 1 and 2) are often summarized to emphasize the main
reaction pathways commonly described as the pseudo-steady state (PSS) approximation. The
steady-state assumption is based on the fact that both reactions are relatively and equally fast so
equilibrium is reached quickly such that the reaction rate is the product of the reaction rate
constant (i.e., k1 or k2) and the concentration of the reactants. So, when steady state is assumed,
the reaction rates are equal and can be expressed as:


                k1[O3][NO] = k2[NO2]                          (Equation B-13)

        where k1 = (15.33/T) exp(–1450/T) ppb–1 sec–1, k2 = 0.0167exp(–0.575/cos(θ)) sec–1, θ
= zenith angle of sun (function of location latitude and time of day), and T is the ambient
temperature in Kelvin (see equations 10 and 11 in Hanrahan, 1999). In the context of modeling
with AERMOD, where NO is released and NOX is conserved, NO can be removed from
Equation B-13 by using the following relationship:


                NO = NOX – NO2                                (Equation B-14)


        Substituting Equation B-14 into Equation B-13 and solving for NO2/NOX gives the
following relationship:


                [NO2] / [NOX] = (k1[O3]) / (k2 + k1[O3])      (Equation B-15)


        Based on Equation B-15, the PSS ratio of NO2/NOX is a function of O3, temperature, and
sunlight. It should be noted that this solution assumes O3 has reached PSS and thus, the reaction
is not O3 limited. The solution of this equation can provide a theoretical maximum for the
NO2/NOX ratio.



                                                    156
       The time-component (t) of this reaction can be computed based on reaction 1 and the
value of k1, giving the expression:


                d[NO]/dt = k1[NO][O3]               (Equation B-16)


       Following the rearrangement and integration of Equation B-16, the relationship becomes:


                ln [NO] = - k1[O3]t + C             (Equation B-17)


where C is some constant. The constant, C, can be determined with the boundary conditions,
when t = 0 and [NO] equals the starting NO concentration, [NO]o. In this case, C = ln [NO]o and
the integrated form for [NO] as a function of time can be expressed as:


                ln [NO] = - k1[O3]t + ln [NO]o      (Equation B-18)


       Following the rearrangement of Equation B-18, the fraction of NO remaining can be
expressed as:


                [NO]/[NO]o = exp (-k1[O3]*t)         (Equation B-19)


       Then the difference of the remaining NO converted to create NO2 can be expressed as
the fraction of NO converted to create NO2 or NOfrac:


                NOfrac = 1 – exp(-k1[O3]*t).        (Equation B-20)

       As shown by Hanrahan (1999) and detailed above, the reaction rate, k1, can be modified
by the ambient temperature (Tamb, degree Kelvin). Equation B-20 provides an expression for
calculating the travel time reaction of the NO-to-NO2 conversion by O3 in modeling systems.




                                                  157
       Figure B-3 illustrates the fractional conversion of emitted NO to NO2 (i.e., NOfrac) as a
function of O3 and travel time for three examples of ambient O3 levels, assuming that the
ambient temperature is 298 K in equation 9 and then substituting k into Equation B-20. Note
that the fraction of NO converted NO2 in the atmosphere needs to account for the in-stack NO2
as well. For example, if the in-stack NO2 is 10% of the emissions and 50% of the remaining NO
is converted to NO2, then the total NO2 fraction is 0.10 + 0.5*0.9 = 0.55.




        Figure B - 3. Conversion fraction of NO-to-NO2 as a function of O3 and travel time.


       It is evident from Figure B-3 that after one minute of travel, much of the NO (for which
there is available ozone) has been converted to NO2, except for cases with low ambient ozone
levels where much of the NO has been converted to NO2 after two minutes of travel. However,
there could be cases with well under a minute or two of travel time for which the reaction time
conversion fraction is an important consideration. It is to be noted that this forward reaction (to
convert NO to NO2) is complemented later by an equilibrium state with another photochemical
reaction that converts NO2 to NO in the opposite direction in the presence of sunlight.


                                                    158
Therefore, the converted NO2 fraction noted above is not an equilibrium result, but it can act as
a lower limit to the conversion in the first minute after fresh NO is emitting from a source.


       The TTRM option determines the initial fraction of NO-to-NO2 conversion (i.e., NOfrac)
in the travel time of each source emissions to each receptor. As shown in Figure B-3 and
Table B-2, this option is most important for refining NO2 modeled impacts under short travel
times and within minutes of emission. Further, the effect of this option is most evident at near-
field receptors or at distances within a kilometer of the emission source (see Table B-2). The
method used in TTRM also assumes that the conversion is capped at an upper limit, which is
typically reached within minutes of plume travel. Once the fraction reaches the upper limit of
the equilibrium fraction (generally 0.9) with the calculation noted above, this reaction for further
travel times is no longer an issue. Beyond these distances, other factors would play a greater
role in NO conversion. Therefore, the implementation of this method could be considered as a
refinement for Tier 3 options, rather than a stand-alone option in AERMOD, with the TTRM
method providing the approach for NO2 concentrations in the near-field and the Tier 3 options
dominant elsewhere in the modeling domain.


   Table B - 2. Conversion Fraction of NO to NO2 as a function of Ozone Concentration,
                            Transport Time, and Wind Speed
                  Conversion Fraction of NO to NO2             Distance from Source [m]
    Travel Time
                  Ozone =    Ozone =      Ozone =    Wind Speed = Wind Speed = Wind Speed =
     (seconds)
                  35 ppb      55 ppb       75 ppb       1 m/s            5 m/s          10 m/s
        10         0.130      0.196        0.257           10           50             100
        20         0.242      0.353        0.448           20          100             200
        30         0.340      0.480        0.590           30          150             300
        40         0.426      0.582        0.696           40          200             400
        50         0.500      0.664        0.774           50          250             500
        60         0.565      0.730        0.832           60          300             600
        120        0.811      0.927        0.972           120         600             1200
        180        0.918      0.98         0.995           180         900             1800
        240        0.964      0.995        0.999           240         1200            2400
        300        0.984      0.999        1.000           300         1500            3000



                                                     159
       The information needed to use this method and calculate the travel time reaction of the
NO-to-NO2 conversion in AERMOD is already available in the model, including the hourly
ambient temperature provided in the meteorological inputs and ambient O3 data when Tier 3
NO2 conversion options are used in AERMOD. The method also involves some assumptions to
calculate the travel time reaction, where the travel time reaction is dependent on which plume
type (i.e., coherent or meander plume) has the larger weighting fraction and the hourly effective
wind speed at plume height. Plume type is important for determining the algorithm that will be
used for calculating the distance relative to the emission source and receptor location and wind
direction. The distance algorithm can be based on a downwind, radial, or crosswind distance in
AERMOD (see Figure B-4). For point and volume sources, the downwind distance is used with
coherent dominant plumes at a receptor and radial distance is used with meander dominant
plumes. The radial distance from the center of the emission source is used for all area sources.
The crosswind distance does not factor in the TTRM calculations. AERMOD also calculates
effective wind speeds based on plume type, where calculations may use the standard effective
wind speed, penetrated effective wind speed, or direct effective wind speed.




Figure B - 4. Distance Terms in AERMOD Relative to Source and Receptor Location and
                                   Wind Direction.




                                                   160
       Once the most appropriate distance and effective wind speeds are determined, the
reaction rate can be calculated and used to determine the NOfrac. As noted above, this fractional
value is capped at the NO2 Equilibrium level that by default in AERMOD is set to 90 percent.
As with the Tier 3 methods, an in-stack NO2 contribution is calculated from the hourly NO value
and the remaining available NO is then scaled by NOfrac and added to the in-stack contribution to
obtain the total NO2 from each source at each receptor. The next step involves selecting the
lower predicted NO2 impact for each source and receptor. Thus, for receptors beyond the
distance needed for the TTRM method to achieve a 90 percent conversion of NO to NO2, the
standard Tier 2 or 3 method would dictate the result. At closer distances, the lower fraction
between the TTRM method and the Tier 2 or 3 method would dictate the result.

       Highly Buoyant Plume (HBP)

       Beginning with AERMOD version 23132, the highly buoyant plume (HBP) option was
added as an ALPHA option. The HBP option was developed to refine AERMOD’s treatment of
hot buoyant plumes that penetrate the top of the convective mixed layer into the stable layer
above. The formulation is based on the work of Weil (2020) and further refined by Warren et
al. (2022) who have demonstrated an overprediction of concentrations at the surface for some
hours due to what the authors describe as a premature mixing of the penetrated plume back into
the mixed layer and down to the surface. Observations show that the penetrated plume resides
in the stable layer until, through surface heating, the depth of the mixed layer grows until the top
of the mixed layer intercepts the plume (Moore, 1988; Weil, 1997).


       In short, Warren et al. (2022) state that the current formulation of AERMOD overstates
vertical plume depth by overestimating sigma-z of the penetrated plume which represents the
vertical dispersion of the plume. This assumes that a portion of the penetrated plume extends
down into the mixed layer prior to when the mixed layer would reach a height to intersect the
plume, resulting in premature mixing of the plume to the ground. The revised approach,
implemented as the HBP alpha option, compares the convective mixing height of the current
hour and next hour to determine how much the convective layer has grown over the hour and

                                                    161
the amount of the plume that has been mixed into the convective layer by the end of the current
hour which determines the amount of plume allowed to mix down to the ground (Warren et al.,
2022).


         Area Meander (AREAMNDR)

         Beginning with AERMOD version 23132, plume meander has been added as an ALPHA
option for area sources (AREAMNDR), including AREA, AREAPOLY, AREACIRC, and
LINE source types. This option applies the plume meander algorithms used with POINT and
VOLUME source types to the area source types listed above, but with limitations.


         Generally, model runtimes associated with area sources (without meander applied), can
be significantly longer than those of POINT and VOLUME, because of the double integral
formulation and numerical integration computations coded for resolving area source plume
dimensions in the alongwind and crosswind directions. Activating plume meander for area
sources, therefore, doubles the number of calculations and computational demands needed to
resolve area source dimensions for both the coherent and meander plumes. As such, processing
time is increased by approximately 30-50% when meander is applied to area sources.


         The current implementation of area meander as an ALPHA option limits the receptors
for which meander is applied. Moreover, meander is not being calculated for upwind receptors
which affects both upwind and downwind concentrations; this is inconsistent with point and
volume source meander. In addition, there remain scientific concerns with applying the
meander algorithms developed for POINT and VOLUME source types to area source types.
Additional research is required.




                                                  162
       Aircraft Plume Rise (ARCFTOPT)

       An ALPHA option to treat plume rise from aircraft emissions, sponsored by the Federal
Aviation Administration (FAA), was added to AERMOD beginning with version 23132.
Aircraft emissions from jet engines experience plume rise from both momentum and buoyancy
but are commonly modeled as AREA and VOLUME source types in AERMOD which do not
account for either momentum or buoyancy. The aircraft plume rise ALPHA option extends the
formulation of AREA and VOLUME sources with additional input parameters based on the
work of Pandey, et. al. (2023). AREA and VOLUME sources that include hourly varying
aircraft input parameters can be identified through the source block keyword ARCFTSRC. The
current implementation of aircraft plume rise requires that an additional seven aircraft
parameters be supplied to AERMOD using an hourly varying emissions input file.




                                                   163
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.
Ahangar, F.E., Heist, D.K., Perry, S.G., and Venkatram, A. (2017). Reduction of air pollution
      levels downwind of a road with an upwind noise barrier. Atmospheric Environment, 155,
      pp. 1-10.
Andre, J. C. and L. Mahrt, 1982: The nocturnal surface inversion and influence of clean-air
       radiative cooling. J.Atmos.Sci., 39, 864-878.
Azzi M. and Johnson G (1992). An Introduction to the Generic Reaction Set Photochemical
      Smog Mechanism. Proc. 11th Clean Air Conf. 4th Regional IUAPPA Conf., Brisbane,
      Australia.
Baerentsen, J. H. and R. Berkowicz, 1984: Monte Carlo simulation of plume dispersion in the
       convective boundary layer. Atmos.Environ., 18, 701-712.
Baldauf, R., Thoma, E., Hays, M., Shores, R., Kinsey, J., Gullett, B., Kimbrough, S., Isakov, V.,
       Long, T., Snow, R. and Khlystov, A. (2008). Traffic and meteorological impacts on near-
       road air quality: Summary of methods and trends from the Raleigh near-road study.
       Journal of the Air & Waste Management Association, 58(7), pp.865-878.
Baldauf, R.W., Heist, D., Isakov, V., Perry, S., Hagler, G.S., Kimbrough, S., Shores, R., Black,
       K. and Brixey, L. (2013). Air quality variability near a highway in a complex urban
       environment. Atmospheric environment, 64, pp.169-178.
Baldauf, R.W., Isakov, V., Deshmukh, P., Venkatram, A., Yang, B. and Zhang, K.M. (2016).
       Influence of solid noise barriers on near-road and on-road air quality. Atmospheric
       Environment, 129, pp.265-276.
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.
Benson, P.E. (1992). A review of the development and application of the CALINE3 and 4
      models. Atmospheric Environment. Part B. Urban Atmosphere, 26(3), pp.379-390.




                                                   164
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.
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, ATDL-
       106, 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.
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.


                                                  165
Carruthers, D. J., and Coauthors, 1992: UK atmospheric dispersion modelling system. Air
       Pollution Modeling and Its Application. Plenum Press, New York.
Carruthers, D.J., Stocker, J.R., Ellis, A., Seaton, M.D. and Smith, S.E. (2017). Evaluation of an
       explicit NOX chemistry method in AERMOD. Journal of the Air & Waste Management
       Association, 67(6), pp.702-712.
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 short-
       term 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,
        390-392.
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.




                                                   166
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.
Fairall, C.W., E.F. Bradley, J.E. Hare, A.A. Grachev, and J.B. Edson, 2003: Bulk
         Parameterization of Air-Sea Fluxes: Updates and Verification for the COARE Algorithm.
         J. Climate, 16, 571-591.
Finn, D., Clawson, K.L., Carter, R.G., Rich, J.D., Eckman, R.M., Perry, S.G., Isakov, V. and
       Heist, D.K. (2010). Tracer studies to characterize the effects of roadside noise barriers on
       near-road pollutant dispersion under varying atmospheric stability conditions.
       Atmospheric Environment, 44(2), pp.204-214.
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.
Hanna, S.R., Schulman, L.L., Paine, R.J., Pleim, J.E., Baer, M., 1985. Development and
       evaluation of the offshore and coastal dispersion model. J. Air Pollut. Control Assoc. 35,
       1039–1047. https://doi.org/10.1080/00022470.1985.10466003
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.
Hanna, S.R., Dicristofaro, D.C., 1988. Development and Evaluation of the OCD/API (Offshore
       and Coastal Dispersion / American Petroleum Institute) Model.
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. 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.

                                                   167
Hanrahan, P.L., 1999. “The plume volume molar ratio method for determining NO2/NOx ratios
      in modeling. Part I: Methodology,” J. Air & Waste Manage. Assoc., 49, 1324-1331.
Hashad, K., Steffens, J.T., Baldauf, R.W., Heist, D.K., Deshmukh, P., Zhang, K.M., 2024.
      Resolving the effect of roadside vegetation barriers as a near-road air pollution mitigation
      strategy. Environ. Sci.: Adv. 3, 411–421. https://doi.org/10.1039/D3VA00220A.
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.
Heist, D.K., Perry, S.G. and Brixey, L.A. (2009). A wind tunnel study of the effect of roadway
        configurations on the dispersion of traffic-related pollution. Atmospheric Environment,
        43(32), pp.5101-5111.
Heist, D., Isakov, V., Perry, S., Snyder, M., Venkatram, A., Hood, C., Stocker, J., Carruthers, D.,
        Arunachalam, S. and Owen, R.C. (2013). Estimating near-road pollutant dispersion: A
        model inter-comparison. Transportation Research Part D: Transport and Environment,
        25, pp.93-105.
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.
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.


                                                   168
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.
Misra, P. K., 1982: Dispersion of non-buoyant particles inside a convective boundary layer.
       Atmos.Environ., 16, 239-243.
Monbureau, E. M., Heist, D. K., Perry, S. G., Brouwer, L. H., Foroutan, H., Tang, W. (2018).
     Enhancements of AERMOD’s building downwash algorithms based on wind tunnel and
     Embedded-LES modeling. Atmospheric Environment, 179, 321-330.
Moore, G. E., L. B. Milich, and M. K. Liu. 1988. Plume behaviors observed using lidar and SF6
      tracer at a flat and hilly site. Atmos. Environ 22 (8):1673–88. 0004-6981.
      doi:10.1016/0004-6981(88)90396-4.
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.



                                                   169
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.
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.
Pandey, G., A Venkatram, and S. Arunachalam, 2023: “Accounting for plume rise of aircraft
      emissions in AERMOD.” Atmospheric Environment, 314.
      https://doi.org/10.1016/j.atmosenv.2023.120106.
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.




                                                  170
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.
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.
Perry, S.G., Heist, D.K., Brouwer, L.H., Monbureau, E.M., and L.A. Brixley (2016).
        Characterization of pollutant dispersion near elongated buildings based on wind tunnel
        simulations, Atmospheric Environment, Vol. 42, 286-295.
Petersen, R.L., 1986. Wind tunnel investigation of the effect of platform-type structures on
       dispersion of effluents from short stacks. J. Air Pollut. Control Assoc. 36, 1347–1352.
       https://doi.org/10.1080/00022470.1986.10466185
Petersen, R.L., A. Kolesnikov, and A. Beyer-Lout. 2014. Evaluation of new wind tunnel,
       AERMOD and CFD methodologies to determine EBD/BPIPPRM in light of the 2011
       EPA Memorandum. Technical Paper 33263. Presented at the 107th Annual A&WMA
       Conference & Exhibition, Long Beach, CA, June 24 27.
Petersen, R. L., Sergio A. Guerra & Anthony S. Bova. (2017). Critical Review of the Building
       Downwash Algorithms in AERMOD. J. Air Waste Management Association Vol. 67,
       Issue 8.
Petersen, R. L. and Guerra, S. A., (2018). PRIME2: Development and evaluation of improved
       building downwash algorithms for rectangular and streamlined structures. Atmospheric
       Environment, 173, 67-78.
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., and J.S. Scire, 1980: Boyant Line and Point Source (BLP) Dispersion Model
      User’s Guide. Final Report. Environmental Research & Technology, Inc. P-7304B. July
      1980.



                                                   171
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.
Schulte, N., Snyder, M., Isakov, V., Heist, D. and Venkatram, A. (2014). Effects of solid barriers
       on dispersion of roadway emissions. Atmospheric environment, 97, pp.286-295.
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.
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.
Snyder, M. G., A. Venkatram, D. K. Heist, S. G. Perry, W. B. Petersen, and V. Isakov, 2013:
       RLINE: A Line source dispersion model for near-surface releases. Atmos.Environ., 77(0),
       748-756.
Stocker, J., M. Seaton, S. Smith, J. O’Neill, K. Johnson, R. Jackson, and D. Carruthers (CERC).
       Evaluation of the Generic Reaction Set Method for NO2 conversion in AERMOD. The
       modification of AERMOD to include ADMS chemistry. Cambridge Environmental
       Research Consultants (CERC) Technical Report. August 8, 2023.
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. EPA, 1992: Protocol for Determining the Best Performing Model. EPA-454/R-92-025, U.S.
       Environmental Protection Agency, RTP, NC.
U.S. EPA, 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. EPA, 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.

                                                   172
U.S. EPA, 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. EPA, 2003a: AERMOD: Latest Features and Evaluation Results. EPA-454/R-03-003, U.S.
       Environmental Protection Agency, Research Triangle Park, North Carolina.
U.S. EPA, 2003b: AERMOD Deposition Algorithms Science Document (Revised Draft). U.S.
       Environmental Protection Agency, Research Triangle Park, North Carolina 27711.
U.S. EPA, 2006: Air Quality Criteria for O3 ane Related Photochemical Oxidants, Volumes I
       and II. U.S. Environmental Protection Agency, Office of Research and Development,
       EPA-600/R-05-004aF. February 2006.
U.S. EPA, 2012: User’s Manual AERCOARE Version 1.0. EPA-910-R-12-008. U.S. EPA,
       Region 10, Seattle, WA.
U.S. EPA, 2015: Technical Support Document (TSD) for NO2-related AERMOD modifications.
       EPA-454/B-15-004. U.S. Environmental Protection Agency, Research Triangle Park,
       North Carolina. July 2015.
U.S. EPA, 2016: Technical Support Document (TSD) for AERMOD/BLP Development and
       Testing. EPA-454/B-16-009. U.S. Environmental Protection Agency, Research Triangle
       Park, North Carolina. December 2016.
U.S. EPA, 2023a: Evaluation of Addition of Terrain Treatment to the RLINE Source Type in
       AERMOD. EPA-454/R-23-012. U.S. Environmental Protection Agency, Research
       Triangle Park, North Carolina. October 2023.
U.S. EPA, 2023b: Evaluation of the Implementation of the Coupled Ocean Atmosphere
       Response Experiment (COARE) Algorithms into AERMET for Marine Boundary Layer
       Environments. EPA-454/R-23-008. U.S. Environmental Protection Agency, Research
       Triangle Park, North Carolina. October 2023.
U.S. EPA, 2023c: Evaluation of Prognostic Meteorological Data in AERMOD Overwater
       Applications. EPA-454/R-23-012. U.S. Environmental Protection Agency, Research
       Triangle Park, North Carolina. October 2023.
U.S. EPA, 2023d: Incorporation and Evaluation of the RLINE source type in AERMOD for
       mobile source applications. EPA-454/R-23-011. U.S. Environmental Protection Agency,
       Research Triangle Park, North Carolina. October 2023.
U.S. EPA, 2024a. Guideline on Air Quality Models: Ehancements to the AERMOD Dispersion
       Modeling System, Final Rule. 40 CFR Part 51, Appendix W.
U.S. EPA, 2024b: Technical Support Document (TSD) for Adoption of the Generic Reaction Set
       Method (GRSM) as a Regulatory Non-Default Tier-3 NO2 Screening Option. EPA-
       454/R-24-005. U.S. Environmental Protection Agency, Research Triangle Park, North


                                               173
       Carolina. October 2024.
U.S. EPA, 2024c: User's Guide for the AERMOD Terrain Preprocessor (AERMAP). EPA-
       454/B-24-008. U.S. Environmental Protection Agency, Research Triangle Park, NC.
U.S. EPA, 2026a User’s Guide for the AERMOD Meteorological Preprocessor (AERMET).
       EPA-454/B-26-005. U.S. Environmental Protection Agency, Research Triangle Park,
       North Carolina. July 2026.
U.S. EPA, 2026b: User's Guide for AERSURFACE Tool. EPA-454/B-26-007. U.S.
       Environmental Protection Agency, Research Triangle Park, North Carolina. July 2026.
U.S. EPA, 2026c: User’s Guide for the AMS/EPA Regulatory Model (AERMOD). EPA-454/B-
       26-001. U.S. Environmental Protection Agency, Research Triangle Park, North Carolina.
       July 2026.
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.
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.
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.



                                                  174
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.
Venkatram A., Karamchandani P., Pai P. & Goldstein R. (1994). The Development and
      Application of a Simplified Ozone Modelling System (SOMS). Atmos. Environ. 28 (22),
      pp 3365-3678. doi:10.1016/1352-2310(94)00190-V
Venkatram, A., Snyder, M.G., Heist, D.K., Perry, S.G., Petersen, W.B. and Isakov, V. (2013).
      Re-formulation of plume spread for near-surface dispersion. Atmospheric environment,
      77, pp.846-855.
Venkatram, A., Heist, D.K., Perry, S.G., and Brouwer, L. (2021). Dispersion at the edges of near
      road noise barriers. Atmospheric Pollution Research, 12, pp. 367–374.
Walcek, C., G. Stensland, L. Zhang, H. Huang, J. Hales, C. Sweet, W. Massman, A. Williams, J,
      Dicke, 2001: Scientific Peer Review of the Report "Deposition Parameterization for the
      Industrial Source Complex (ISC3) Model." The KEVRIC Company, Durham, North
      Carolina.
Warren, J. C., Robert J. Paine, Jeffrey A. Connors, Carlos Szembek & Eladio
      Knipping: 2022. Evaluation of a revised AERMOD treatment of plume dispersion in the
      daytime elevated stable layer, Journal of the Air & Waste Management
      Association, 72:9, 1040-1052, DOI: https://doi.org/10.1080/10962247.2022.2094031.
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.
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.




                                                  175
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.
Weil, J. C.: 2020. New dispersion model for highly-buoyant plumes in the convective boundary
        layer. Modeling report to the Western Australia department of environmental
        conservation. January 2020.
Wesely, M.L, P.V. Doskey, and J.D. Shannon, 2002: Deposition Parameterizations for the
      Industrial Source Complex (ISC3) Model. Draft ANL report ANL/ER/TRB01/003,
      DOE/xx-nnnn, Argonne National Laboratory, Argonne, Illinois 60439.
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.
Yang, B., Gu, J., & Zhang, K. M., 2020. Parameterization of the building downwash and
       sidewash effect using a mixture model. Building and Environment, 172, 106694.
Zilitinkevich, S. S., 1972: On the determination of the height of the Ekman boundary layer.
        Bound.Layer Meteor., 3, 141-145.




                                                   176
United States              Office of State Air Partnerships         Publication No. EPA-454/B-26-003
Environmental Protection   Air Quality Assessment Division          July 2026
Agency                     Research Triangle Park, NC




                                                              177
