## Introduction

Lagrangian models have some attractive features over Eulerian models when diffusion and transport of airborne materials from a point source are considered. Eulerian models obtain numerical solutions of the finite-difference forms of the mass conservation equations. The finite-difference method is known to introduce large numerical diffusivity, the magnitude of which could be much larger than the physical diffusivity. The magnitude of numerical diffusivity is proportional to the grid spacing used in the finite-difference equations. Typical grid spacings used in air quality models are 1–10 km, which are very large in comparison with the typical size of a point source. Thus, Eulerian models are not appropriate for predicting concentration distributions that originate from a point source.

On the other hand, Lagrangian models are free from numerical diffusivity and maintain mass conservation. Lagrangian models have been used extensively for prediction of transport and diffusion of airborne materials from a point source over homogeneous boundaries. Gaussian plume and puff models are the most widely used models in this category. Lagrangian models also have been applied over inhomogeneous boundaries such as complex terrain. In such applications, pseudoparticles are released to mimic the movement of airborne materials. Trajectories of particles are computed by using mean and turbulence velocities and the concept of random walk.

The Higher-Order Turbulence Model for Atmospheric Circulation–Random Puff Transport and Diffusion (HOTMAC–RAPTAD) modeling system (Yamada and Bunker 1988) has been applied successfully to simulate wind, turbulence, and distributions of pollutant concentration over complex terrain where conventional methods had failed. HOTMAC is a three-dimensional mesoscale model, which is based on a set of second-order turbulence-closure equations. RAPTAD is a three-dimensional Lagrangian puff model, which is based on a random displacement method for which mean and turbulence components of winds are provided by HOTMAC. Successful simulations of concentration distributions depend greatly on the accuracy of wind directions, wind speeds, and atmospheric turbulence used in pollutant transport and diffusion models.

Previous RAPTAD applications were limited to scalar transport during which plumes follow the ambient air motion exactly. There are cases for which plume density differs considerably from the surrounding air density. Under these conditions, plumes do not follow the ambient air motion. When the plume density is less than the surrounding air density, it is called a buoyant plume. On the other hand, if the plume density is greater than the surrounding air density, it is called a dense (heavier than air) gas situation.

The purpose of the current study was to incorporate the nonneutrally buoyant plume algorithm developed by Van Dop (1992) into RAPTAD and to evaluate the model results by using the data and the statistical parameters described by Hanna et al. (1991, 1993). The motivation of the study was to improve the forecasting capability of the transport and dispersion of airborne hazardous materials at Vandenberg Air Force Base (VAFB). On occasion, it is necessary to abort a rocket flight prematurely. Destroying a rocket in the atmosphere releases hazardous materials at different altitudes. The source distributions could extend over a kilometer in the vertical direction. Chemical transformation of nitrogen tetroxide, which is used in rocket propellants at VAFB, is a complex process and is beyond the subject of this study. HOTMAC and RAPTAD, however, can be used to forecast concentration distributions if source distributions and characteristics are specified when initial chemical transformation is completed.

Except for a few models, dispersion models available today are “specialized” or “tuned” to simulate either a neutral, a buoyant, or a dense gas plume. Many regulatory models used for environmental impact statements are limited to a homogeneous environment (constant winds and turbulence). The authors have been developing a modeling system that can be used to forecast concentration distributions of nonneutrally buoyant plumes in an inhomogeneous environment, such as airflows over complex terrain.

One of the most difficult tasks for a model developer is to identify observations that can be used for model performance evaluations. A comprehensive set of meteorological and buoyant tracer gas data fortunately have become available. Project MOHAVE (Measurement of Haze and Visual Effects) conducted an extensive field campaign over the arid southwestern United States (Pitchford et al. 1997). The study area included Mohave Power Project, which is a large coal-fired power plant in southern Arizona, and the Grand Canyon National Park. The HOTMAC–RAPTAD modeling system with the van Dop buoyant plume algorithm discussed in this paper was applied to simulate perfluorocarbon tracer gas concentrations at eight sampling locations whose distances from MPP varied from a few tens of kilometers to over 200 km. The model performance was evaluated statistically (Yamada 2000).

