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
2015, March
2018, February
2020, April
 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.5 Terrain following dry hydrostatic pressure-Sigma 4
2.5 Hybrid Absolute Pressure Sigma 5
2.6 Vertical Motion 5
2.7 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 10
4.4 Horizontal Mixing Coefficient 11
4.5 Particle Dispersion 11
4.6 Turbulent Velocity Variances 13
4.7 Convection 15
4.8 Puff Dispersion 16
4.9 Air Concentration Calculation 17

5. DEPOSITION 18

5.1 Gravitational Settling 19

5.2 Explicit Dry Deposition 19
5.3 Wet Removal 20
5.4 Radioactive Decay 21
5.5 Dry Deposition via Resistance Method 21
5.6 Pollutant Resuspension 24

6. SUMMARY 24

7. ACKNOWLEDGMENTS 24

8. REFERENCES 25
 DESCRIPTION OF THE HYSPLIT_4 MODELING SYSTEM

ABSTRACT. The HYSPLIT (HYbrid Single-Particle Lagrangian Integrated

Trajectory) model 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 cell-
average 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 Terrain following dry hydrostatic pressure-Sigma

Terrain following dry hydrostatic pressure-sigma coordinates are used in the Weather
Research and Forecasting (WRF ) model. This vertical coordinate system is defined as:

η = pdh – pdht / μd , where μd = pdhs – pdht , (13)

where pdh, pdhs, and pdht are the dry hydrostatic pressure at the location of interest, surface, and
model top, respectively. The change in height from one level to the next is defined by:
 Δz = - μd Δη αd / g , (14)

where αd is the layer average dry inverse density. Virtual temperature, relative humidity, and
saturation mixing ratio is calculated as described in the other coordinate systems above.

2.5 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, (15)

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.6 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), (16)

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.7 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 τ, (17)

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), (18)

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, (19)

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

and the final position is

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

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

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 and from these the Obukhov length is
calculated.

First the boundary layer depth (Zi) is computed at each grid point. There are five options
for estimating the boundary layer depth. The user may choose an option with the KMIXD
namelist variable in the [Link] file.
1) The default is to use the value provided by the meteorological model.
2) The boundary height may also be estimated from the temperature profile. In this case 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.
3) Boundary layer depth is calculated from the TKE profile provided by the
meteorological model. Zi is estimated to be at the height at which the TKE either decreases by a
factor of two or falls to a value of less than 0.21.
4) Boundary layer depth is set to a constant value. Users may wish to do this to match
external observations of the mixing depth.
5) Boundary layer depth is calculated using a modified Richardson number approach that
includes excess temperature for convective cases as a function of virtual potential temperature,
friction velocity, and convective vertical velocity following Vogelezang and Holtslag (1996).

If option 1 or 3 is chosen but the meteorological model does not provide the data
necessary, HYSPLIT will automatically switch to using option 2.

For option 5, mixing height is estimated by a modified bulk Richardson number (Ri)
approach that includes excess temperature for convective cases,

Ri = [g z (θv - θvs) ] / { θvsc [ (uz - us)2 (vz - vs)2 ] }, (23)

where θv and θvs are the virtual potential temperature at height z and the surface, respectively.
Under stable conditions, θvsc= θvs and under unstable conditions θvsc includes parcel excess
temperature:

θvsc = θvs + [8.5 u* T* (ρsl / ρz)] / [u*3 + 0.6 W*3]1/3 (24)

where ρsl and ρz is the density at the top of the surface layer and at height z, respectively. The
mixing height is set to height z where Ri equals a critical Richardson number (Ricr) of 0.25.

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, (25)

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

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

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

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

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, (29)

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, (30)

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, (31)

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, (32)

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

[2 β (1-βRb )]-1, 0≥Rb<0.08, (33)
z/L = (0.005 s +41.2) Rb2 + (1.18 s - 1.5 v - 1.37) Rb , Rb ≥ 0.08, (34)

and where

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

t = ln (z/Zh + 1.0), and (36)
v = ln (Zo/Zh). (37)

Here Zh is the roughness length for heat, Zo is the roughness length for momentum, and
Z0 / Zh = 10. The parameter β = 5. Equation (33) 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, (38)

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

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) ], (40)
φh = Prn {1 + z/L [ a ( 1+ a b z/L)½ + b exp(-d z/L) (1 - d z/L + c) ] }, (41)

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

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

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, (44)
wh = u* φh-1, z/Zi ≤ 0.1, (45)
wh = wm Pr-1, 1 > z/Zi > 0.1. (46)

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

