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 L aboratory Silver Spring, Maryland

December 1997

Revised: 1998, August 2002, September
2003, October 2004, January
2004, November
2010, December 2014, SeptemberNOTICE

Mention of a commercial company or product does not constitute

an endorsement b y 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 Surfa ce 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 Disper sion 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 22DESCRIPTION OF THE HYSPLIT_4 MODELING SYSTEM

ABST RACT. 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 use s 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 base d 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 numeric al 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 o f 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 a ccording 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 t o use differ ent 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,σ = (Z top - Z msl ) / (Z top - Z gl ), (1)

where al l the heights are expressed relative to mean

- sea level and where Z top 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 = ak 2 + bk +c, (2)

where a=30, b=

- 25, and c=5. This relation results in decreasing resolution away from the
surface, w ith the first level (k=1) at 10 m, the second level at 75 m, the third at 200 m, while by

the 20 th 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 b etween 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 f our 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 horizo ntal wind components), T (temperature), Z (height) or P (pressure), and the pressure at the surface, P 0 . Moisture and vertical motion are optional, however the vertical motion may be computed based upon how the vertical coordinate is defined. If wet depo sition 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

Whe n the input data are given on pressure surfaces, the only surface value required is pressure. Heights of the pressure surfaces (Z p ) are assumed to be relative to msl (Mean Sea Level). Model output fields on absolute pressure surfaces have usually been in terpolated to those surfaces. Moisture output is expected as relative humidity (RH) while temperature is assumed to be a virtual temperature (T

v ). Surface values are computed if not otherwise provided. For instance, the surface temperature (T 0 ) is built down adiabatically from the lowest data level using the input data level

- 1 temperature (T 1 ) and pressure (P 1 ),

T 0 = T 1 (P 1 /P 0 ) - 0.286

. (3)

The height of ground surface (Z g ) is interpolated directly from the profile if the ground level fallsbetwe en input data levels. Otherwise the ground surface is estimated from the equation of state,

Z g = Z 1 - R d T 01 ln (P 0 /P 1 ) g - 1

, (4)

where R d is the gas constant for dry air (287.04 J kg - 1 K - 1 ), T 01 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 (U model ) are estimated at all model levels (Z model ) below the lowest data level (Z data ) using a logarithmic profile for
neutral conditions,

U model = U data ln (Z model /Z 0 ) / ln (Z data /Z 0 ), (5)

where Z 0 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 m odel 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,

ρ = k 1 P T - 1 R d - 1

, (6)

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

p otential temperature

θ = T (1000/P) 0.286

. (7)

If some measure of vertical motion (W p ) is provided (usually in units of hPa s - 1 ), it is converted
to (W

σ ) sigma coordinates (s - 1 ),

W

σ = k 1 W p ρ - 1 g - 1 Z top - 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 interpola ted 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 a ssumed to be a dry temperature. Therefore the virtual temperature is first computed from the specific humidity (Q)

T

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

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

pressure

- s igma 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(P 1 /P 2 ) R d T v12 g - 1

. (10)

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

RH ≈ Q P ( 0.622 E

s ) - 1

, (11)

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

E

s = 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 conversi ons. 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 sp atial 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

coord inate 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 p ressure

- sigma data as described previously. Input data are defined

on hybrid pressure levels given by the relation

P k = A k + B k P 0 , (13)where k indicates the index of the level, P the pressure, P 0 the surface pressure, A a pressure
offset, a nd 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 ve rtical 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 sy stem 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 th ere 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, 1 972) 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 e t 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 computation al 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 ha ve 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 t he 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 int egration 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

th e grid spacing. The maximum transport velocity U

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

U

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

duratio n 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 motio n 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 thesecond 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 a s 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 meteorolo gical 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 dis tances 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 (Z i ) 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 n ear 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 = ρC D u where C D 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 (ρ C p u * ) - 1

, (23)

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

W * = ∣ g u * T * Z i T - 1 ∣ 1/3

, (24)

where Z i 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 = Z 2 k g T * (u * 2 T 2 ) - 1

, - 2 ≤ z /L ≤ 10, (25)and where Z 2 indicates height of 2 nd 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,

R b ’ = 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 R b ’ is adjusted to account for those meteorological data files that have very coarse vertical spacing, where the m eteorological values at the top of the surface layer (Z 2 ) have been extrapolated from input data levels (Z d ) at much greater heights. In those special situations