We are not aware of observations that can be used to test the proposed dense gas algorithm for transport and diffusion over complex terrain. Only limited documentation, consisting mostly of visual observations, is available for accidental releases associated with rocket abortions. Visual observation data are useful for qualitative evaluations of modeled plume locations, but concentration measurements are essential for statistical evaluation of model performance. This requirement is why we have chosen to use the dense gas concentration data reported in the Modelers’ Data Archives (MDA; Hanna et al. 1991, 1993), although MDA was collected over flat terrain, and sampling stations were located relatively close to a source.

## Model

### RAPTAD

#### Transport

Highly buoyant plumes, such as a fire plume, modify the wind and turbulence distributions of the ambient flow. Similarly, highly dense gas released near the ground under calm conditions (no winds or turbulence) produces its own motion, converting downward motion to horizontal spreading at the ground. Furthermore, dense gas could alter density stratification of the atmosphere, resulting in modification of atmospheric turbulence.

It is almost impossible to parameterize or express such modifications without deploying a dynamic plume model. A physically correct way would be to incorporate plume density into a hydrodynamic model. However, a dynamic plume model requires considerable computer time and is not practical for operational uses.

An alternative approach, to modify a dispersion model to include plume density, has the advantage of being faster. The disadvantage is that it does not consider the effects of plume density on the ambient airflows. In other words, the ambient airflows are considered to be unchanged even if plume density modifies both wind and temperature distributions. This approach is justifiable only when ambient air modifications are confined to a small area, and the magnitude of the modifications is small. Such an approximation is considered here.

*x*

_{i}

*t*

*t*

*x*

_{i}

*t*

*U*

_{pi}

*t,*

In the above expressions, *U*_{pi} is the puff velocity in the *x*_{i} direction, *U*_{i} is the mean velocity, *u*_{i} is the turbulence velocity, *ζ* is a random number from a Gaussian distribution with zero mean and unit variance, *t* is time, Δ*t* is the time increment, *t*_{Lxi}*u*_{i}, and *σ*_{ui} is the variance of *u*_{i}. The mean velocity and *σ*_{ui} are obtained from HOTMAC (section 2b).

Gifford (1982) suggested that the horizontal timescales should be *O*(*f*^{−1}), where *f* is the Coriolis parameter (10^{−4} s^{−1} at 40°N). Thus, the Lagrangian timescales *t*_{Lx}*t*_{Ly}*t*_{Lx}*t*_{Ly}

Eddies in the vertical direction in the atmospheric boundary layer are considered to be much smaller than the eddies in the horizontal directions. The eddy size in the convective (unstable) boundary layer may be limited by the mixed-layer height (a few km at most), and the eddy size in the stable boundary layer may be suppressed by the density stratification and become very small (less than 10 m). Thus, eddy size varies considerably from daytime (convective) to nighttime (stable).

Turbulence intensities also exhibit diurnal variations. Turbulence become large during the convective period, but turbulence is suppressed during the nocturnal period. If we assume that the eddy size represents a length scale, and that turbulence kinetic energy represents a velocity scale, a timescale may be proportional to the ratio of the length scale to the velocity scale. In this way, the Lagrangian timescale in the vertical direction *t*_{Lz}

Legg and Raupach (1982) introduced an additional term (1 − *a*)*T*_{Lz}(∂*σ*^{2}_{w}/∂*z*)*σ*^{2}_{w}

For a neutrally buoyant plume, a simple reflection algorithm has been used when a puff hits the ground. The same boundary condition was used for the dense gas.

#### Plume density

*W*

*B*

*Z*

In Eq. (9), the subscripts *p* and *a* indicate plume and ambient air, respectively. Gravitational acceleration is *g,* *T* is temperature, ϒ is potential temperature, *z* is height, and *t* is time. Here *T*_{W} and *T*_{B} are relaxation timescales for *W**B*

*ρ*

_{p}is the puff density, and

*ρ*

_{a}is the ambient air density. The original application of van Dop’s (VD) model was limited to buoyant plumes. As far as the methodology is concerned, however, VD should be applicable not only to buoyant plumes but also to dense gas releases.

*U*

_{i}and turbulence

*u*

