A major difficulty encountered in studying the Khamisiyah Pit release is the absence of any on-site detection of chemical warfare agent concentrations, combined with the absence of conventional meteorological observations near the release. Consequently, this analysis has completely relied on predictions by numerical weather models based on data-sparse background conditions. Moreover, dispersion models estimate the possible low-level health hazards due to the release. These models’ ability to simulate the dispersion patterns correctly varies due to uncertainties in both the models’ input data as well as their formulation. The input uncertainties exist in the characterization of the release and in the meteorological input data, consisting in this study of the predicted fields (winds, temperatures, mixing depths, heat flux, water vapor, clouds, precipitation, etc.) produced by mesoscale numerical weather prediction models.

To some extent, the Dugway and Edgewood tests and corroborative modeling have helped reduce the uncertainties in the proportion of total agent released and the relative importance of the release mechanisms (e.g., aerosolization, bulk spills in soil, and saturation of wood rocket crates) that dominated the pit release. However, uncertainties in the release characteristics still remain. Mesoscale meteorological fields also can contain errors contributing to uncertainties in simulating agent concentrations. Errors in the dynamic models can come from uncertainties in the initial fields and lateral boundary conditions, model numerics, and the need to parameterize sub-grid scale phenomena. These errors lead to uncertainties in agent concentrations predicted by transport and dispersion models. The transport and dispersion models’ formulation also provide sources of uncertainties in predicted concentration fields. For example, the atmospheric boundary layer is characterized by turbulence, which by definition cannot be precisely described. Uncertainties also occur in parameterizing various physical processes, such as deposition and decay. As a result, it is crucial to take all these uncertainties into account.

Recently, ensemble methods have been applied to quantify uncertainties in meteorological input data and these uncertainties’ impact on chemical warfare agent concentration predictions. An ensemble mesoscale forecast melds a large number of individual forecasts generated by perturbing one or more baseline simulations. Various approaches have been used to construct the ensemble members, including perturbing the initial global fields the mesoscale forecast uses (Molteni et al., 1996). These perturbations may be constructed from singular vectors, which determine the most unstable phase space directions in the early part of the forecast period. Straume et al. (1998) compared puff predictions from a dispersion model using such an ensemble forecast input and found a large spread in the predicted puff evolution. Hamill and Colucci (1998) mixed perturbation methodologies by combining the singular vectors (i.e., bred initial conditions) with independent objective analyses to initialize two different regional scale models. They used the corresponding results to test the skill of various probabilistic precipitation forecasts. Hanna et al. (1998) use the Monte Carlo method to estimate the uncertainties in photochemical model predictions of ground-level ozone concentrations due to the combined effects of uncertainties in the full set of input variables. Stensrud et al. (2000) describe ensembles created by two distinct approaches—perturbations of both initial conditions and model physics—to simulate mesoscale convective systems and show the latter approach may be just as effective in revealing inherent model variability.

To account for uncertainties due to model inputs and algorithms, the DoD and CIA modeling team followed the recommendations of the IDA panel (Section IV.A) and considered various meteorological and dispersion models and global data sources to assess the possible dispersion of chemical warfare agent from the Khamisiyah Pit detonation. This approach compared credible models that use different initializations and model parameterizations. Estimating the possible hazard based on several independent simulations using valid data and well-tested models produces: (1) meteorological data sets describing a reasonable, plausible range of weather conditions and (2) dispersion patterns yielding a set of credible estimates of concentrations and dosages from the release. Depending on the questions a given modeling effort is to address, the model runs’ results may be combined in various ways. Threshold contours may be superimposed to envelope regions unconfirmed by any other combination of meteorological and dispersion models.

Based on several criteria, the DoD team used a collection of atmospheric models run by their developers to assess the possible dispersion of chemical warfare agent from the Khamisiyah Pit detonation. The IDA panel recommended basic criteria for model selection, including using high-resolution mesoscale meteorological models and transport and dispersion models that accept temporally and spatially varying meteorological fields. The IDA panel also recommended DoD use models currently sponsored by various organizations under DoD and Department of Energy to perform additional modeling analyses.