wm = (u*3 + 0.6 W*3)1/3, (48)

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

Kh = - Cs u* T* (∂θ/∂z)-1, (50)

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 = l2 ∣ ∂V/∂ (ℓ/Lo)-1,
z ∣ hφ (51)

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

l-1 = kz-1 + 150-1, (52)

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

ℓ/Lo = 1.0893 Rib (53)

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

ℓ/Lo = a1 + Rib a2 + Rib [a3 + Rib (a4 + a5 Rib )] , (54)

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, (55)

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, (56)

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

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, (58)

where Δz is the vertical grid spacing, and hence

Δt = (Δz)2 / ( 8 σw2 TLw ). (59)

The vertical velocity variance, σw2, 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, (60)

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 , (61)

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 ), (62)

such that TLi= [TLw or TLu ], TLw (vertical) can be set to a constant, typically 5 s, and TLu
(horizontal) is set to a constant, typically 10800 s (∼1/f). These values result in a random walk
(R≈0) vertical dispersion for most of the longer time steps. TLw can also be computed to vary in
space in time following Hanna (1982). For stable conditions following Hanna (1982),

TLw = 0.1 zi/w’ (z/zi)0.8. (63)

For unstable conditions following Hanna (1982),

if z / zi < 0.1 and z - z0 > -L, TLw = 0.1 z / w’ [ 0.55 + 0.38 ( z - z0 ) / L ], (64)
if z / zi < 0.1 and z - z0 < -L, TLw = 0.59 z/w’, (65)
otherwise, TLw = 0.15 zi / w’ [ 1 – exp( -5 z / zi ) ], (66)

where z0 is the aerodynamic roughness length.

 The Gaussian random component U" or W" comes from the computer generated random
number, such that for any directional component,

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

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. 44, 50, 51) and horizontal diffusivities (Khor - Eq. 55),

σi = (Ki / TLi )0.5, (68)

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

HYSPLIT can also use the STILT dispersion scheme (Lin et al., 2003), which utilizes a very thin
layer above the mixed layer, a fractional timestep, and the Thomson (1997) transmission /
reflection scheme. Within this option, a very thin layer is created just above the mixed layer to
prevent particles from getting trapped in a low-turbulence layer. A fractional timestep is
employed where the timestep stops and a new one begins every time a particle is transported
vertically to a model level interface. When a particle is at the model level interface, the Thomson
(1997) transmission / reflection scheme is used to preserve a well-mixed distribution of particles
across the interface. For an upward moving particle that reaches a model interface with a vertical
turbulent velocity W’i, it is transmitted with velocity W’tup if the ratio, αup, is greater than a
random number between 0 and 1, otherwise it is reflected from the interface with vertical
velocity –W’i, where,

W’tup = W’i [ σw(zi+) / σw(zi-) ], αup = [ σw(zi+) ρ(zi+) ] / [ σw(zi-) ρ(zi-) ], (69)

where ρ is density and zi+ and zi- represent the interface above and below the current grid cell.
For a downward moving particle that reaches a model interface, the particle is transmitted with
velocity W’tdn if the ratio, αdp, is greater than a random number between 0 and 1, otherwise it is
reflected with a velocity –W’i, where:

W’tdn = W’i [ σw(zi-) / σw(zi+) ], αdn = [ σw(zi-) ρ(zi-) ] / [ σw(zi+) ρ(zi+) ]. (70)

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. 68). This permits the
use of the meteorological model’s TKE (turbulent kinetic energy field) if it is available.
Turbulent vertical velocities can also be computed following Kantha Clayson (2000) or Hanna
(1982).

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 = 1.7 u*2 (1 – z/zi)3/2, (71)
 u’2 = 4.0 u*2 (1 – z/zi)3/2, (72)
v’2 = 5.0 u*2 (1 – z/zi)3/2. (73)

In the stable or neutral surface layer,

w’2 = 1.7 u*2, (74)
u’2 = 4.0 u*2, (75)
v’2 = 5.0 u*2. (76)

If the user defined nighttime TKE turbulence anisotropy factor (TKERN) is greater than the
daytime turbulence anisotropy factor (TKERD),

TKER{D|N} = w’2 / (u’2 v’2), (77)

then the vertical velocity variance is increased by 40% during stable nighttime conditions.

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

If the user defined TKE turbulence anisotropy factors are set to a value greater than zero, then
velocity variances are defined as,

w’2 = 2.0 E / (1.0 + 1.0 / TKER), (79)

u’2 = v’2 = 0.5 TKER w’2. (80)

