NOAA Technical Memorandum ERL ARL-224

DESCRIPTION OF THE HYSPLIT_4 MODELING SYSTEM

Roland R. Draxler

Air Resources Laboratory

Silver Spring, Maryland

G.D. Hess

Bureau of Meteorology Research Centre

Melbourne, Australia

Air Resources Laboratory

Silver Spring, Maryland

December 1997

Revised: 1998, August

2002, September

2003, October

2004, January

2004, November

2010, December

2014, September
NOTICE

Mention of a commercial company or product does not constitute 

an endorsement by NOAA Environmental Research Laboratories. 

Use of information from this publication concerning proprietary 

products or the test of such products for publicity or advertising 

purposes is not authorized.
CONTENTS

ABSTRACT 1

1. INTRODUCTION 1

2. METEOROLOGICAL DATA FIELDS 1

2.1 Pressure-Absolute 2

2.2 Pressure-Sigma 4

2.3 Terrain-Sigma 4

2.4 Hybrid Absolute Pressure Sigma 4

2.5 Vertical Motion 5

2.6 Solar Radiation at the Earth’s Surface 5

3. ADVECTION 6

4. DISPERSION CALCULATION 7

4.1 Stability 7

4.2 Estimating Heat and Momentum Fluxes 8

4.3 Vertical Mixing Coefficient 9

4.4 Horizontal Mixing Coefficient 10

4.5 Particle Dispersion 11

4.6 Turbulent Velocity Variances 12

4.7 Puff Dispersion 14

4.8 Air Concentration Calculation 15

5. DEPOSITION 16

5.1 Gravitational Settling 16

5.2 Explicit Dry Deposition 17

5.3 Wet Removal 17

5.4 Radioactive Decay 18

5.5 Dry Deposition via Resistance Method 19

5.6 Pollutant Resuspension 21

6. SUMMARY 22

7. ACKNOWLEDGMENTS 22

8. REFERENCES 22
DESCRIPTION OF THE HYSPLIT_4 MODELING SYSTEM

ABSTRACT. The HYSPLIT_4 (HYbrid Single-Particle Lagrangian Integrated 

Trajectory) model (Version 4) is a complete system for computing simple 

trajectories to complex dispersion and deposition simulations using either puff or 

particle approaches. The model uses previously gridded meteorological data on one 

of three conformal map projections (Polar, Lambert, Mercator). Air concentration 

calculations associate the mass of the pollutant species with the release of either 

puffs, particles, or a combination of both. The dispersion rate is calculated from the 

vertical diffusivity profile, wind shear, and horizontal deformation of the wind field. 

Air concentrations are calculated at a specific grid point for puffs and as cellaverage concentrations for particles.

1. INTRODUCTION

The HYSPLIT model has evolved over several stages during the last decade and a half. 

The initial version of the model (Draxler and Taylor, 1982) used only rawinsonde observations 

and the dispersion was assumed to consist of uniform mixing during the daytime and no mixing 

at night. Dispersion due to wind shear was introduced by splitting up the daytime mixed layer 

into smaller layers each night. In the next revision (Draxler and Stunder, 1988), variable strength 

mixing was introduced based upon a temporally and spatially varying diffusivity profile. In the 

last revision (HYSPLIT_3 - Draxler, 1990; 1992), the use of rawinsonde data was replaced by 

gridded meteorological data from either analyses or short-term forecasts from routine numerical 

weather prediction models.

Several significant new features have been added to the current version of the model that 

will be described in more detail below. In particular, the advection algorithms have been 

updated to include temporal interpolation. Although the transport and dispersion of a pollutant 

can still be calculated by assuming the release of a single puff, the puff can now be defined with 

either a Gaussian or top-hat horizontal distribution. As in previous versions, a single released 

puff will expand until its size exceeds the meteorological grid cell spacing, and then it will split 

into several new smaller puffs. Further, a three-dimensional particle dispersion routine has been 

added that computes air concentrations from the dispersal of an initial fixed number of particles. 

All the equations used to compute the strength of the vertical mixing have been revised based 

upon the more recent literature, and the rate of horizontal dispersion is no longer assumed to be 

constant but can vary according to the deformation of the wind field. Each of these features will 

be discussed in more detail in the following sections.

2. METEOROLOGICAL DATA FIELDS

Meteorological model output fields usually cannot be directly used by the dispersion 

model without some pre-processing, primarily because the data may have been interpolated to a 

variety of different vertical coordinate systems prior to output. To maintain a larger degree of 

flexibility within the dispersion model structure, i.e., the ability to use different meteorological 

data sources for input, the meteorological profiles at each horizontal grid point are linearly 

interpolated to an internal dispersion model terrain-following (σ) coordinate system, 
σ = (Ztop - Zmsl) / (Ztop - Zgl), (1)

where all the heights are expressed relative to mean-sea level and where Ztop is the top of the 

dispersion model’s coordinate system. The model’s internal heights above ground level can be 

chosen at any interval, however a quadratic relationship between height (z) and model level was 

specified, such that each level’s height with respect to the model’s internal index, k, is defined by

z = ak2 + bk +c, (2)

where a=30, b=-25, and c=5. This relation results in decreasing resolution away from the 

surface, with the first level (k=1) at 10 m, the second level at 75 m, the third at 200 m, while by 

the 20th level, at 11500 m, the difference between levels is about 1200 m. Any vertical resolution 

can be specified by altering the constants in Eq. (2), however the model’s internal resolution 

should be at the same or better vertical resolution than the input data.

The dispersion model’s horizontal grid system is designed to be identical to that of the 

meteorological data. Three different conformal map projections are supported, Polar 

Stereographic, Mercator, and Lambert Conformal, using a set of universal mapping 

transformation routines (Taylor, 1997). 

Gridded fields of meteorological variables are also required at regular temporal intervals. 

The time interval between fields should be constant for each defined grid, i.e., the fine grid 

regional data may be available at 3 h intervals while the coarser grid global model fields may be 

available only every 6 h. 

Meteorological data fields may be provided on one of four different vertical coordinate 

systems: pressure-sigma, pressure-absolute, terrain-sigma, or a hybrid absolute-pressure-sigma 

(used by ECMWF - the European Centre for Medium-Range Weather Forecasting). At a 

minimum the model requires U, V (the horizontal wind components), T (temperature), Z (height) 

or P (pressure), and the pressure at the surface, P0. Moisture and vertical motion are optional, 

however the vertical motion may be computed based upon how the vertical coordinate is defined. 

If wet deposition processes for soluble gases or for particles are to be included, the model also 

requires the rainfall field. The data conversions required for each of these inputs are described in 

more detail in the following sections.

2.1 Pressure-Absolute

When the input data are given on pressure surfaces, the only surface value required is 

pressure. Heights of the pressure surfaces (Zp) are assumed to be relative to msl (Mean Sea 

Level). Model output fields on absolute pressure surfaces have usually been interpolated to 

those surfaces. Moisture output is expected as relative humidity (RH) while temperature is 

assumed to be a virtual temperature (Tv). Surface values are computed if not otherwise provided. 

For instance, the surface temperature (T0) is built down adiabatically from the lowest data level 

using the input data level-1 temperature (T1) and pressure (P1),

T0 = T1 (P1/P0)

-0.286. (3)

The height of ground surface (Zg) is interpolated directly from the profile if the ground level falls 
between input data levels. Otherwise the ground surface is estimated from the equation of state, 

Zg = Z1 - Rd T01 ln (P0/P1) g-1

, (4)

where Rd is the gas constant for dry air (287.04 J kg-1 K-1), T01 is the average temperature 