_{i}. The vertical velocity used in RAPTAD for a nonneutrally buoyant plume was obtained by adding the HOTMAC velocity (

*W*

_{HT}) and VD’s velocity (

*W*

_{VD}), that is,

*W*

_{puff}

*W*

_{HT}

*W*

_{VD}

*W*

_{VD}is the solution of Eqs. (6)–(8). Equations (6)–(8) were solved numerically since

*N*

^{2}in Eq. (10) was not a constant. Otherwise the solution of Eqs. (6)–(8) can be obtained in analytical form.

#### Concentration

Various functional forms can be assumed to express the concentration distribution in the puff. One of the simplest ways is to assume a Gaussian distribution, where variances are determined as the time integration of the velocity variances encountered over the history of the puff. The concentration level at a given time and place is calculated as the sum of the concentrations from each puff.

*C*at (

*X, Y, Z*) is computed by using the following expression:where (

*x*

_{k},

*y*

_{k},

*z*

_{k}) is the location of the

*k*th puff;

*σ*

_{xk}

*σ*

_{yk}

*σ*

_{zk}

*z*

_{g}is the ground elevation. The variances are estimated based on Taylor’s (1921) homogeneous diffusion theory. For example,

*σ*

_{y}is obtained fromwhere a correlation function

*R*(

*ξ*) =

*ξ*/

*t*

_{Ly})

It may appear that the current method counts the effect of turbulence twice: once in random motion of puffs [Eq. (3)] and again in broadening the kernel bandwidth [Eq. (17)]. The readers are referred to Yamada and Bunker (1988) in which detailed discussions on this subject were presented. Our reasoning is based mainly on numerical experiments: the current method reproduced well the analytical Gaussian plume solution with a relatively small number of puffs.

### HOTMAC

HOTMAC is a three-dimensional mesoscale model based on a second-order turbulence-closure equations (Mellor and Yamada 1974, 1982). The governing equations are conservation equations for mass, momentum, potential temperature, mixing ratio of water vapor, turbulence kinetic energy, and turbulence length scale. The readers are referred to Yamada and Bunker (1988) for a detailed description of HOTMAC.

HOTMAC was used to produce wind and turbulence distributions, which provided RAPTAD with inputs such as *U*_{i}, *u*_{i}, and *σ*_{ui} in Eqs. (2) and (3). HOTMAC also provided potential temperature (or air density) distributions, which were used to compute *N*^{2} in Eq. (10) and *B* in Eq. (14).

According to Hanna et al. (1991), sampling sites were located relatively close to a source (mostly from 50 to 800 m) and experimental sites were flat. Thus, wind and turbulence distributions may be assumed to be homogeneous in the horizontal directions. Hanna et al. (1991) listed MDA files in appendix A of their report. From the MDA files input parameters used in HOTMAC simulations, such as date, roughness heights, wind speeds, and air temperatures, were estimated. Reported wind speeds were measured in the surface layer (less than 10 m above the ground). HOTMAC required wind speeds in the free atmosphere as boundary conditions to compute winds in the boundary layer. Wind speeds in the free atmosphere were estimated by using the observed surface wind speeds, roughness heights, and a logarithmic wind profile for a neutral stratification.

HOTMAC was run for 5 h to obtain quasi–steady state solutions with 20 × 19 × 16 (vertical) grid points. A large horizontal grid spacing of 10 km was used to reduce computer time. The magnitude of the horizontal grid spacing should not affect the results because only horizontally homogeneous solutions are of interest. Vertical grid spacing varied with height: 4 m for the first 20 m above the ground and increasing to over 500 m near the top of the computational domain (5000 m above the ground).

The horizontal extent of the HOTMAC computational domain was reduced artificially to 1.9 km × 1.8 km for RAPTAD simulations. Again, this resolution should not affect the wind and turbulence distributions because horizontal homogeneity is assumed. The vertical grid structure was unchanged.

## Verification of model results

Zapert et al. (1991) conducted an extensive evaluation of several air quality simulation models: DEGADIS, SLAB, AIRTOX, CHARM, FOCUS, SAFEMODE, and TRACE. These models were designed to simulate dosage/concentration distributions over a relatively short distance (<1 km) from a source over flat terrain. Input data vary for each model, but all use wind velocity at a single point.