Otherwise, velocity variances for the TKE scheme are defined consistent with the TKE
turbulence parameterization, as described below. Velocity variances in the stable or neutral
boundary or surface layer are defined following constant relationships to the TKE:

w’2 = 0.32 E, (81)

u’2 = 0.74 E, (82)
v’2 = 0.85 E. (83)

Vertical velocity variances calculated using the Hanna scheme (1982) are derived from u*, w*,
and Zi. For stable conditions using the Hanna (1982) scheme,

w’2 = [ 1.3 u* (1 - z/zi) ]2, (84)

For unstable conditions following Hanna (1982),

if z/zi < 0.3, w’2 = [ 0.96 w* (3.0 z/zi – L/zi)1/3 ]2, (85)
if z/zi > 0.3 and z/zi < 0.4, w’2 = {w* × MIN[0.96 (3.0 z/zi – L/zi )1/3,
0.763 (z/zi)0.175]}2, (86)
if z/zi >= 0.4 and z/zi < 0.96, w’ = [0.722 w* (1 – z/zi)0.207 ]2,
2
(87)
if z/zi > 0.96, w’2 = 0.37 w*2, (88)

When the Hanna (1982) scheme is used to calculate vertical velocity variance, horizontal
velocity variance is defined as,

u’2 = v’2 = 0.5 w’2, (89)

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), (90)

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

In the unstable surface layer,

w’2 = 1.74 u*2 (1 – 3 z/L)2/3, (92)

u’2 = v’2 = 0.36 w*2. (93)

In the unstable surface or boundary layer, if the TKE field is available, then from the definition
(78) and Eqs. 90-93,

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

u’2 = v’2 = 0.5 Kx w’2, (95)

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

and in the unstable boundary layer,

Kx = 0.72 / [ (z/zi)2/3 (1 – z/zi)2/3 (1 + 0.5 R 2/3) ]. (97)

Unlike the diffusivity (Eq. 44) approach, the vertical velocity variance is not averaged for
a single boundary layer value. The turbulence varies by height according to Eqs. 71-95. The un-
averaged turbulent velocities are much more realistic for short-range dispersion simulations.
Note that the horizontal turbulent velocity also replaces the deformation based diffusvity (Eq.
55). 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 Convection

Vertical transport due to convection can be employed with HYSPLIT using one of three
options: 1) CAPE threshold method; 2) extreme convection method; or 3) Grell convection
scheme. With the CAPE threshold method, enhanced mixing due to convection is employed
when CAPE exceeds a user defined threshold. CAPE enhanced mixing results in particles in the
cloud-layer being randomly redistributed within the cloud-layer if they reside in a grid cell where
CAPE exceeds the user defined minimum CAPE threshold. This scheme is based on the extreme
convection described below, but only takes place when CAPE exceeds the user defined
threshold.

The extreme convection method vertically mixes all particles in grid cells with positive
CAPE upward throughout the entire unstable layer defined by the limit of convection (Gerbig et
al., 2003). This scheme randomly assigns a vertical position between the surface and the limit of
convection altitude (ZLOC) to each particle with a vertical position below ZLOC. The random
redistribution is weighted by density. This scheme assumes that updrafts nad downdrafts are
large enough to leave a perfectly well mixed column behind each convective event. This scheme
is designed to provide an upper limit of the impact of subgrid-scale vertical redistribution due to
convective cloud transport. A lower limit can be achieved by running HYSPLIT with no
convective parameterization.

The Grell convection scheme within HYSPLIT ingests convective mass fluxes from the
meteorological input files that were generated from Weather Research and Forecasting (WRF)
model output files that were run with the Grell et al. (1994) or Grell and Devanyi (2002)
convection schemes. These include vertical profiles of updrafts, downdrafts, detrainment, and
entrainment fluxes. Vertical profiles of up- and downward vertical velocity are derived from the
flux profiles and grid cell fractional coverage of the up- and downdrafts. The vertical profiles of
mass fluxes of updrafts, downdrafts, entrainment, and detrainment are used to compute the
probability of particles being located within the updrafts, downdrafts, or outside of the
convective system.

4.8 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σz2/dt = 2 σw2 TL, (98)

and there are two options for the horizontal growth rate,

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

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

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,

σzt2 (t+Δt) = σz2(t) + 2 σwt2(t) TL Δt, and (100)