Thus, the Naval Research Laboratory (NRL) teamed with the Naval Surface Warfare Center (NSWC) to link the COAMPS meteorological model and the VLSTRACK dispersion model. The Lawrence Livermore National Laboratory (LLNL) Airborne Release Advisory Capability (ARAC) operated the Mass-Adjusted Three-Dimensional Wind Field (MATHEW) diagnostic meteorological model linked with the Atmospheric Dispersion by Particle-in-cell (ADPIC) dispersion model. Finally, DTRA ran the OMEGA prognostic meteorological model linked to the HPAC/SCIPUFF dispersion model. In addition, responding to the IDA panel’s suggestion to include an established civilian mesoscale model to provide comparative results, the NRL obtained 48 hours of meteorological reconstruction generated by the MM5 mesoscale model from NCAR. Comparisons among MM5, COAMPS, and OMEGA indicated that these models produced similar reconstructions of the meteorology at Khamisiyah. The IDA Panel viewed LLNL’s MATHEW model as less capable because it modeled atmospheric phenomena with less fidelity, so it was not further considered.

As described below, the effort to model the Khamisiyah Pit release of chemical warfare agent consisted of several phases involving active modeling and subsequent external review. Using the above collection of model combinations throughout this effort helped bound the uncertainties in the model results.

The primary features of the mesoscale modeling systems used in the modeling are described briefly below.

From November 1996 to March 1997, NRL’s Marine Meteorology Division reanalyzed the meteorology surrounding the March 10, 1991, demolition at the Khamisiyah Pit. NRL used the NOGAPS forecast model for global and COAMPS for mesoscale reanalysis at the computing facility of the US Navy FNMOC with participation of FNMOC personnel. Some model highlights appear below. See Westphal et al. (1999) for more details.

a. Model Numerics and Physics. COAMPS is a mesoscale, primitive equation, non-hydrostatic model using a leapfrog-based, split time (e.g., Klemp and Wilhelmsen, 1978) integration scheme to predict dynamic and thermodynamic variables. Its grids use a terrain-following coordinate system to calculate flow over an irregular surface. The two-parameter, turbulence kinetic energy and mixing length scheme (Mellor and Yamada, 1974) is used to establish parameters for the subscale mixing to close turbulent flow equations. The mixing-layer height, a crucial input in dispersion models, is diagnosed based on the bulk Richardson number.

b. Observations Used. NRL prepared conventional and satellite observational data for the COAMPS reanalysis, collecting and storing data archived at the NRL, FNMOC, and other sources. These data included the operational set available in March 1991, the delayed observations, and the since declassified observations. NRL’s analysis used these types of observations:

The following are the primary differences between the operational data set available in March 1991 and the one used in the reanalysis. The reanalysis includes observations that arrived at FNMOC after the deadline for the operational forecast run. Conversely, a small number of secret observations were used operationally, but have not been declassified and were not used in the re-analysis. Finally, Special Sensor Microwave / Imager dataset marine wind speeds became operational in September 1990, but were not archived until May 1991, and therefore were not used in the reanalysis. Despite the absence of these two types of observations, the global reanalysis is superior to the operational reanalysis due to improved data assimilation and forecast models, as indicated by the current and historical analysis errors.

c. Convective Parameterization. The model includes large-scale and convective precipitation modules and incorporates evaporation of falling precipitation. A deep cumulus convection scheme is based on the Kuo (1974) scheme, which allows deep moist convection to occur when the moisture convergence exceeds a specified threshold. Cloud temperatures and mixing ratios follow a moist adiabat, thereby representing a non-entraining cloud. Explicit cloud microphysics is applied for grid-scale precipitation, including all the condensation on the fine grid and the large-scale precipitation on the outer grids. The fine grid uses the Kain and Fritsch (1993) convective parameterization for sub-grid-scale condensation. COAMPS cloud forecasts predicted substantial cloudiness over the area of Kuwait, southern Iraq, and northern Saudi Arabia over March 11–13, 1991. Despite the general trend of increasing cloudiness, the modeled accumulated precipitation exceeded 5 mm in only a few locations, with heavy amounts in the mountains north of the Tigris and Euphrates Rivers.

d. Grid Structure. For the Khamisiyah analyses, COAMPS used a Lambert conformal projection with standard parallels at 20� and 60� north latitude. NRL configured COAMPS to run with three nested grids described in Table A-2 and depicted in Figure A-16. The nested grid configuration uses an outermost domain corresponding to the boundary of a 45-km grid encompassing the entire Persian Gulf region, extending from eastern Africa to Afghanistan and the Black Sea to the Arabian Sea. The Khamisiyah site is enclosed by the two inner nested grids.

NRL chose a grid spacing of 5 km for the innermost grid in order to resolve local features (e.g., local terrain and land-water interfaces) around Khamisiyah and the high-resolution wind fields and temperature profiles in south central Iraq. The intermediate grid is a 15-km grid showing medium-resolution atmospheric conditions over the rest of the primary area. COAMPS uses a one-way, interactive, nested-grid treatment to propagate large-scale information from coarse to fine grids. The interface between grids is based on Perkey and Kreitzberg’s (1976) tendency modification method. Each grid uses a stretched vertical coordinate, with a vertical grid spacing of 10 meters at the surface to 6650 m near the top at approximately 31 km above ground level. The first 10 vertical levels located at roughly 10, 30, 55, 90, 140, 215, 330, 500, 750, and 1100 m above the ground had vertically averaged winds.