The databases chosen for the model evaluation were the Desert Tortoise Liquid Ammonia Experiment in Nevada (Goldwire et al. 1985), the Burro Liquid Natural Gas (LNG) Spill Tests in California (Koopman et al. 1982), and the Goldfish Anhydrous Hydrofluoric Acid Spill Experiments (Blewitt et al. 1987). Three tests from each experiment were selected for evaluation. Model performance was evaluated by comparing observed and simulated concentrations by using measures such as bias, scatter, and correlation.

Hanna et al. (1993) conducted an even more extensive evaluation of 15 models (AFTOX, AIRTOX, Britter and McQuaid, CHARM, DEGADIS, FOCUS, GASTAR, GAUSSIAN PLUME, HEGADAS, HGSYSTEM, INPUFF, OB/DG, PHAST, SLAB, and TRACE) by using the observations from the eight field experiments (Burro, Coyote, Desert Tortoise, Goldfish, Hanford, Maplin Sands, Prairie Grass, and Thorney Island). As did Hanna et al. (1993), the following five groups of datasets are used here.

- Group 1: All continuous release, dense gas datasets, for short averaging times—that is, minimum time resolution in the data (Burro, Coyote, Desert Tortoise, Goldfish, Maplin Sands, and Thorney Island continuous).
- Group 2: Same as group 1 but for longer averaging times, approximately equal to the duration of release (several minutes).
- Group 3: All continuous release, neutral buoyancy, passive gas datasets (Prairie Grass and Hanford continuous).
- Group 4: All instantaneous release, dense gas datasets (Thorney Island instantaneous).
- Group 5: All instantaneous release, neutral buoyancy, passive gas datasets (Hanford instantaneous).

A total of 96 simulations were conducted. For each simulation, HOTMAC was run for 5 h to obtain quasi–steady state wind, potential temperature, and turbulence distributions. HOTMAC input parameters were specified according to the values reported in the MDA files, as mentioned earlier in section 2b.

If the release was in gas phase, the plume density was computed by using the molecular weight listed in the MDA files and the equation of state for ideal gas. For the releases that include aerosols, the plume density estimated by Spicer and Havens (1988) was used. Buoyancy [Eq. (14)] was initialized by using the plume density and the air density. The vertical velocity *W*

Sampling sites reported in the MDA were located mostly 50 to 800 m from the source over flat terrain. To detect a peak concentration, RAPTAD sampling sites were placed on the arcs whose radii correspond to the sampling distance reported in the MDA. Concentration averaging time for a peak concentration was varied from 1 to 600 s. RAPTAD simulation time varied from 4 to 30 min.

*C*

_{o}is an observed concentration and

*C*

_{p}is the corresponding predicted concentration.

For the neutral buoyancy, passive gas case (group 3 and group 5), the geometric mean biases were close to 1, and geometric variances were small. These results indicate that the present HOTMAC–RAPTAD models have a good predicting capability for the neutral buoyancy, passive gas case for both instaneous and continuous releases.

For the dense gas case (group 1, group 2, and group 4), MG values were within a factor-of-2 difference although the VG (=4.03) was relatively large for the continuous release, dense gas case with longer averaging time (group 2). The models overestimated the concentration values for the instantaneous release, dense gas case (group 4) and underestimated the values slightly for the continuous release, dense gas case (group 1). The bottom-right panel in Fig. 1 summarizes the models’ overall performance for all five groups. The solid line indicates the relationship between perfect model results and perfect observations. The area between the dashed lines indicates the values of MG that are within a factor-of-2 difference. All predicted concentration means are within a factor-of-2 difference of the observations, with small geometric variances (except for group 2, which has a moderate VG of 4.03). By comparing the results with those evaluated by Hanna et al. (1993), it is found that RAPTAD’s overall performance is at least as good as those of “better” models in their study. The results shown here and in Yamada (1999) indicate that the HOTMAC–RAPTAD has the capability of simulating neutral, buoyant, and dense gas plumes by using the single algorithm developed by van Dop.