we assume that

R b ≈ R b ’ (Z 2 /Z d ) 2

, (27)

to compensate for the increase of R b ’ 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 est imated from R b . 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 = R b (s 2 /t

- 0.50), R b ≤ 0, (28)

z/L = [

- t+2s βR b + (t 2

- 4s t βR b + 4s 2 βR b ) 1⁄2 ]

[2 β (1

- βR b )] - 1

, 0≥R b <0.08, (29)

z/L = (0.005 s +41.2) R b 2 + (1.18 s

- 1.5 v

- 1.37) R b , R b ≥ 0.08, (30)

and where

s = ln (z/Z

o + 1.0), (31)

t = l n (z/Z h + 1.0), and (32)

v = ln (Z

o /Z h ). (33)

Here Z h is the roughness length for heat, Z

o is the roughness length for momentum, and

Z 0 / Z h = 10. The parameter β = 5. Equation (29) is an exact solution for a log

- linear profile (Choudhu ry, et al. 1986). The friction velocity and temperature are then given by

u * = k Z 2 Δu (φ

m ΔZ)

- 1

, (34)T * = k Z 2 Δθ (φ 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 = Pr

n {1 + z/L [ a ( 1+ a b z/L) 1⁄2 + b exp(

- d z/L) (1

- d z/L + c) ] }, (37)

where Pr

n is the turbulent Pra ndtl 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

K h = k w h z (1

- z/Z i ) 2

, (40)

w h = u * φ h - 1

, z/Z i ≤ 0.1, (41)

w h = w m Pr - 1

, 1 > z/Z i > 0.1. (42)

The Prandtl number, with φ h and φ

m evaluated at z/Z i = 0.1, is given by

Pr = (φ h / φ

m ) + 7.2 k ( z/Z i ) (W * /w

m ), where (43)

w m = (u * 3 + 0.6 W * 3 ) 1/3

, (44)

and W * =0 fo r 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 = Z i during convective conditions (W * >0) is computed based upon the surface flux parameters and the strength of the inversion, where

K h = C

s H ( ρ C p ∂θ/∂z )

- 1 , or (45)

K h =

- C

s u * T * (∂θ/∂z)

- 1

, (46)and C

s = 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

verti cal diffusivity for heat using mixing length theory where

K h = l 2 ∣ ∂V /∂z ∣ φ h (l/L

o ) - 1

, (47)

where l is a Blackadar

- type mixing length in meters,

l - 1 = kz - 1 + 150 - 1

, (48)

and L

o 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/L

o = 1.0893 Ri b (49)

during near

- neutral conditions (0 ≤ Ri b ≤ 0.001) and otherwise ( 0.001 ≤ Ri b ≤ 20 )

l/L

o = a 1 + Ri b a 2 + Ri b [a 3 + Ri b (a 4 + a 5 Ri b )] , (50)

and where a 1 =0.2828x10 - 3 , a 2 =0.8049, a 3 =1.6583, a 4 =0.5090x10 - 2 , and a 5 =

- 1.0063x10 - 3

.

4.4 Horizontal Mixing Coefficient

The subgrid

- scale horizontal mixing coefficient is computed from the velocity deformation (Smagorinsky, 1963; Deardorff, 1973),

K hor = 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 Dis persion

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 approa ch 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),

X final (t+Δt) = X

mean (t+Δt) + U'(t+Δ t) Δt G, (52)

Z final (t+Δt) = Z

mean (t+Δt) + W'(t+Δt) Δt Z top - 1

, (53)

where the horizontal and vertical positions are given in grid and sigma units, respectively, whilethe turbulent velocity components are in m s - 1 . G and Z top are the required u nit 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. F ull 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,

Δz p < 0.5 Δz, (54)

where Δz is the vertical grid spacing, and hence

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

w 2 T Lw ). (55)

The vertical velocity variance, σ

w 2 , and the Lagrangian time scale, T Lw , are discussed in more

detail below.

The horizontal turbulent v elocity 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

+ T Lw (1

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

and as defined by Wils on 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 / T Li ), (58)

such that T li = [T Lw or T Lu ], and wher e T lw (vertical)=100 s and T lu (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 t he 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 estima ted from the previously calculated vertical
(K h - Eqs. 40, 46, 47) and horizontal diffusivities (K hor - Eq. 51),σ i = (K i / T Li ) 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 Z i following Kantha and Clayson (2000). This method does not use the diffusivity and hence no

assumptions are required about turbulent scales. Ei ther 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/z i ) 3/2