e. Surface Boundary Conditions. COAMPS calculates surface fluxes from similarity equations (Louis, 1979) and solves prognostic equations for ground wetness and temperature, where tendencies include the effect of the calculated surface fluxes. Both NOGAPS and COAMPS use the FNMOC global data sets of terrain height, surface albedo, climatological deep soil temperature, climatological ground wetness, and surface roughness, after interpolation to their respective model grids. These parameters’ spatial resolution is approximately one degree longitude and latitude; except for terrain, which originates from a master data set, the Defense Mapping Agency Digital Terrain Elevation Data Level 1 data, with a 100-meter resolution. The sea surface temperatures come from the daily archived sea surface temperatures fields from the Navy’s ocean analysis system OTIUS (Cummings et al., 1997), where the primary observations are multi-channel sea surface temperatures. The Naval Oceanographic Office produces multi-channel sea surface temperatures retrievals from high-resolution satellite infrared imagery, with additional buoy, ship, and bathymetry data.

g. Global Analysis. As described above, to provide three-dimensional dynamic and thermodynamic fields for the COAMPS initial and lateral boundary conditions, NRL performed a complete global reanalysis with the NOGAPS system, where data are assimilated every six hours with the operational configuration (0.75 degree horizontal resolution). The reanalysis used all meteorological data cited above, including the declassified field data from the reanalysis period. Data assimilation consists of quality control and MVOI, in which meteorological observations distributed in time and space are optimally combined with forecasts from previous analyses.

h. Data Assimilation. COAMPS applies intermittent dynamic data assimilation to correct model errors by incorporating observations into the solutions during the course of simulation. COAMPS assimilates observations every 12 hours, using its previously forecasted fields as the first guess, where observations within the � 6-hour window are assimilated into the first-guess fields. In other words, the model stops every 12 hours during the numerical integration, uses the model fields as a background to generate a new objective analysis, and restarts the model for the next integration period. The COAMPS quality control algorithms for observational data screen the data for integrity and representativeness, eliminate redundant observation, check the observations against climatological limits, check the vertical soundings for hydrostatic consistency, check the vertical wind shears or speed and direction, and check rawinsondes against the first-guess and neighboring observations. An MVOI analysis technique maps observations to model grids. Each restart incorporates fresh data to limit errors, but such data also causes spin-up errors due to imbalances in the re-initialized fields. NRL currently applies a digital filter to reduce temporal discontinuity in the surface winds at restart.

OMEGA constitutes a complete mesoscale dispersion modeling system, combining a prognostic mesoscale meteorological model with an embedded Lagrangian Atmospheric Dispersion Model. The DTRA used only OMEGA’s mesoscale meteorological component. Some model features are summarized below. Bacon et al. (2000) contains more details.

a. Model Numerics and Physics. OMEGA discretizes the conservative form of the Navier-Stokes equations. OMEGA solves its prognostic equations with a finite volume solver based on Smolarkiewicz’s positive definite advection algorithm (Smolarkiewicz, 1984) that naturally applies the finite-volume integral formulation to its unstructured grid. Time steps are adapted locally to the grid structure to satisfy a local Courant-Friedrichs-Lewy constraint, thereby increasing computational efficiency. The model does not separately resolve slow and fast mode behavior using distinct time steps. Applying an implicit solver to the acoustic terms attenuates vertically propagating waves.

b. Observations Used. Among the meteorological data OMEGA assimilated were WMO rawinsondes, NCAR pibal reports, from Saudi Arabia that were originally delayed, since declassified observations from USAF special operations and Navy ship rawinsondes and pibals.

c. Boundary Layer Parameterization. OMEGA accounts for the viscous sublayer and transition layer in its surface flux parameterization scheme. The 2.5-level closure technique (Mellor and Yamada, 1974) is used to define vertical turbulent exchange in the latter. OMEGA version 3.7 has adopted improved methods for flux parameterizations over land surfaces, including using closure functions for momentum and heat transfer coefficients that apply more generally to all stability conditions during diurnal cycles (Beljaars and Holtslag, 1991). A relatively simple surface layer model permits the use of an extensive set of look-up tables containing detailed non-uniform surface characteristics. OMEGA currently diagnoses the Planetary Boundary Layer (PBL) height based on the turbulence kinetic energy level; previous versions (2.0) used a prognostic equation.