between the two lowest data levels, and g is the acceleration of gravity.

Normally ground-level (≤10m) winds are available from the meteorological model. In the 

rare circumstances when these fields are not provided, the low-level winds (Umodel) are estimated 

at all model levels (Zmodel) below the lowest data level (Zdata) using a logarithmic profile for 

neutral conditions, 

Umodel = Udata ln (Zmodel /Z0) / ln (Zdata /Z0), (5)

where Z0 is the roughness length and the subscripts data and model refer to the values at the 

lowest level in the input data file and the internal model coordinate system, respectively. Neutral 

conditions are assumed because stability data have not yet been processed at this stage of the 

calculations.

Upper-level fields (such as U,V,T,RH) are linearly interpolated from the data level heights 

to the model levels. Several additional dependent variables are also computed at each model 

level such as the local density, 

ρ = k1 P T-1 Rd

-1

, (6)

where P is the local pressure and k1 (= 100) is a conversion factor from hPa to N m-2 , and the 

potential temperature

θ = T (1000/P)0.286. (7)

If some measure of vertical motion (Wp) is provided (usually in units of hPa s-1

), it is converted 

to (Wσ) sigma coordinates (s-1

),

Wσ = k1 Wp ρ-1 g-1 Ztop-1

. (8)

Equation (8) neglects the terrain slope because the output fields archived from most 

meteorological models, although typically output on pressure surfaces, are usually the same 

fields as calculated on the meteorological model’s internal sigma coordinate system, only 

interpolated to pressure surfaces.
2.2 Pressure-Sigma

When the input data are given on pressure-sigma surfaces it is assumed that these surfaces 

are the native grid of the meteorological model, hence moisture is expected as specific humidity 

and temperature is assumed to be a dry temperature. Therefore the virtual temperature is first 

computed from the specific humidity (Q)

Tv = T ( 1 + 0.61 Q). (9)

before conversion to potential temperature. Heights are computed for each data level because 

pressure-sigma data files do not contain heights. The hypsometric equation is integrated from 

the surface to obtain the height of each data level based upon the layer average virtual 

temperature. Therefore the height increment between levels 1 and 2 becomes,

Δz = ln(P1/P2) Rd Tv12 g-1. (10)

The specific humidity is converted to relative humidity fraction (0 to 1),

RH ≈ Q P ( 0.622 Es )

-1

, (11)

from the specific humidity (Q) and from the saturation mixing ratio,

Es = exp ( 21.4 - [5351 T-1

] ), (12)

expressed in hPa. Relative humidity is required in calculations, described in more detail in later 

sections, relating to determining cloud levels for wet removal processes, incident solar radiation 

for dry deposition, and chemical conversions. Surface pressure, density, vertical motion, and 

other variables are computed as described in the previous section.

2.3 Terrain-Sigma

Terrain-following sigma coordinates are typically used by highly detailed mesoscale 

models covering a limited spatial domain. As in the previous sections, surface values are 

computed if needed, however the height of the surface is not required because the input data are 

already terrain following. Upper-level data are processed in a manner similar to the other 

coordinate systems, but the terrain-sigma coordinate system of the input data needs to be 

converted to the terrain-sigma system of the dispersion model following the same procedure as 

with internal model data, except using the meteorological model's top scaling height.

2.4 Hybrid Absolute Pressure Sigma

A hybrid coordinate system, typical of output fields available from the ECMWF (1995), 

consists of an absolute pressure added to the pressure-sigma level. The conversion process is 

almost identical to normal pressure-sigma data as described previously. Input data are defined 

on 

hybrid pressure levels given by the relation

Pk = Ak + Bk P0, (13)
where k indicates the index of the level, P the pressure, P0 the surface pressure, A a pressure 

offset, and B a sigma value (0-1).

2.5 Vertical Motion

In most circumstances the input meteorological data will contain a vertical motion field, 

usually in pressure units, and regardless upon which vertical coordinate system these input data 

are provided, the vertical velocity field is almost always relative to the meteorological model’s 

native terrain-following sigma coordinate system. The trajectory and dispersion model 

calculations can use these data fields directly because the model’s internal coordinate system will 

always be terrain following regardless of the form of the input data. This is one of the primary 

reasons that the input data need to be remapped to a common vertical coordinate system.

When the vertical motion fields are missing, or perhaps there are some special conditions 

required for a simulation, the dispersion model has an option to replace these fields with an 

internally calculated vertical velocity based upon an assumption that the pollutant parcel is 

transported on some other surface. The input data can be remapped to various surfaces by 

computing the velocity (Wη ) required to maintain a parcel on the selected (η) surface,

Wη = (-∂η/∂t -u ∂η/∂x -v ∂η/∂y) / (∂η/∂z), (14)

given the slope of the surface and its local rate of change and where the surfaces, η, can be either 

isobaric (p), isosigma (σ), isopycnic (ρ), or isentropic (θ).

2.6 Solar Radiation at the Earth's Surface

Some special applications require the additional parameter of solar radiation at the earth's 

surface. If it is not available in the input data file it may be computed at each meteorological 

grid point based upon the cloud-cover and sine of the solar elevation angle such that

G = sin(α) S τ, (15)

where S is the solar constant incident at the top of the cloud layer and τ is the fraction transmitted 

through the clouds. Using the normal top-of-the-atmosphere solar constant to be about 1380 W 

m-2 and assuming that on average about 20% of the radiation is absorbed or reflected back into 

space (Sellers, 1972) by a clear atmosphere, then S = 1100 W m-2. The effect of cloud-cover can 

be treated as a linear function of cloud-cover by defining a critical transmissivity (τc) to be the 

limiting factor of 50% transmission at 100% cloud coverage as used by Koch et al. (1985), such 

that

Τ = 1 - A (1-τc), (16)

where A is the fractional cloud-cover. Following a procedure similar to that used in the Nested 

Grid Model (NGM - Tuccillo, 1988), A is determined from the average fractional relative 

humidity (RH - 0 to 1) of the three model sigma layers about the layer with the maximum RH, 

so that at each grid location

A = [ 5 (RH-0.80) ]2

, (17)
and where the RH is limited to a range of 0.80 to 1.0, which yields cloud amounts of 0 to 1. For 

computational simplicity the effects of cloud-depth and height are not considered. However if 

more detailed cloud information and solar fluxes are available from the meteorological data 

archive, the above computations can be replaced directly by the values from the meteorological 

input.

3. ADVECTION

The basis of any Lagrangian model is that the dispersion is computed following the 

particle or puff. That is, the advection of a particle is computed independently. Hence once the 

basic (U,V,W) meteorological data have been processed and interpolated to the internal model 

grid, trajectories (the integrated advection term of a particle) can be computed to test the 

advection components of the model. 

The advection of a particle or puff is computed from the average of the three-dimensional 

velocity vectors for the initial-position P(t) and the first-guess position P'(t+Δt). The velocity 

vectors are linearly interpolated in both space and time. The first guess position is

P'(t+Δt) = P(t) + V(P,t) Δt (18)

and the final position is