Another way of evaluating model performance is through the use of residual plots where residual is defined as the ratio of a predicted concentration to the corresponding observed concentration. A perfect model should not exhibit any trend with variables such as wind speed and distance from the source and should not exhibit large deviation from unity (unity implies a perfect agreement between the modeled and observed values).

Figures 2 through 6 show residual box plots as functions of the distance downwind, wind speed, Pasquill–Gifford stability (PG class), and BINIT [initial buoyancy value, Eq. (14)] for the five groups of data discussed earlier. The horizontal lines for each box indicate the 2d, 16th, 50th, 84th, and 98th percentiles of the cumulative distribution function. The total number of points *N* used for statistical analyses is indicated directly below each box plot.

For the neutral cases, group 3 (Fig. 4) and group 5 (Fig. 6), the model residuals show very little trend with the distance downwind and wind speed and remain close to unity. The model residuals in group 3 (Fig. 4), however, show a dependency on PG class, with overestimates in unstable conditions and underestimates in stable conditions.

For the dense gas cases, group 1 (Fig. 2) and group 2 (Fig. 3), the model results exhibit a consistent performance with respect to the distance downwind, wind speed, and PG class. Most of the mean values are within the factor-of-2 boundary. For the instantaneous release, dense gas case, group 4 (Fig. 5), residual plots showed overestimates for distances close to the source, large wind speed, and slightly unstable stratification. Nevertheless, the overall performance of the model still is satisfactory, with MG = 0.57 and VG = 2.14.

## Summary

A total of 96 simulations were conducted to examine the performance of the nonneutral plume model algorithm developed by van Dop, which was incorporated into the HOTMAC–RAPTAD modeling system. The overall performance of the model in terms of geometric mean biases, geometric variances, and residual plots was found to be at least as good as those of the “better” models examined by Hanna et al. (1993).

Recently, Sykes and Cerasoli (1998) evaluated the performance of the Second-Order Closure Integrated Puff model (SCI PUFF; Sykes et al. 1993, 1996) with the MDA. They assumed a Gaussian form for the radial velocity distribution and adjusted empirical constants in their turbulent entrainment formula on the basis of comparison with observations. Their results indicated MG = 1.02 and VG = 1.30 for the passive releases, and MG = 1.10 and VG = 1.29 for the dense gas releases.

Although the results in this study were not as good as those in Sykes and Cerasoli (1998), our models can be used for both positive and negative buoyant plumes over complex terrain (Yamada 2000). In other words, the current modeling system is very versatile and is not tuned to any specific conditions.

## Acknowledgments

Ms. Ruey-in Chiou and Dr. Danny Lu assisted with the simulations. This work was supported by the U.S. Air Force Research Laboratory Airbase and Environmental Technology Division Contract F08637-95-C-6039 and the Small Business Innovation Research program. Dr. Steven Hanna and Mr. Joe Chang of Earth Tech Corporation kindly provided their report, computer software for statistical analysis, and archived data on diskettes. We greatly appreciate their help. Dr. Hanna is now with George Mason University, Virginia.

## REFERENCES

Blewitt, D. N., J. F. Yohn, R. P. Koopman, and T. C. Brown, 1987: Conduct of anhydrous hydrofluoric acid spill experiments.

*Int. Conf. on Vapor Cloud Modeling,*New York, NY, Amer. Inst. of Chem. Eng., 1–38.Gifford, F. A., 1982: Horizontal diffusion in the atmosphere: A Lagrangian-dynamical theory.

*Atmos. Environ.,***16,**505–515.Goldwire, H. C., T. G. McRae, G. W. Johnson, D. L. Hipple, R. P. Koopman, J. W. McClure, L. K. Morris, and R. T. Cederwall, 1985: Desert Tortoise series data report—1983 pressurized ammonia spills. Lawrence Livermore National Laboratory Report UCID-19075, Vol. 1, 243 pp. [Available from Lawrence Livermore National Laboratory, Publication Services, L-658, P.O. Box 808, Livermore, CA 94550.].

Hanna, S. R., D. G. Strimaitis, and J. C. Chang, 1991: Hazard response modeling uncertainty (a quantitative method). Vol. 2, Evaluation of commonly used hazardous gas dispersion models. Sigma Research Corporation for AFESC, Tyndall AFB, FL, and API, Report Nos. 4545, 4546, and 4547, 338 pp. [Available from API, 1220 L St. NW, Washington, DC 20005.].