d. Radiation and Convective Parameterization. OMEGA’s longwave and shortwave radiation schemes enable the cloud and precipitation fields to interact with each other, as well as with the surface. The model includes large-scale convective and precipitation modules and treats large-scale processes explicitly. A range of options exists to specify how to model the microphysics, where mixed-phase microphysics is allowed. Options for deep cumulus convection include parameterizations based on the Kuo (1974) and modified Arakawa-Schubert (1974) schemes.

e. Grid Structure. The OMEGA model grid adapts to underlying surface features (e.g., terrain, land-water interfaces) and, in its fully interactive mode, dynamically adapts to atmospheric phenomena as they evolve. In addition, grid resolution may be specified pointwise about a location of interest. The grid elements are made up of triangular prisms, with each prism’s axis aligned along the earth’s radius through its center (Figure A-17). The planar view of the grid therefore is a mesh of triangles, whose top and bottom faces can have any slope so that any terrain features can fit the grid. The bottom faces of the lowest row of cells describe the topography; hence, the bottom layer follows the terrain. The layers’ thicknesses can be adjusted so the layers relax to a horizontal position at the top of the domain. This grid structure allows surface terrain resolution down to less than 1 km.

f. Initial and Lateral Boundary Conditions. OMEGA used archived analysis produced in March 1991 by the NCEP Global Data Assimilation and Forecast System (GDAS, Kanamitsu, 1989) for the initial and lateral boundary conditions. The GDAS analysis products, generated in March 1991 at six-hour intervals, included horizontal wind components, geopotential heights, temperature, and humidity on a 180-km horizontal grid spacing with 12 vertical layers. At initialization, OMEGA implements a complete data assimilation module combining raw observational data with a dynamically consistent description of the atmospheric state. The data assimilation module includes an automated quality control procedure followed by a multivariate optimum interpolation scheme.

g. Surface Boundary Conditions. OMEGA uses eight types of surface characteristic databases: land and water distribution, terrain elevation, soil type, climatological vegetation index, land use and cover, and climatological sea surface temperature, subsurface temperature, and soil moisture. The first two databases are determined from global data sets and have been widely used in numerical modeling. The others rely on newer data sets from various sources. Table A-4 lists OMEGA’s available surface characteristics databases.

h. Data Assimilation and Global Analysis. As described for the initial and lateral boundary conditions, OMEGA used archived analyses from the NCEP GDAS for gridded first-guess data. In 1997, OMEGA was used in a pure forecast mode; but in 2000, it was used in dynamic data assimilation mode.

Some of the model’s features are highlighted below. See Warner and Sheu (2000) for more details.

a. Model Numerics and Physics. MM5 is a regional-scale primitive equation model that can be configured in the hydrostatic or non-hydrostatic mode. The Khamisiyah analyses used the non-hydrostatic mode. MM5 uses a terrain-following coordinate in pressure and solves its finite difference equations in the horizontal with a time-split scheme using the explicit leapfrog algorithm. A semi-implicit vertical differencing scheme allows a short time step independent of the model’s vertical resolution; thus, the time step depends only on the horizontal grid interval. Its boundary layer physics package can be a scheme based on either the revised Blackadar (1979) PBL model or the Troen and Mahrt (1986) PBL model as implemented in the Medium-Range Forecast (MRF) model.

b. Observations Used. NCAR used NCEP’s Automatic Data Processing system global upper-air and surface observations obtained operationally from the standard WMO stations, with additional surface observations obtained from other sources. From its own research archive and the FNMOC archive, NRL provided NCAR with raw data (surface observations, rawinsondes, and a few cloud track wind observations). These data sets, which included declassified and delayed data described in Appendix A, Section II.B above, provided roughly five additional soundings for use in the MM5 reanalysis. Figure A-19 shows the locations of the surface and radiosonde stations (a combination of NCEP Automatic Data Processing and NRL observations) at 1200 UTC March 10, 1991.

c. Radiation and Convection Parameterization. Both longwave and shortwave atmospheric radiation schemes are included. The model includes large-scale and convective precipitation modules. The large-scale processes are treated explicitly. Options for deep cumulus convection include parameterizations based on the Kuo (1974) and modified Arakawa-Schubert (1974) schemes.

d. Grid Structure. NCAR used the Lambert conformal map projection with standard parallels at 30� and 60� N. The MM5 grid has moving, overlapping nesting capability with one- or two-way interaction. The nesting ratio is predefined at 3:1. Table A-5 shows the nested grid characteristics for the Khamisiyah analysis. As with COAMPS, the MM5 horizontal grid configuration includes an outermost grid encompassing two finer grids.