P(t+Δt) = P(t) + 0.5 [ V(P,t) + V(P',t+Δt) ] Δt, (19)

The integration method is very common (e.g. Kreyszig, 1968) and has been used for trajectory 

analysis (Petterssen, 1940; Draxler, 1996) for quite some time. Higher order integration methods 

were investigated and rejected because as long as the data observations are linearly interpolated 

from the grid to the integration point, higher order methods will not yield greater precision. 

Trajectories are terminated if they exit the model top, but advection continues along the surface 

if trajectories intersect the ground. The integration time step (Δt) can vary during the simulation. 

It is computed from the requirement that the advection distance per time-step should be less than 

the grid spacing. The maximum transport velocity Umax is determined from the maximum 

particle/puff transport speed during the previous hour. Time steps can vary from 1 minute to 1 

hour and are computed from the relation, 

Umax (grid-units min-1

) Δt (min) < 0.75 (grid-units). (20)

4. DISPERSION CALCULATION

A Lagrangian model can compute air concentrations through either of two assumptions. 

In a puff model, the source is simulated by releasing pollutant puffs at regular intervals over the 

duration of the release. Each puff contains the appropriate fraction of the pollutant mass. The 

puff is advected according to the trajectory of its center position while the size of the puff (both 

horizontally and vertically) expands in time to account for the dispersive nature of a turbulent 

atmosphere. In a Lagrangian particle model, the source can be simulated by releasing many 

particles over the duration of the release. In addition to the advective motion of each particle, a 

random component to the motion is added at each step according to the atmospheric turbulence 

at that time. In this way a cluster of particles released at the same point will expand in space and 

time simulating the dispersive nature of the atmosphere. In a homogeneous environment the size 

of the puff (in terms of its standard deviation) at any particular time should correspond to the 
second moment of the particle positions. An hybrid approach, developed by Hurley (1994), is 

incorporated into this version of the model, in which the calculation uses particle dispersion in 

the vertical direction and puff dispersion in the horizontal. Regardless of which approach is 

used, stability and mixing coefficients need to be computed from the meteorological data. These 

computations, as well as the details of the puff and particle models are discussed in the following 

sections.

4.1 Stability

There are two options to estimate the boundary layer stability. The preferred method is to 

use the fluxes of heat and momentum provided by the meteorological model, if available. 

Otherwise the temperature and wind gradients of each grid-point sounding are used to estimate 

stability. This latter situation may not be ideal if meteorological data aloft are only provided on 

mandatory surfaces with large distances between levels near the ground. In either situation the 

friction velocity and temperature are estimated first, from which the Obukhov length is 

calculated.

First the boundary layer depth (Zi) is computed at each grid point. It is assumed to equal 

the height at which the potential temperature first exceeds the value at the ground by 2K. The 

temperature profile is analyzed from the top down to determine the boundary layer depth. The 

top-down approach reduces the influence of shallow stable layers near the ground. In addition, a 

minimum depth of 250 m is assumed for all hours. The height was chosen to correspond with 

the minimum height resolution typical of the meteorological input data. Night-time depths are 

probably overestimated. 

When surface fluxes are available from the meteorological model, and depending upon the 

variables available in the model output, the friction velocity is computed either from the scalar 

exchange coefficient (E),

u* = (E u / ρ)0.5, (21)

where the product E u is equivalent to the stress (E = ρCD u where CD is the drag coefficient), or 

it is computed directly from the vector momentum fluxes (F),

u* = ( ∣-F∣ / ρ)0.5. (22)

The vector momentum fluxes (N m-2) are converted to a scalar before computation of the friction 

velocity. The friction temperature is always computed from the sensible heat flux (H),

T* = -H (ρ Cp u*)

-1

, (23)

where Cp is the heat constant (1005 J kg-1 K-1 ) for dry air. The convective velocity scale is then 

W* = ∣*T* g u Zi T-1 ∣1/3, (24)

where Zi is the height of the convective boundary layer. At this point the Obukhov length L is 

computed from the friction values to be consistent with model derived flux fields, such that

z/L = Z2 k g T* (u*

2 T2)

-1

, -2 ≤ z/L ≤ 10, (25)
and where Z2 indicates height of 2nd model level, assumed to be the depth of the surface layer, k 

is Von Karman’s constant (0.40), and g is the acceleration of gravity (9.8 m s-2

)

4.2 Estimating Heat and Momentum Fluxes

When no fluxes are provided by the meteorological model, a typical situation especially 

when using certain data archives, z/L is estimated from wind and temperature profiles. The 

meteorological sounding is used first to compute the bulk Richardson Number,

Rb’ = g Δθ ΔZ {θ12 [(Δu)2

+(Δv)2]}-1, (26)

where Δ indicates a gradient between levels 1-2 and θ12 is the layer-average potential 

temperature. The Rb’ is adjusted to account for those meteorological data files that have very 

coarse vertical spacing, where the meteorological values at the top of the surface layer (Z2) have 

been extrapolated from input data levels (Zd) at much greater heights. In those special situations

we assume that

Rb ≈ Rb’ (Z2/Zd)

2

, (27)

to compensate for the increase of Rb’ with increasing layer depth. The bulk Richardson number 

can be shown to be proportional to the square of the geometric mean height between the bottom 

and top of the layer (Golder, 1972). Then the stability parameter z/L, where L is the Obukhov 

length, is estimated from Rb. This is done by using empirical interpolation formulas fit to 

Monin-Obukhov profile functions derived for a surface layer height of 75 m (similar equations 

have been derived by Launiainen (1995) and Abdella and McFarlane (1996)), 

z/L = Rb (s2

/t - 0.50), Rb ≤ 0, (28)

z/L = [-t+2s βRb + (t2 - 4s t βRb + 4s2 βRb )

1⁄2 ]

[2 β (1-βRb )]-1, 0≥Rb<0.08, (29)

z/L = (0.005 s +41.2) Rb

2 + (1.18 s - 1.5 v - 1.37) Rb ,Rb ≥ 0.08, (30)

and where

s = ln (z/Zo + 1.0), (31)

t = ln (z/Zh + 1.0), and (32)

v = ln (Zo/Zh). (33)

Here Zh is the roughness length for heat, Zo is the roughness length for momentum, and 

Z0 / Zh = 10. The parameter β = 5. Equation (29) is an exact solution for a log-linear profile 

(Choudhury, et al. 1986). The friction velocity and temperature are then given by

u* = k Z2 Δu (φm ΔZ)-1

, (34)
T* = k Z2 Δθ (φh ΔZ)-1

, (35)

where k is von Kármán's constant (k ≈ 0.4), and the normalized profiles (φ) for heat (h) and 

momentum (m) are from Beljaars and Holtslag (1991) for a stable surface layer (0≤ z/L≤10), 

φm = 1 + z/L [ a + b exp(-d z/L) (1 - d z/L + c) ], (36)

φh = Prn {1 + z/L [ a ( 1+ a b z/L)1⁄2 + b exp(-d z/L) (1 - d z/L + c) ] }, (37)

where Prn is the turbulent Prandtl number (0.923) during neutral conditions, and a=1, b= 2/3, 

c=5, and d=0.35. In an unstable surface layer (-2 ≤ z/L ≤ 0) Betchov and Yaglom (1971) and 

Kadar and Perepelkin (1989) propose

φm = { [1+0.625 (z/L)2] / (1- 7.5 z/L)}1/3, and (38)

φh = 0.64 { [3 - 2.5 z/L] / [1 - 10 z/L + 50 (z/L)2 ] }1/3. (39)

4.3 Vertical Mixing Coefficient

Pollutant vertical mixing is assumed to follow the coefficients for heat. Within the 

Boundary Layer (BL), vertical mixing coefficients are computed following Troen and Mahrt 

(1986) and Holtslag and Boville (1993), where

Kh = k wh z (1 - z/Zi)

2

, (40)

wh = u* φh

-1

, z/Zi ≤ 0.1, (41)

wh = wm Pr-1, 1 > z/Zi > 0.1. (42)

The Prandtl number, with φh and φm evaluated at z/Zi= 0.1, is given by

Pr = (φh / φm ) + 7.2 k ( z/Zi ) (W*/wm ), where (43)

wm = (u*

3 + 0.6 W*3

)

1/3, (44)

and W*=0 for neutral and stable conditions. 

Once the K profile has been computed, a single average value for the entire BL is 

computed from the profile and that value replaces all the values within the BL. Each horizontal 

grid point will have a different value.

Following Beljaars and Betts (1993), mixing through the inversion layer at z = Zi during 

convective conditions (W*>0) is computed based upon the surface flux parameters and the 

strength of the inversion, where

Kh = Cs H ( ρ Cp ∂θ/∂z )-1, or (45)

Kh = - Cs u* T* (∂θ/∂z)-1, (46)
and Cs = 0.4. During stable conditions the mixing is calculated as described below.

Above the BL, pollutant mixing in the remainder of the atmosphere is defined by the 

vertical diffusivity for heat using mixing length theory where

Kh = l

2 ∣ ∂V /∂z h∣ φ(l/Lo)-1

, (47)

where l is a Blackadar-type mixing length in meters,

l

-1 = kz-1 + 150-1

, (48)

and Lo is the local Obukhov length. The stability function φh is given by Eq. (37), but with z/L 

equal to the mixing length ratio 

l/Lo = 1.0893 Rib (49)

during near-neutral conditions (0 ≤ Rib ≤ 0.001) and otherwise ( 0.001 ≤ Rib ≤ 20 )

l/Lo = a1 + Rib a2 + Rib [a3 + Rib (a4 + a5 Rib )] , (50)

and where a1=0.2828x10-3

, a2=0.8049, a3=1.6583, a4=0.5090x10-2, and a5=-1.0063x10-3.

4.4 Horizontal Mixing Coefficient

The subgrid-scale horizontal mixing coefficient is computed from the velocity 

deformation (Smagorinsky, 1963; Deardorff, 1973),

Khor = 2-0.5 (c Χ)2 [ (∂v/∂x + ∂u/∂y)2 + (∂u/∂x - ∂v/∂y)2 ]

0.5, (51) 

where Χ is the meteorological data grid size, and c is 0.14. Grid scale dispersion is simulated by 

the particles and puffs horizontally spreading into different wind regimes.

4.5 Particle Dispersion

Both the puff and particle dispersion equations are formulated in terms of the turbulent 

velocity components. These velocity components are a function of the turbulent diffusivities 

computed in the previous section. In the particle implementation of the model, the dispersion 

process is represented by adding a turbulent component to the mean velocity obtained from the 

meteorological data. The particle model can be applied in either the vertical, horizontal, or both 

directions. The specific approach used follows the one described by Fay et al. (1995).

After computation of the new position at a time step due to the mean advection of the 

wind, a turbulent component is added to the mean particle positions (X,Z), 

Xfinal(t+Δt) = Xmean(t+Δt) + U'(t+Δt) Δt G, (52)

Zfinal(t+Δt) = Zmean(t+Δt) + W'(t+Δt) Δt Ztop-1

, (53)

where the horizontal and vertical positions are given in grid and sigma units, respectively, while 
the turbulent velocity components are in m s-1

. G and Ztop are the required unit conversion 

factors. The contribution of the turbulent wind components (U' - horizontal, W' - vertical) are 

added to the "mean" position (due only to the mean flow) to give a final position from which the 

advection at the next time step is computed. Full reflection is assumed for particles that intersect 

the ground or model-top. The integration time step is computed from the requirement that the 

change in the vertical plume dimension,

Δzp < 0.5 Δz, (54)

where Δz is the vertical grid spacing, and hence

Δt = (Δz)2 / ( 8 σw

2 TLw ). (55)

The vertical velocity variance, σw

2

, and the Lagrangian time scale, TLw, are discussed in more 

detail below.

The horizontal turbulent velocity components at the current time U'(t+Δt) are computed 

from the turbulent velocity components at the previous time U'(t), an auto-correlation coefficient 

(R) that depends upon the time step, the Lagrangian time scale, and a computer generated 

random component ("); therefore

U'(t+Δt) = R(Δt) U'(t) + U" ( 1-R(Δt)2 )

0.5, (56)

and for the vertical turbulent velocity,

W'/σw (t+Δt) = R(Δt) (W'(t)/σw(t) ) + (W"/σw(t) ) ( 1-R(Δt)2 )

0.5

+ TLw (1-R(Δt)) ∂σw(t)/∂z , (57)

and as defined by Wilson et al. (1983),

σw (t+Δt) = σw(t) + W'(t) Δt ∂σw(t)/∂z.

The velocity variance gradient term on the vertical turbulent velocity is applied to prevent 

accumulation of particles in low turbulence regions (Legg and Raupach, 1982). The importance 

of this latter term is somewhat reduced due to the averaging of the diffusivity profile within the 

BL. In computing the turbulence an exponential velocity autocorrelation is assumed, 

R(Δt) = exp ( -Δt / TLi ), (58)

such that Tli= [TLw or TLu ], and where Tlw (vertical)=100 s and Tlu (horizontal)=10800 s ( ∼1/f) 

are assumed to be constant for convenience. These values result in a random walk (R≈0) vertical 

dispersion for most of the longer time steps. The Gaussian random component U" or W" comes 

from the computer generated random number, such that for any directional component, 

U" (or W") = σi λ, (59)

where λ is a Gaussian random number with mean of 0 and a standard deviation of 1. The 

standard deviation of the turbulent velocities is estimated from the previously calculated vertical 

(Kh - Eqs. 40, 46, 47) and horizontal diffusivities (Khor - Eq. 51),
σi = (Ki / TLi )

0.5, (60)

where i represents the appropriate directional component (u or w). 

4.6 Turbulent Velocity Variances

The turbulent velocity variance can also be obtained directly from the stability functions 

instead of through the intermediate step of computing a diffusion coefficient (Eq. 60). This 

permits the use of the meteorological model’s TKE (turbulent kinetic energy field) if it is 

available. The boundary layer velocity variances are defined as a function of u*, w*, and Zi

following Kantha and Clayson (2000). This method does not use the diffusivity and hence no 

assumptions are required about turbulent scales. Either the turbulent velocity or the diffusivity 

approach can be selected for a simulation. In the stable or neutral boundary layer, 

w’2 = 3.0 u*

2 (1 – z/zi)

3/2, (61)

u’2 = 4.0 u*

2 (1 – z/zi)

3/2, (62)

v’2 = 4.5 u*

2 (1 – z/zi)

3/2. (63)

In the stable or neutral surface layer,

w’2 = 3.0 u*

2

, (64)

u’2 = 4.0 u*

2

, (65)

v’2 = 4.5 u*

2

. (66)

If the TKE field is available, then the velocity variances can be computed from its definition,

E = 0.5 (u’2 + v’2 + w’2

), (67)

and equations 61-66, resulting in the following constant relationships to the TKE:

w’2 = 0.52 E, (68)

u’2 = 0.70 E, (69)

v’2 = 0.78 E. (70)

In the unstable boundary layer, again following Kantha and Clayson (2000), 

w’2= w*

2 (z/zi)

2/3 (1 – z/zi)2/3 (1 + 0.5 R2/3), (71)

where R is the ration of the heat flux at the inversion to the flux at the surface. It is assumed to be 

constant at about 0.2. From Garratt (1992), the horizontal variances are simply a constant 

function of the convective velocity scale,

u’2 = v’2 = 0.36 w*

2

. (72)

In the unstable surface layer, 

w’2 = 1.74 u*

2 (1 – 3 z/L)2/3, (73)
u’2 = v’2 = 0.36 w*

2

. (74)

If the TKE field is available, then from the definition (67) and Eqs. 71-74,

w’2 = 2 E / (1+Kx), (75)

u’2 = v’2

 = 0.5 Kx w’2, (76)

and where in the unstable surface layer Kx is defined to be 

Kx = 0.41 (w*

2

/ u*

2 ) / (1 – 3 z/L)2/3, (77)

and in the unstable boundary layer,

Kx = 0.72 / [ (z/zi)

2/3 (1 – z/zi)2/3 (1 + 0.5 R 2/3) ]. (78)

Unlike the diffusivity (Eq. 40) approach, the vertical velocity variance is not averaged for 

a single boundary layer value. The turbulence varies by height according to Eqs. 61-76. The unaveraged turbulent velocities are much more realistic for short-range dispersion simulations. 

Note that the horizontal turbulent velocity also replaces the deformation based diffusvity (Eq. 

51). If the TKE field is available, the mixed layer depth can also be computed from the TKE 

profile instead of the temperature profile. In this case the mixed layer depth is assumed to be the 

height above ground at which the TKE drops by more than half of its previous value. 

4.7 Puff Dispersion

Puff dispersion is treated in two domains, when the puff is smaller than the meteorological 

model grid size and when it is larger. In the latter case it is assumed that the meteorological 

model is capable of resolving turbulent motions on that scale. Gaussian and "top-hat" puffs are 

treated almost identically.

When puffs have dimensions less than the meteorological grid spacing, the rate of vertical 

puff growth is assumed to be, 

dσz

2

/dt = 2 σw

2 TL, (79)

and there are two options for the horizontal growth rate,

dσh/dt = σu , (Growth linear with time) (80)

dσh/dt = σu (0.5 TL / t)0.5 (Growth proportional to square-root of time) (80a)

Turbulent velocity variances are computed as discussed in the particle model section. Additional 

variations to the growth rates are achieved through the puff-splitting process as the puff grows 

into different mixing regions.

In the vertical direction, a puff distribution is always assumed to be “top-hat”, that is, a 

constant value inside the puff and zero outside. The edge of a top-hat puff is assumed to be at an 

abscissa value of 1.54σz, the point on a Gaussian distribution where the areas above and below 

the corresponding ordinate are equal. Vertical growth rates are computed at the top (t) and 

bottom (b) of the puff (±1.54σ) such that,
σzt

2 (t+Δt) = σz2

(t) + 2 σwt

2

(t) TL Δt, and (81)

σzb

2 (t+Δt) = σz2

(t) + 2 σwb

2

(t) TL Δt. (82)

The final vertical puff standard deviation is then just the average of the two,

σz(t+Δt) = 0.5 [σzt(t+Δt) + σzb(t+Δt)]. (83)

At each time step, for puffs near the model domain limits, the σz is truncated so that a puff 

cannot grow below the ground surface nor the model-top. The puff horizontal standard deviation

σh(t+Δt) = σh(t) + 2 σu(t) Δt, (84)

is evaluated from the turbulent velocity component at the puff center position.

When a puff expands to cover several meteorological grid points, a top-hat puff splits 

horizontally into four puffs, each with 25% of the mass, when 1.54σh > Lh, at positions 

P(x±0.5σh, y±0.5σh), and where Lh is the meteorological model grid size. A large Gaussian puff 

splits into five smaller puffs when 3.0σh > Lh, at the same positions as the top-hat and with an 

additional puff at the center position. The center puff gets 60% of the mass while the outside 4 

puffs get 10% each. Puffs split vertically into “n” components when ±1.54σz > 2ΔZ, where ΔZ 

is the vertical grid size and “n” is the number of ΔZ layers within ±1.54 σz . Each new puff has a 

position, Pn (z), at the center height of “n” layers, each of depth ±1.54 σz /n.

Puff splitting can quickly exceed array dimension space. There are three mechanisms to 

remove excessive puffs. (1) Every hour, puffs are spatially sorted so that puffs near to each other 

are in contiguous array locations. Puffs whose centers are within 1.0 σh and 0.5 σz , and their σh

are within 0.1 of each other, are merged together. The new puff's dispersion coefficients are a 

mass-weighted sum of the contributing individual puffs. (2) Every six hours all the puffs are 

sorted by mass, and those puffs whose accumulated mass is less than 10% of the total mass are 

again sorted by position and merged with less restrictive criteria; i.e., centers are within 1.75 σh 

and 2.0 σz, and σh within 0.20. (3) In addition, limits can be set as to the maximum age of a 

pollutant and the minimum mass that any one puff is permitted to retain. These limits may be 

modified according to the problem under consideration.

4.7 Air Concentration Calculation 

Puff distributions may be defined in either the vertical and horizontal directions, or only in 

the horizontal direction. For each puff, concentrations are summed at each time step to all grid 

points that fall within the puff extent defined for top-hat distributions (±1.54 σi), where i 

indicates z or h, or Gaussian distributions (±3.0 σh). Vertical distributions are always defined as 

top-hat while horizontal distributions may be either. The incremental concentration contribution 

by each puff of mass m to a grid point is computed as follows for a top-hat puff, 

Δc = m (π r2 Δz)-1, (85)

where the vertical extent Δz = 3.08 σz, and the horizontal radius r = 1.54 σh. All grid-nodes 

within the puff extent receive the same Δc. The incremental concentration contribution for a 
Gaussian puff is, 

Δc = m (2 π σh

2 Δz)-1

 exp(-0.5 x2/σh

2

), (86)

where x is the distance from the puff center to grid-node, and the other terms are as previously 

defined.

Particle calculations can be performed in either the vertical or both the vertical and 

horizontal directions. Calculations with a vertical particle distribution may have either a top-hat 

or Gaussian puff horizontal distribution. However, particle calculations are summed into a gridcell rather than computed at a grid-point. A cell is defined at the center of the node and has an 

area corresponding to the half-way distance to adjacent nodes. The incremental concentration 

contribution to a cell by a single particle of mass m is defined for a 3D particle,

Δc = m (Δx Δy Δz)-1, (87)

where Δx,Δy,Δz are the grid-cell dimensions. For a particle with a horizontal top-hat the 

incremental concentration is the same as Eq. (85), but with Δz defined as grid-cell height. If the 

horizontal distribution is Gaussian then the incremental concentration is the same as Eq. (86), but 

again with Δz defined as the grid-cell height. The incremental concentrations are added to each 

grid cell or node each advection time step for all particles or puffs that intersect that point. The 

final average concentration is the incremental sum divided by the number of time steps in the 

concentration averaging period. To avoid the situation where particles or puffs might skip a grid 

point due to large advection time steps, the time step computed in Eq. (20) uses the concentration 

grid spacing instead of the meteorological grid spacing when the model is used to compute air 

concentrations rather than just trajectories.

5. DEPOSITION

There are three different removal mechanisms available: dry deposition, wet depletion, 

and radioactive decay. Dry deposition is either explicitly defined as a deposition velocity, or for 

particles it may be computed as being the equivalent to the gravitational settling velocity, or it 

may be computed using the resistance method and information about the nature of the surface. 

Computation of particle settling velocity requires the particle diameter and density. Wet removal 

can be defined for soluble gases by specifying its Henry's Law constant. Gaseous wet removal 

only occurs for the fraction of the pollutant below the cloud top. Particle wet removal is defined 

by a scavenging ratio (l/l) within the cloud and by an explicit scavenging coefficient (s-1

) for 

pollutants below the cloud base.

For computational simplicity, the total deposition from both dry and wet removal 

processes is expressed in terms of (reciprocal) time constants. The time constants can be added 

and hence the total deposition over a time step becomes 

Dwet+dry = m {1-exp[-Δt (βdry+βgas+βinc+βbel ) ] }, (88)

where m is the pollutant mass of either the particle or puff. The pollutant mass is then reduced 

by the deposition amount. Each of these time constants, for dry deposition (βdry), wet removal 

for gases (βgas), in-cloud wet removal of particles (βinc), and below-cloud wet removal of 
particles (βbel), will be discussed in more detail in the following sections.

In addition to the mass removal option for dry deposition, a probability based approach is 

also available. In this method any particle that is within the deposition layer may stick to the 

surface and loose all of its mass to deposition. The particle is assigned a random number (0-1) 

every time step, and if that number is less than Δt βdry, the particle will be deposited.

5.1 Gravitational Settling

Particle settling is computed following Van der Hoven (1968), where the settling velocity 

(Vg) is calculated for a spherical particle from the particle diameter (dp), air density (ρ), and 

particle density (ρg),

Vg = dp

2 g (ρg - ρ) (18 μ)-1

, (89)

where μ is the dynamic viscosity of air (0.01789 g-1 

s

-1

). The settling velocity is then adjusted for 

a slip correction (Cc) for small particles and a dynamic shape factor (α) to account for nonspherical particles. Hence the final settling velocity is

Vs = Vg Cc α-1, (90)

where α can vary between 1.0 to 2.0 and the Cunningham slip correction is given by

Cc = 1 + 2 (Λ/d) { 1.26 + 0.4 exp [-1.1 d / (2 Λ) ]}. (91)

The molecular mean free path Λ at ambient conditions is approximated from the value Λstp at 

STP (6.53 x 10-8 m) through the relation,

Λ = Λ stp (ρstp / ρ). (92)

The settling velocity is then applied to the pollutant's vertical position each time step to permit 

the gradual sinking of particles and (non-gaseous) puffs. 

5.2 Explicit Dry Deposition

Dry removal is computed when the bottom of the puff or the particle center position is 

within the surface layer, usually defined internally to the model as the second meteorological 

data level. Then the mass deposited by dry removal 

Ddry = Vd C (93)

can be calculated by assuming a uniform vertical concentration distribution in the deposition 

layer. This is implicit in the puff calculation, because a top-hat distribution is assumed in the 

vertical, and the depth of the pollutant layer, ΔZp, equals ±1.54 σz. In the particle calculation, 

ΔZp defaults to the depth of the surface layer. The deposition velocity is converted to a time 

constant of the form

βdry = Vd ΔZp

-1

, (94)
where the deposition velocity may be directly specified in the input or calculated as the settling 

velocity. Dry deposition velocity may also be computed through the resistance method which is 

discussed in a later section.

An alternative to removing a fraction of the particle’s mass is to compute the probability 

that a particle will deposit all of its mass and then only a fraction of the particles in the 

deposition layer will be deposited each time step. In this situation, if R is a random number from 

0 to 1, then a particle will deposit if R < βdry Δt.

5.3 Wet Removal

Wet deposition (Hicks, 1986) is divided into two processes, those in which the polluted air 

is continuously ingested into a cloud from a polluted boundary layer and those in which rain falls 

through a polluted layer. For particulate pollutants, the simplifying assumption of a scavenging 

ratio is assumed for pollutants located within a cloud layer and the scavenging coefficient is used 

for pollutant removal in rain below a cloud layer. At the grid points where it is raining, the cloud 

bottom is defined at the level when the RH first reaches 80% and the cloud top is reached when 

the RH drops below 60%. All removal amounts are adjusted by the fraction of the pollutant mass 

that is within the cloud layer by defining the fraction of the pollutant layer that is below the cloud 

top (Ft) and the fraction of the pollutant layer that is above the cloud bottom (Fb).

For the wet removal of particles by within-cloud processes a scavenging ratio, which is the 

ratio of the pollutant's concentration in water to its concentration in air, is expressed as a wet 

deposition velocity,

Vinc = S P, (95)

where the precipitation rate is given by P. The time constant for within-cloud removal, 

βinc = Ft Fb Vinc ΔZp

-1

, (96)

where the average scavenging ratio is S=5.0x105 by volume, and as before, ΔZp is the depth of 

the pollutant layer for puffs and the depth of the cloud layer for particles. Different scavenging 

ratios can be defined for different pollutants. Another option for within-cloud removal is to use 

the rate constant directly, similar to the below-cloud procedure, but only when the removal 

coefficient is defined to be less than one, 

βinc = Ft

 Fb 8x10-5 P 0.79. (96a)

Below-cloud removal is defined directly as a rate constant, independent of the precipitation rate

when using Eq. 96 or a function of P (mm/h) when using Eq. 96a for within-cloud removal. The 

below-cloud removal constant (s-1

) is given by,

βbel = 8x10-5 (1.0-Fb), or (97)

βbel = 8x10-5 P 0.79 (1.0-Fb). (97a)

The wet deposition of gases depends upon their solubility and for inert non-reactive gases 

it is a function of the Henry's Law constant (Molar atm-1 ), the ratio of the pollutant's equilibrium 
concentration in water to that in air. A gaseous wet deposition velocity can be defined as

Vgas = H R T P, (98)

where R is the universal gas constant (0.082 atm M-1 K-1), T is temperature, and hence the 

gaseous wet removal time constant, 

βgas = FtVgas ΔZp

-1

. (99)

Note that the wet removal of gases is applied at all levels from the ground to the top of the cloudlayer.

5.4 Radioactive Decay

Although radioactive decay, by itself, does not result in deposition, deposited radioactive 

pollutants do decay, and hence deposition amounts are adjusted for radioactive decay each time 

step. The decay constant for radioactive processes is defined by the half-life T1⁄2,

βrad = ln 2 / T1⁄2 (100)

and the radioactive decay of the pollutant's mass, either in the air or that has been deposited 

becomes

m2 = m1 exp(-βrad Δt). (101)

5.5 Dry Deposition via Resistance Method

Rather than explicitly defining a dry deposition velocity for a pollutant, the total 

deposition velocity can alternately be computed from the sum of various resistances (Hicks, 

1986) and the settling velocity for particles such that the total deposition velocity

Vd = [ Ra + Rb + Rc + Ra Rb Vg ]

-1 + Vg, (102)

where the subscripts for the resistances R, represent the atmospheric layer (a), the quasi-laminar 

sublayer (b), and the canopy layer (c) which represents the bulk resistance of various surfaces. 

Gravitational settling, Vg, is zero for gases and Rc is zero for particles. The resistance 

components depend upon meteorological conditions as well as the properties of the surface. The 

surface properties used in these computations are obtained from external sources. Each of the 

resistance components will be discussed in more detail in the following sections; however the 

method very closely follows that proposed by Wesely (1989) as incorporated into the RADM 

model (Chang, 1990), and updated by Walmsley and Wesely (1996).

Atmospheric resistance parameterizes the limiting role of atmospheric turbulence and is 

incorporated for gases and particles through the aerodynamic resistance (Wesely and Hicks, 

1977) 

Ra = Prn (ln z/Zo - ψh) / k u*. (103)

Over land the aerodynamic roughness length Zo is constant and determined from the land-use 
and vegetative cover; over water Charnock's (1958) relation, as modified by Smith (1988), 

Zo = 0.011 u*

2

/g +υ/(9.1u*), (104)

is used to define the aerodynamic roughness length, where υ is the kinematic viscosity of air 

(υ=μ/ρ). The second term is added to account for light wind cases (u*>0). The stability 

correction for heat is calculated from 

ψh = ∫[(1 - φh/Pr)/(z/L)]d(z/L) (105)

where φh was given by Eqs. (37) and (39), and where z is evaluated at the top of the surface 

layer. For unstable conditions 

ψh = - 2.7283 z/L, -0.001 ≤ z/L ≤ 0 (106)

ψh = a1 + z/L (a2 + z/L (a3 + z/L (a4 + a5 z/L))), -2 ≤ z/L ≤ -0.001 (107)

where a1 = 0.1164x10-4, a2 = -2.7188, a3 = -2.1551, a4 = -0.9859, and a5 = -0.1765. For stable 

conditions

ψh = - (1 + ab z/L)3/2 - b (z/L - c/d) exp(-dz/L) - bc/d + 1, 0 ≤ z/L ≤10 (108)

where a = 1, b = 2/3, c = 5, and d = 0.35.

The quasi-laminar sublayer resistance incorporates the effects of the laminar layer just 

above the surface. Over water the sub-layer resistance is assumed to be small and only limited 

by the atmospheric resistance (Slinn and Slinn, 1980); however over land the resistance for 

gaseous deposition is parameterized through the Schmidt (Sc) number following Wesely and 

Hicks (1977)

Rb = Pr (d1/ku*) Scd2. (109)

In the above relations the constants d1=2, d2=2/3, and with the other constants as defined 

previously. The resistance for particulates is computed from the same relationship with an 

additional impaction term (Raupach, 1993)

Rb = Pr { [ (d1 / ku*) Sc

d2 ]

-1 + u*[ St / (St+p) ]q }-1

, (110)

where the constants p=0.8 and q=2.0 (Peters and Eiden, 1992), and the Stokes number St is 

computed from 

St = (2 Vg u*) (l* g Cc)

-1

, (111)

where l* is the laminar layer length scale (l*= υ/u*) and the Schmidt number is given by

Sc = υ / D. (112)

The diffusivity of a specific gaseous pollutant, D, is related to the ratio of the molecular 

weights of the pollutant (Wp) and air (Wa) through Graham's Law,
D = υ (Wa/Wp)

1⁄2. (113)

For particulate pollutants, the diffusion rate (Seinfeld, 1986, p. 324) is given by

D = kb T Cc ( 3 π μ dp)

-1

, (114)

where kb equals Boltzmann's constant (1.38x10-20 g

2 K-1 

s

-2

) and the other symbols have been 

previously defined.

The canopy resistance depends primarily upon a number of plant physiological and 

ground surface characteristics that control the uptake of gases into plants and act in parallel. As 

noted earlier, Rc is zero for particles. The following procedure follows the equations and 

notation as outlined by Wesely (1989) for the total canopy resistance,

Rc = [ 1/(Rs+Rm) + 1/Rlu + 1/(Rdc+Rcl) + 1/(Rac+Rgs) ]-1

, (115)

and is comprised of the stomatal (s), mesophyll (m), upper canopy leaf cuticles (lu), gas-phase 

transfer by convection (dc), surfaces within the lower canopy (cl), canopy height and density 

factor (ac), and ground surface (gs) resistances. The total canopy resistance is limited to a 

minimum value of 10 s m-1 to prevent unrealistic high deposition velocities under certain 

conditions or terrain.

The stomatal resistance primarily depends upon the solar radiation and pollutant species 

and can be expressed as,

Rs = Ri Dhx [1 + {200/(G+0.1)}2

] [400 / (Ts (40-Ts))], (116)

where G is the solar irradiation in W/m2

, Dhx is the ratio of the diffusivity of water vapor to that 

of the pollutant (Table 2 in Wesely, 1989), Ts is the ambient temperature in degrees Celsius, and 

Ri is the minimum resistance for water vapor (Table 1 in Wesely, 1989), which depends upon 

season and land-cover. For temperatures outside the limits of 0 to 40, respectively, Rs is set to a 

very large value. The incident solar radiation is computed at each meteorological grid point 

based upon the cloud-cover and sine of the solar elevation angle, as described in section 2.6. The 

other resistances depend primarily upon the solubility and reactivity of the pollutant and a simple 

fractional expression can be used for each:

Rm = [ H*

/3000 + 100 fo ]

-1

, (117)

Rlu = Rlu1 [ H* 10-5 + fo ]

-1

, (118)

Rdc = 100 [ 1+1000/(G+10) ], (119)

Rcl = [ H* 10-5

/Rcl1 + fo/Rcl2 ]

-1

, (120)

Rgs = [ H* 10-5

/Rgs1 + fo/Rgs2 ]

-1

, (121)

where H* is an effective (relative to SO2) Henry's constant and fo is a pollutant specific reactivity 

parameter. All constants (including Rac in Eq. 115) are defined in tables given by Wesely (1989) 

and Chang (1990), and they will not be repeated here. These scaling constants depend upon 
land-use and season. The relationship for Rdc has been simplified by assuming that the slope of 

the terrain is zero and the relation for Rlu applies only to dry surfaces. In addition the term (sm1

)

1000 exp(-Ts -4), (122)

is added to Rlu, Rcl, and Rgs, to account for cold-temperature resistance increases. 

5.6 Pollutant Resuspension

Under certain conditions pollutants that have deposited can be resuspended into the 

atmosphere if the winds are sufficiently strong and the material is not bound to the surface. 

Pollutant resuspension is parameterized (INSRP, 1993) using a pollutant resuspension factor that 

represents the ratio of the pollutant concentration in air C, to the amount deposited on the surface 

S, such that

K = C / S, (123)

where K has units of m-1 with typical values on the order 10-6 . The resuspension factor can 

also be expressed as a flux,

K = (R/S) dS/dt, (124)

where R is related to the atmospheric resistance. If we assume that R = (k u*)

-1

, then the upward 

directed resuspension flux is simply

dS/dt = k u* K S. (125)

This process is applied in the model over land, when the deposition process is turned on. At each 

time step, for those concentration grid points which contain a non-zero deposition values, a puff 

or particle with mass computed according to Eq. (125) is emitted. The deposition total at that 

cell is reduced accordingly.

6. SUMMARY

A detailed description and equations of HYSPLIT_4, a Lagrangian model that can be used 

to calculate trajectories and air concentrations has been presented. The configuration of the 

model is very generic, in that it could be set up to perform a variety of different scenarios. In 

general Lagrangian models are well suited for quick calculations from pollutant point sources 

and such a modeling approach is ideal for situations where quick turnaround is essential; i.e., 

about 1 to 3 min CPU time per simulation day on a Pentium workstation.

7. ACKNOWLEDGMENTS

The authors would like to thank Gang Liang, BMRC, for development of the HYSPLIT 

graphical user interface; and Les Logan, BMRC, for developing some of the NCAR graphics 

display routines and an interface for the BoM meteorological data.

8. REFERENCES
Abdella, K., and N.A. McFarlane, 1996. Parameterization of the surface-layer exchange 

coefficients for atmospheric models. Boundary-Layer Meteorol., 80, 233-248.

Beljaars, A.C.M., and A.K. Betts, 1993. Validation of the boundary layer representation in the 

ECMWF model. Proceedings, ECMWF Seminar on Validation of Models over Europe, 

Vol. II, 7-11 September, 1992, European Centre for Medium-Range Weather Forecasts, 

Shinfield Park, UK, 159-195.

Beljaars, A.C.M., and A.A.M. Holtslag, 1991. Flux parameterizations over land surfaces for 

atmospheric models. J. Appl. Meteorol., 30, 327-341.

Betchov, R., and A.M. Yaglom, 1971. Comments on the theory of similarity as applied to 

turbulence in an unstably stratified fluid. Atmos. Oceanic Phys,, 7, 829-832.

Chang, J.S., 1990. Appendix E, NAPAP Report 4, The regional acid deposition model and 

engineering model. In Acid Deposition: State of Science and Technology, Vol I, 

Emissions, Atmospheric Processes, and Deposition, U.S. Government Printing Office, 

Washington D.C., 20402-9325.

Charnock, H., 1958. A note on empirical wind-wave formulae. Q.J.R. Met. Soc., 84, 443 - 447.

Choudhury, B. J., R.J. Reginato, and S.B. Idso, 1986. An analysis of infrared temperature 

observations over wheat and calculation of latent heat flux. Agric. Forest Meteorol., 37, 

75-88.

Deardorff, J.W., 1973. Three-dimensional numerical modeling of the planetary boundary layer. 

In Workshop on Micrometeorology, D.A. Haugen (Ed.), American Meteorological 

Society, Boston, 271-311

Draxler. R.R., 1990. The calculation of low-level winds from the archived data of a 

regional primitive equation model. J. Appl. Meteorol., 29, 240-248.

Draxler, R.R., 1992. Hybrid single-particle Lagrangian integrated trajectories (HY-SPLIT): 

Version 3.0 -- User’s guide and model description. NOAA Tech. Memo. ERL ARL-195, 

26 pp. and Appendices. [Available from National Technical Information Service, 5285 

Port Royal Road, Springfield, VA 22161.]

Draxler, R.R., 1996. Trajectory optimization for balloon flight planning. Weather and 

Forecasting, 11, 111-114.

Draxler, R.R., and B.J.B. Stunder, 1988. Modeling the CAPTEX vertical tracer concentration 

profiles. J. Appl. Meteorol., , 27, 617-625.

Draxler, R.R., and A.D. Taylor, 1982. Horizontal dispersion parameters for long-range transport 

modeling. J. Appl. Meteorol., 21, 367-372.

ECMWF, 1995. The description of the ECMWF/WCRP level III-A global atmospheric data 

archive. European Centre for Medium-Range Weather Forecasts, Reading, Berkshire, 

England.
Fay, B., H. Glaab, I. Jacobsen, and R. Schrodin, 1995. Evaluation of Eulerian and Lagrangian 

atmospheric transport models at the Deutscher Wetterdienst using ANATEX surface tracer 

data. Atmos. Environ., 29, 2485-2497.

Garratt, J.R., 1992. The Atmospheric Boundary Layer, Cambridge University Press, Cambridge, 

United Kingdom, 316 pp.

Golder, D., 1972. Relations among stability parameters in the surface layer. Boundary- Layer 

Meteorol., 3, 47-58.

Hicks, B.B., 1986. Differences in wet and dry particle deposition parameters between North 

America and Europe. In Aerosols: Research, Risk Assessment, and Control Strategies, 

Lewis Publishers, Chelsea, MI, 973-982.

Holtslag, A.A.M., and B.A. Boville, 1993. Local versus nonlocal boundary-layer diffusion in a 

global climate model. J. Climate, 6, 1825-1842.

Hurley, P., 1994. PARTPUFF - A Lagrangian particle-puff approach for plume dispersion 

modeling applications. J. Appl. Meteorol., 33, 285-294.

INSRP, 1993. The role of resuspension of radioactive particles in nuclear assessments. 

Proceedings, Interagency Nuclear Safety Review Panel (INSRP), Technical Interchange 

Meeting, Cocoa Beach, FL, 21-23 September, 32 pp., NOAA Air Resources Laboratory, 

Silver Spring, MD 20910.

Kadar, B.A., and V.G. Perepelkin, 1989. Effect of the unstable stratification on wind and 

temperature profiles in the surface layer. Atmos. Oceanic Phys., 25, 583-588.

Kantha, L.H. and C.A. Clayson, 2000. Small Scale Processes in Geophysical Fluid Flows, Vol. 

67, International Geophysics Series, Academic Press, San Diego, CA, 883 pp.

Koch, S.E., W.C. Skillman, P.J. Kocin, P.J. Wetzel, K.F. Brill, D.A. Keyser, and M.C. 

McCumber, 1985. Synoptic scale forecast skill and systematic errors in the MASS 2.0 

model. Mon. Weather Rev., 113, 1714-1737.

Kreyszig, E., 1968. Advanced Engineering Mathematics. 2nd Ed., J. Wiley and Sons, New 

York, 898 pp.

Launiainen, J., 1995. Derivation of the relationship between the Obukhov stability parameter 

and the bulk Richardson number for flux-profile studies. Boundary- Layer Meteorol., 76, 

165-179.

Legg, B.J., and M. Raupach, 1982. Markov chain simulation of particle dispersion in 

inhomogeneous flows: the mean drift velocity induced by a gradient in Eulerian velocity 

variance. Boundary- Layer Meteorol., 24, 3-13.

Peters, K., and R. Eiden, 1992. Modelling the dry deposition velocity of aerosol particles to a 

spruce forest. Atmos. Environ., 26A, 2555-2564.
Petterssen, S., 1940. Weather Analysis and Forecasting. McGraw-Hill Book Company, New 

York, 221-223.

Raupach, M.R., 1993. Dry deposition of gases and particles to vegetation. Clean Air, 27, 200-

203.

Seinfeld, J.H., 1986. Atmospheric Chemistry and Physics of Air Pollution. John Wiley and 

Sons, New York, 324 pp.

Sellers, W.D., 1972. Physical Climatology. The University of Chicago Press, Chicago, IL 

60637, 272 pp.

Slinn, S.A., and W.G.N. Slinn, 1980. Predictions for particle deposition on natural waters. 

Atmos. Environ., 14, 1013-1016.

Smagorinsky, J., 1963. General circulation experiments with the primitive equations: 1. The 

basic experiment. Mon. Weather Rev., 91, 99-164.

Smith, S.D., 1988. Coefficients for sea surface wind stress, heat flux and wind profiles as a 

function of wind speed and temperature. J. Geophys. Res., 93, C12, 15467-15472.

Taylor, A.D., 1997. Conformal map transformations for meteorological modelers. Computers 

and Geosciences, 23, 63-75.

Troen, I., and L. Mahrt, 1986. A simple model of the atmospheric boundary layer: sensitivity to 

surface evaporation/ Boundary- Layer Meteorol., 37, 129-148.

Tuccillo, J.J., 1988. Parameterizations of physical processes in NMC's Nested Grid Model. 

Preprints, 8th Conf. on Numerical Weather Prediction, Feb. 22-26, Baltimore MD, 238-

243.

Van der Hoven, I., 1968. Deposition of particles and gases. In Meteorology and

Atomic Energy. D. Slade (Ed.), TID-24190, NTIS, Springfield, VA, 445 pp..

Walmsley, J. L., and M.L. Wesely, 1996. Modification of coded parameterizations of surface 

resistances to gaseous dry deposition. Atmos. Environ., 30, 1181-1188.

Wesely, M.L., 1989. Parameterizations of surface resistances to gaseous dry deposition in 

regional-scale numerical models. Atmos. Environ., 23, 1293-1304.

Wesely, M.L., and B.B. Hicks, 1977. Some factors that affect the deposition rates of sulfur 

dioxide and similar gases on vegetation. J. Air. Poll. Control Assoc., 27, 1110-1116.

Wilson, J.D., B.J. Legg, and D.J. Thomson, 1983. Calculation of particle trajectories in the 

presence of a gradient in turbulent velocity variance, B. Layer Meteorol., 27, 163-169.