Hanna, S. R., J. C. Chang, and D. G. Strimaitis, 1993: Hazardous gas model evaluation with field observations.

*Atmos. Environ.,***27A,**2265–2285.Kao, C.-Y. J., and T. Yamada, 1988: Use of the CAPTEX data for evaluation of a long-range transport numerical model with a four-dimensional data assimilation technique.

*Mon. Wea. Rev.,***116,**293–206.Koopman, R., and Coauthors, 1982: Burro series data report LLNL/NWC—1980 LNG Spill Tests. Lawrence Livermore National Laboratory report UCID-19075, Vol. 1, 286 pp. [Available from Lawrence Livermore National Laboratory, Publication Services, L-658, P.O. Box 808, Livermore, CA 94550.].

Legg, J., and M. F. Raupach, 1982: Markov-chain simulation of particle dispersion in inhomogeneous flows: The mean drift velocity induced by a gradient in Eulerian velocity variance.

*Bound.-Layer Meteor.,***24,**3–13.Mellor, G. L., 1974: A hierarchy of turbulence closure models for planetary boundary layers.

*J. Atmos. Sci.,***31,**1791–1806.Mellor, G. L., 1982: Development of a turbulence closure model for geophysical fluid problems.

*Rev. Geophys. Space Phys.,***20,**851–875.Pitchford, M., M. Green, H. Kuhns, and R. Farber, 1997: Characterization of regional transport and dispersion using Project MOHAVE tracer data.

*Proc. Visual Air Quality: Aerosols and Global Radiation Balance,*Bartlett, NH, Air and Waste Manage. Assoc., 181–200.Spicer, T. O., and J. Havens, 1988: Development of vapor dispersion models for nonneutrally buoyant gas mixtures—Analysis of TFI/NH3 test data. Chemical Engineering Department, University of Arkansas, Tech. Rep. ESL-TR-87-72, 115 pp. [Available from Engineering and Services Laboratory, Air Force Engineering and Service Center, Tyndall Air Force Base, FL 32403.].

Sykes, R. I., and C. P. Cerasoli, 1998: Incorporating dense gas effects into a Lagrangian puff model.

*10th Joint Conf. on the Applications of Air Pollution Meteorology with the A&WMA,*Phoenix, AZ, Amer. Meteor. Soc., 312–316.Sykes, R. I., S. F. Parker, D. S. Henn, and W. S. Lewellen, 1993: Numerical simulation of ANATEX tracer data using a turbulence closure model for long-range dispersion.

*J. Appl. Meteor.,***32,**929–947.——, D. S. Henn, S. F. Parker, and R. S. Gabruk, 1996: SCIPUFF—A generalized hazard dispersion model. Preprints,

*Ninth Joint Conf. on the Applications of Air Pollution Meteorology,*Atlanta, GA, Amer. Meteor. Soc., 184–188.Taylor, G. I., 1921: Diffusion by continuous movements.

*Proc. London Math. Soc.,***26A,**196–211.Van Dop, H., 1992: Buoyant plume rise in a Lagrangian framework.

*Atmos. Environ.,***26A,**1335–1346.Yamada, T., 2000: Numerical simulations of airflows and tracer transport in the southwestern United States.

*J. Appl. Meteor.,***39,**399–411.Yamada, T., and S. Bunker, 1988: Development of a nested grid, second-moment turbulence-closure model and an application to the 1982 ASCOT brush creek data simulation.

*J. Atmos. Sci.,***27,**562–578.Yamada, T., S. Bunker, and M. Moss, 1992: Numerical simulations of atmospheric transport and diffusion over coastal complex terrain.

*J. Appl. Meteor.,***31,**565–578.Zapert, J. G., R. J. Londergan, and H. Thistle, 1991: Evaluation of dense gas simulation models. EPA Contract No. 68-02-4399, Environmental Consultants, Inc., EPA-450/4-90-018, 112 pp. [Available from U.S. Environmental Protection Agency, Research Triangle Park, NC 27711.].