e. Initial and Lateral Boundary Conditions. NCAR compared the sensitivity of the mesoscale analysis to the global first-guess fields by performing several different MM5 simulations using both the NCEP global reanalyses (the NCEP Aviation model in 1997 and the NCEP MRF model in 1999) and the ECMWF reanalysis fields, both available on a 2.5� x 2.5 degree grid. Linear temporal interpolations between 12-hourly NCEP or ECMWF global reanalyses define the lateral boundary conditions for the outer grid, Grid 1. To apply lateral boundary conditions, MM5 implements a relaxation technique that relaxes or nudges the model’s predicted variables toward a large-scale analysis, a synthesis of the first-guess field, and observations.

f. Data Assimilation. As described above, FDDA techniques provide a method to assimilate current observations into the meteorological model results over an extended analysis period. MM5 applies Newtonian relaxation, which gradually drives the model variables toward a gridded analysis by including an extra forcing term in each governing equation. For the MM5 Khamisiyah runs, this method nudged the model toward a synthesis of the first-guess NCEP and ECMWF global fields and synoptic data that are updated every 12 hours using a Cressman-type method. The objective analysis is performed at both the mandatory and non-standard levels. This method then interpolates output defined at pressure levels to sigma levels for assimilation into the model fields. Data are assimilated on the coarse grid only.

g. Surface Boundary Conditions. The surface energy budget is treated according to Zhang and Anthes (1982). A major enhancement in MM5 is a new, high-resolution analysis of the surface characteristics that strongly influences the boundary-layer winds and stability through the heating and roughness effects.

OMEGA, COAMPS, and MM5 have much in common: all are three-dimensional, primitive-equation, mesoscale models solving the non-hydrostatic, compressible form of the dynamic equations and use many of the same parameterizations of physical processes (e.g., surface fluxes and moist convection).

However, these models have different features. For example, COAMPS and OMEGA are used in an operational setting, so operational constraints balance features related to data input/output considerations and objectives such as physical fidelity and numerical accuracy. As an example, COAMPS and OMEGA process observational data and perform quality control in a fully automated fashion. Conversely, MM5 is mostly used in research applications and thus contains numerous optional physical algorithms.

MM5 and COAMPS are the established mesoscale meteorological models. MM5 is widely used in research communities, COAMPS is the operational prediction model for the Navy and DoD. The basic equations of both models are based on Klemp and Wilhelmson (1978). Both models use a staggered grid both horizontally and vertically. Grid nesting efficiently treats a wide range of temporal and spatial scales. On the other hand, the OMEGA grid is unstructured horizontally and adapts to both underlying surface features and dynamically evolving atmospheric phenomena. This approach achieves local accuracy of the numerical solution with a single, non-uniform grid and does not require communication between separate nesting grids.

The PBL is a critical factor in controlling mesoscale weather systems. Because of the large fluxes of heat, moisture, and momentum near the earth’s surface, there is generally an agreement on the need for high-resolution treatment of the physics of this layer (Blackadar, 1979). However, the three models apply different approaches to modeling the PBL. COAMPS and OMEGA apply a fine vertical resolution to resolve the PBL, including the stable boundary layer. In addition, they apply the level 2.5 PBL model developed by Mellor and Yamada (1974). The crucial phenomenon to resolve is the transport of mass and momentum in the PBL by large energetic eddies. Traditional local-gradient methods cannot adequately treat such a well-mixed atmosphere. Mellor and Yamada’s (1974) higher-order closure methods, though computationally expensive, are capable of representing a well-mixed boundary layer. On the other hand, the lowest MM5 model computation level is approximately 40 m above ground level, with increasing layer depths above, so it is difficult for the model to properly resolve the shallow nocturnal PBL. Local-gradient theory may fail because it does not account for the influence of large eddy transports and does not treat entrainment effects. MM5 uses non-local atmospheric boundary layer schemes that are more effective for coarser grids.

The PBL’s spatial variability can result from a range of mechanisms, including topographic elevation variation, land-and-sea breeze circulation, and local contrasts in physical properties at the desert surface. The OMEGA and MM5 2000 reanalyses have used high-resolution surface databases. The March 97 COAMPS reanalysis used 1 degree resolution data; whereas the COAMPS 2000 and the two month continuous reanalysis were based on 1 km database.

Since the model simulations’ objective was to best analyze the area’s meteorological conditions, use of FDDA or hindcasting was crucial. Although grid structures, numerical solvers, and PBL parameterizations all contribute to different model features, the most significant difference among the three mesoscale models is probably in their data assimilation strategies, described below.