σzb2 (t+Δt) = σz2(t) + 2 σwb2(t) TL Δt. (101)

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)]. (102)

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, (103)

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.9 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, (104)

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 π σh2 Δz)-1 exp(-0.5 x2/σh2), (105)

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 grid-
cell 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, (106)

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. (104), but with Δz defined as grid-cell height. If the
horizontal distribution is Gaussian then the incremental concentration is the same as Eq. (105),
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. (22) 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 (ℓ/ℓ) 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 ) ] }, (107)

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 = dp2 g (ρg - ρ) (18 μ)-1, (108)

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 non-
spherical particles. Hence the final settling velocity is

Vs = Vg Cc α-1, (109)

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 Λ) ]}. (110)

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 / ρ). (111)

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

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, (113)

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, (114)

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

βinc = Ft Fb Vinc ΔZp-1, (115)

where the average scavenging ratio is S= (5x105 to 1x106 ) by volume, and Δ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 (when the removal coefficient is defined to be less than one), similar to
the below-cloud procedure,

βinc = Ft Fb (4x10-5 to 8x10-5 ) P0.79. (115a)

Below-cloud removal is defined directly as a rate constant, independent of the precipitation rate
when using Eq. 115 or a function of P (mm/h) when using Eq. 115a for within-cloud removal.
The below-cloud removal constant (s-1) is given by when using Eq. 115,

βbel = 1x10-6 (1.0-Fb), or when using Eq. 115a, (116)

 βbel = (4x10-5 to 8x10-5 ) P 0.79 (1.0-Fb). (116a)

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, (117)

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 = Ft Vgas ΔZp-1. (118)

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

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 T½,

βrad = ln 2 / T½ (119)

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

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

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, (121)

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*. (122)

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*), (123)

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

where φh was given by Eqs. (41) and (43), 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 (125)

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

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

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

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*) Scd2 ]-1 + u*[ St / (St+p) ]q }-1, (129)

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*) (ℓ* g Cc)-1, (130)

where ℓ* is the laminar layer length scale (ℓ*= υ/u*) and the Schmidt number is given by
 Sc = υ / D. (131)

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)½. (132)

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

D = kb T Cc ( 3 π μ dp)-1, (133)

where kb equals Boltzmann's constant (1.38x10-20 g2 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, (134)

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))], (135)

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, (136)

Rlu = Rlu1 [ H* 10-5 + fo ]-1, (137)
Rdc = 100 [ 1+1000/(G+10) ], (138)
Rcl = [ H* 10-5/Rcl1 + fo/Rcl2 ]-1, (139)
Rgs = [ H* 10-5/Rgs1 + fo/Rgs2 ]-1, (140)
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. 134) 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 (sm-
1
)

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

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, (142)

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, (143)

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

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. (144) is emitted. The deposition total at that
cell is reduced accordingly.
6. SUMMARY

A detailed description and equations of HYSPLIT, 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.

Gerbig, C.J., J.C. Lin, S.C. Wofsy, B.C. Daube, A.E. Andrews, B.B. Stephens, P.S. Bakwin, and
C.A. Grainger, 2003. Toward constraining regional-scale fluxes of CO2 with atmospheric
observations over a continent: 2. Analysis of COBRA data using a receptor-oriented
framework. J. Geophys. Res., 108, 4757, [Link]

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

Grell, G.A., J. Dudhia, and D.R. Stauffer, 1994. A description of the fifth generation Penn
State/NCAR mesoscale model (MM5), NCAR Tech Note NCAR/TN-398+STR, pp. 64-
72, NCAR, Boulder, CO, [Link]

Grell, G.A. and D. Devenyi, 2002. A generalized approach to parameterizing convection

combining ensemble and data assimilation technizues. Geophys. Res. Lett., 29, 1693,
[Link]

Hanna, S.R., 1982. Applications in air pollution modeling. Atmospheric Turbulence and Air
Pollution Modelling. F.T.M. Nieuwstadt and H. Van Dop, Eds., D. Reidel, 358 pp.

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.

Lin, J.C., C. Gerbig, S.C. Wofsy, B.C. Daube, A.E. Andrews, K.J. Davis, and C.A. Grainger,
2003. A near-field tool for simulating the upstream influence of atmospheric observations:
The Stochastic Time-Inverted Lagrangian Transport (STILT) model. J. Geophys. Res.,
108, 4493.

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.
Thomson, D.J., W.L. Physick, and R.H. Maryon, 1997. Treatment of interfaces in random walk
dispersion models. J. Appl. Meteor., 36, 1284-1295.

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

Vogelezang, D.H.P. and A.A.M. Holtslag, 1996. Evaluation and model impacts of alternative
boundary-layer height formulations. B. Layer Meteorol., 81, 245-269.

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.