, (61)

u’ 2 = 4.0 u * 2 (1

– z/z i ) 3/2

, (62)

v’ 2 = 4.5 u * 2 (1

– z/z i ) 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 (20 00),

w’ 2 = w * 2 (z/z i ) 2/3 (1

– z/z i ) 2/3 (1 + 0.5 R 2/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 definitio n (67) and Eqs. 71

- 74,

w’ 2 = 2 E / (1+K

x ), (75)

u’ 2 = v’ 2 = 0.5 K

x w’ 2

, (76)

and where in the unstable surface layer K

x is defined to be

K

x = 0.41 (w * 2 / u * 2 ) / (1

– 3 z/L) 2/3

, (77)

and in the unstable boundary layer,

K

x = 0.72 / [ (z/z i ) 2/3 (1

– z/z i ) 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 un

averaged t urbulent 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 comp uted 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 ar e 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 T L , (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 T L / 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 ar e 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 to p

- 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σ) su ch that,σ zt 2 (t+Δt) = σ

z 2 (t) + 2 σ wt 2 (t) T L Δt, and (81)

σ zb 2 (t+Δt) = σ

z 2 (t) + 2 σ wb 2 (t) T L Δ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 tur bulent 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 > L h , at positions P(x±0.5σ h , y±0.5σ h ), and where L h is the meteorological model grid size. A large Gaussian puff splits into five smaller puffs when 3.0σ h > L h , 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, P

n (z), at the center height of “n” l ayers, 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 a re 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, limit s 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 def ined 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 e xtent 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 contributi on by each puff of mass m to a grid point is computed as follows for a top

- hat puff,

Δc = m (π r 2 Δ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 sam e Δc. The incremental concentration contribution for aGaussian puff is,

Δc = m (2 π σ h 2 Δz)

- 1 exp(

- 0.5 x 2 /σ h 2 ), (86)

where x is the distance from the puff center to grid

- node, and the other terms are as previously

defined.

Particle calculatio ns 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 in to 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

, (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. I f 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 pa rticles 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 la rge 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 become s

D wet+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 gase s (β gas ), in

- cloud wet removal of particles (β inc ), and below

- cloud wet removal ofparticles (β 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

als o 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 p article will be deposited.

5.1 Gravitational Settling

Particle settling is computed following Van der Hoven (1968), where the settling velocity
(V g ) is calculated for a spherical particle from the particle diameter (d p ), air density (ρ), and particle de nsity (ρ g ),

V g = d p 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 (C

c ) for small particles and a dynamic shape factor (α) to account for non

sph erical particles. Hence the final settling velocity is

V

s = V g C

c α - 1

, (90)

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

C

c = 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 me teorological data level. Then the mass deposited by dry removal

D dry = V d 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, ΔZ p , equals ±1.54 σ

z . In the particle calculation,

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

constant of the form

β dry = V d ΔZ p - 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 a lternative 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 laye r 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 r ain 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 pollu tant mass that is within the cloud layer by defining the fraction of the pollutant layer that is below the cloud top (F t ) and the fraction of the pollutant layer that is above the cloud bottom (F b ).

For the wet removal of particles by within

- cloud process es 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,

V inc = S P

, (95)

where the precipitation rate is given by P . The time constant for withi n

- cloud removal,

β inc = F t F b V inc ΔZ p - 1

, (96)

where the a verage scavenging ratio is S=5.0 x10 5 by volume, and as before, ΔZ p is the depth of the pollutant layer for puffs and the depth of the cloud layer for particles . Different scavenging

rat ios 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 = F t F b 8 x10 - 5 P 0.79

. (96 a )

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 = 8 x10 - 5 (1.0

- F b ), or (97)

β bel = 8 x10 - 5 P 0.79 (1.0

- F b ). (97 a )

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 equilibriumconcentration in water to that in air. A gaseous wet deposition velocity can be de fined as

V gas = 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 = F t V gas ΔZ p - 1

. (99)

Note that the wet removal of gases is applied at al l 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 deca y each time step. The decay constant for radioactive processes is defined by the half

- life T 1⁄2 ,

β rad = ln 2 / T 1⁄2 (100)

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

becomes

m 2 = m 1 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 settl ing velocity for particles such that the total deposition velocity

V d = [ R

a + R b + R

c + R

a R b V g ] - 1 + V g , (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, V g , is zero for gases and R

c is zero for particles. The resistance components depend upon meteorological conditions as well as the properties of the surface. The

surfac e 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 incorporate d 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 Hic ks,
1977)

R

a = Pr

n (ln z/Z

o - ψ h ) / k u *

. (103)

Over land the aerodynamic roughness length Z

o is constant and determined from the land

- useand vegetative cover; over water Charnock's (1958) relation, as modified by Smith (1988),

Z

o = 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 = a 1 + z/L (a 2 + z/L (a 3 + z/L ( a 4 + a 5 z/L))),

- 2 ≤ z/L ≤

- 0.001 (107)

where a 1 = 0.1164x10 - 4 , a 2 =

- 2.7188, a 3 =

- 2.1551, a 4 =

- 0.9859, and a 5 =

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

c ) number following Wesely and Hicks (1977)

R b = Pr (d 1 /ku * ) S

c d 2

. (109)

In the above relations the constants d 1 =2, d 2 =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)

R b = Pr { [ (d 1 / ku * ) S

c d2 ] - 1 + u * [ S t / (S t +p) ] q } - 1

, (110)

where the constants p=0.8 and q=2.0 (Peters and Eiden, 1992), and the Stokes number S t is computed from

S t = (2 V g u * ) (l
* g C

c ) - 1

, (111)

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

S

c = υ / D. (112)

The diffusivity of a specific gaseous pollutant, D, is related to the ratio of the molecular weights of the pollutant (W p ) and air (W

a ) through Graham's Law,D = υ (W

a /W p ) 1⁄2

. (113)

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

D = k b T C

c ( 3 π μ d p ) - 1

, (114)

where k b 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 c haracteristics that control the uptake of gases into plants and act in parallel. As
noted earlier, R

c is zero for particles. The following procedure follows the equations and notation as outlined by Wesely (1989) for the total canopy resistance,

R

c = [ 1/(R

s +R

m ) + 1/R lu + 1/(R dc +R cl ) + 1/(R

ac +R gs ) ] - 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

fac tor (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,

R

s = R i D hx [1 + {200/(G+0.1)} 2 ] [400 / (T

s (40

- T

s ))], (116)

where G is the solar irradiation in W/m 2 , D hx is the ratio of the diffusivity of water vapor to that of the pollutant (Table 2 in Wesely, 1989), T

s is the ambient temperature in degrees Celsius, and

R i 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, R

s 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 solubilit y and reactivity of the pollutant and a simple fractional expression can be used for each:

R

m = [ H * /3000 + 100 f

o ] - 1

, (117)

R lu = R lu1 [ H * 10 - 5 + f

o ] - 1

, (118)

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

R cl = [ H * 10 - 5 /R cl1 + f

o /R cl2 ] - 1

, (120)

R gs = [ H * 10 - 5 /R gs1 + f

o /R gs2 ] - 1

, (121)

where H * is an effective (relative to SO 2 ) Henry's constant and f

o is a pollutant specific reactivity parameter. All constants (including R

ac 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 uponland

- use and season. The relationship for R dc has been simplified by assuming that the slope of

the terrain is zero and the relation for R lu applies only to dry surfaces. In addition the term (sm 1 )

1000 exp(

- T

s - 4), (122)

is added to R lu , R cl , and R gs , to account for cold

- temperature resistance increases.

5.6 Pollutant Resuspension

Under certain conditions pollutants that have deposite d 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 pollut ant 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 proce ss 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. SUM MARY

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 varie ty 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 simu lation 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 int erface for the BoM meteorological data.

8. REFERENCESAbdella, 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.

Bel jaars, 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 unstab ly 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 Processe s, 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 L agrangian 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.]

D raxler, 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.

Dr axler, 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 Me dium

- 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. Bovill e, 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 Processe s 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 systema tic errors in the MASS 2.0

model. Mon. Weather Rev. , 113, 1714

- 1737.

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

Launiainen, J., 1995. Derivation of the relationship between the Obukhov stabili ty 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 g radient 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

func tion 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 atmo spheric 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

- 2 6, 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 parameter izations 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.

Wesel y, 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 traject ories in the presence of a gradient in turbulent velocity variance, B. Layer Meteorol., 27 , 163

- 169.