Bret A. Schichtel and Rudolf B. Husar

Center for Air Pollution and Trend Analysis (CAPITA)

Campus Box 1124

Washington University

St. Louis, MO 63130

September 25, 1995


A Monte Carlo model for the simulation of regional scale transport, transformation, and dry and wet removal is presented. The model was newly re-designed in a modular framework, separating the emissions, transport and kinetics calculations. The transport module employs a quantized Monte Carlo technique for the simulation of atmospheric boundary layer physics. Kinetic processes are simulated using rate equations where the rate coefficients are dependent upon meteorological variables, and thus vary in space and time. The rate coefficient equations are determined via a tuning process in which simulated values are compared to observed measurements. Results from simulations of SO42- and SO42- over the eastern US are presented for 1992. Comparisons of simulated daily SO42- concentrations to observations had r2 = 0.35 - 0.83 depending on season and location, while for weekly SO42- wet deposition rates r2 = 0.5 - 0.95. The utility of the model to investigate the source receptor relationship is demonstrated by explicitly examining the role of emission rates, transport, and kinetic processes in the attribution of sulfur dioxide and sulfate at a receptor in Massachusetts during the summer of 1992.



During the 1970's, air quality measurements demonstrated that pollutant emissions could be transported over large distances, significantly influencing distant receptor concentrations. The realization of regional scale transport of pollutants led to numerous air quality networks, models, and large scale field studies to relate the contributions from near as well as distant sources to receptor concentrations, i.e. to quantify the source receptor relationship. Many of the initial studies focused on sulfur oxides, and were designed to attain better understandings of the physical/chemical processes governing regional scale transport. For example, in Europe, the OECD established the Long Range Transport of Air Pollutants (LRTAP) program in 1972 (1), and in the US, the Sulfate Regional Experiment (SURE) was conducted (2), along with numerous plume studies.

Associated with these studies were the development of models. The initial models tended to be simple, using aggregated and parameterized meteorology and rate coefficients to simulate the physical/chemical processes responsible for regional transport. These included trajectory models (3, 4), statistical models (5, 6), as well as simple diffusion and box models (7, 8). As the understanding of the atmospheric processes increased, these types of models were generalized to include processes such as non-linear chemistry and atmospheric diffusion (9, 10).

During the 1980's, with the concern over acid rain, substantial effort was devoted to the development of elaborate Eulerian models, such as the Regional Acid Deposition Model (RADM) (11) and the Acid Deposition and Oxidant Model (ADOM) (12). These models attempt to simulate atmospheric processes, such as in-cloud oxidation, with great detail. However, the data and computer resource needs have restricted their use to research activities with runs limited in time and space. Also, it is difficult to separate the contributions from the various physical/chemical processes involved in the source receptor relationship (13). As a result these models tend to be operated in a "black box" approach.

In this paper, an IBM-PC based Monte Carlo model for the simulation of regional scale transport, transformation, and removal of atmospheric pollutants is presented. This model is an extension of the CAPITA Monte Carlo model developed in the 1980's (9, 14). The model is "data driven," and was designed to take advantage of multiple meteorological, emission, concentration, and deposition databases that are available from various data suppliers. Beyond the calculation of ambient concentrations and deposition fluxes, the model was designed to be used as a diagnostic tool for the investigation and quantification of the physical/chemical processes governing the source receptor relationship, as well as for the estimation of unknown emission fields given ambient concentration data.

The paper presents the simulation methodology, along with the underlying physical and chemical processes. The suitability of the model for regional scale simulation is shown through the simulation of sulfate and the total sulfate wet deposition over the Eastern US. Also, the utility of the model to investigate the source receptor relation is demonstrated by examining the role of emission rates, transport, and kinetic processes in the attribution of sulfur dioxide and sulfate at a receptor in Massachusetts.
Back to Contents


Monte Carlo Approach

In the Monte Carlo approach to air pollution simulation, individual quanta or particles are subjected to the same physical and chemical processes and events in the computer as in the physical world. It can be thought of as a direct simulation of atmospheric pollutants. This approach contrasts with more customary numerical simulations in that it does not require explicit conservation equations and their solution as the basis for the simulation.

The Monte Carlo approach to the simulation of atmospheric pollutants has the following key characteristics. First, each emitted quantum contains a fixed quantity of mass of various pollutants based on the source's emission rate. Individual quanta are then subjected to transport, transformation, and removal processes. Mass conservation is maintained at the quantum level by accounting for the mass that has been chemically transformed and physically removed as an assembly of quanta travel through three dimensional space (Figure 1). An airmass is depicted as an ensemble of multiple quanta. The dispersion is represented by the spread of the quanta arising from the integration of the equations of motion, where horizontal and vertical mixing are based on mean and random wind components. This approach can be considered to be a Lagrangian simulation since the trajectories and the fate of individual quanta are followed. In the limit of infinite quanta, this approach results in a stochastic solution of the dispersion equation, yielding results equivalent to an Eulerian solution to the governing equations.

 Figure 1. A depiction of the Monte Carlo approach simulating the physical and chemical processes applied to a quantum.

In the Monte Carlo approach, the simulation of the meteorological dispersion of pollutants proceeds independently from the nature and chemistry of the pollutants. This is based upon the assumption that meteorology is not significantly influenced by trace constituents in the atmosphere; however, it is understood that kinetic processes are dependent upon the pollutant transport. Consequently, the calculation of the dispersion of the pollutants can proceed independently from the calculation of kinetic processes. This allows for the separate examination of the influence of the transport and kinetics processes upon receptor concentrations.

Monte Carlo simulations of atmospheric processes were developed in the 1970's (15, 16). However, this technique lost favor in the 1980's as other approaches were pursued. Due to the multiple benefits of this modeling approach, it regained popularity during the 1990's for mescoscale simulations (17, 18, 19), and is believed to offer considerable improvements over conventional Gaussain-plume models (20). The use of the Monte Carlo technique for regional scale simulations has received less attention with little development since the initial studies were conducted by Patterson et al. (14) and Husar (9).

Back to Contents

Data Flow

The independence of transport from chemical processes allowed for the separate development of a meteorological module, for the calculation of pollutant dispersions, and a kinetic module, for the calculation of pollutant concentrations and deposition rates. The flow of raw meteorological, emission, chemical, etc. data through the meteorological and kinetic modules is presented in Figure 2. The meteorological module transforms meteorological data into a Lagrangian airmass history database, by computing the transport of conservative species through three dimensional space. The first step in the processing of the meteorological data, involves converting raw data from a random spatial distribution to a grid. Various institutions convert meteorological state variables collected by the National Weather Service to a gridded format using primitive equation models. Currently, we use data derived from the National Meteorological Centers Nested Grid Model, NGM, (21, 22) and Colorado State University's Regional Atmospheric Modeling System, RAMS (23). The meteorological database is then converted into Lagrangian airmass histories by the Eulerian/Lagrangian transformer. Along each particle's calculated trajectory, the particle's positions and meteorological variables, such as temperature, humidity, precipitation, etc., are saved. The resulting Lagrangian database, can be either source oriented, simulating emissions forward in time from sources to receptors, or receptor oriented, simulating the transport of airmasses backwards in time.

 Figure 2. The Monte Carlo Model data flow chart.

The kinetic module combines a source oriented Lagrangian airmass history database with other data, such as emissions and surface meteorological data, for the calculation of the kinetic processes. This is accomplished by first assigning species mass weights to the quantum when it is released from the source. The weighting is dependent on each source's emission rate. Kinetic rate equations are then integrated along each quantum's trajectory, accounting for the transformation and removal of the pollutants with time. This creates additional Lagrangian variables containing species mass weights. Receptor concentration fields can then be created from the resulting database by defining a receptor volume and summing the species mass for all quanta located in the volume. Deposition fields can also be created in a similar manner, but only a receptor area needs to be defined. The concentration and deposition fields can be created during the processing in the kinetic module, or as a post processor on the resulting kinetic Lagrangian variables.

Normally, during the Eulerian to Lagrangian transformation only three dimensional variables like humidity are stored along with the particles position. Within the kinetic module, quanta airmass histories can be enriched with two dimensional Eulerian variables, such as precipitation and cloud cover, required for the kinetic calculations. The data used to enrich the airmass histories can be located in multiple databases. However, each database must use the same Eulerian grid and time step. This is a flexible and efficient method allowing for the inclusion or removal of meteorological variables and data sources for the calculation of kinetic processes.

Back to Contents

Transport and Kinetic Processes for Regional Scale Simulations

In a simulation of atmospheric transmission of a pollutant from its source to receptor, the temporal and spatial scale of the simulation must match the atmospheric residence time and transport distance of the pollutant. In this paper, we are interested in synoptic/regional scale pollutants. Consequently, the scales of interest are on the order of 5000 - 10000 km in space, and between 1 and 10 days in time. However, the understanding of the regional scale phenomena lies in the proper simulation of the smaller scale processes that comprise the synoptic picture.

Regional Scale Transport

Regional scale dispersion is the result of the coupling of advective and random processes in the atmosphere. The random processes are caused by the stochastic nature of turbulence in the atmospheric boundary layer (ABL). The ABL experiences diurnal (Figure 3) and seasonal cycles driven by the solar radiation. During the daytime, the ABL is unstable and well mixed. This is a result of strong heating of the surface from solar radiation which produces temperature gradients in the lower layers of the atmosphere. The temperature gradients produce discrete convective elements called thermals. These are buoyant parcels of air 50-1000 m in diameter which develop above a surface that is considerably warmer than the overlying air. The vertical turbulent structure of the ABL then becomes organized in a pattern of updraughts and downdraughts. (24).

 Figure 3. The dynamics of the planetary boundary layer, and pollutants being emitting into this layer as simulated by the CAPITA Monte Carlo Model. The mixing heights were taken from the NGM Eulerian database (21, 22), averaged over the Southwest, during the summer of 1992.

This mixing process is efficient enough that generally the potential temperature, relative humidity, and pollutant concentrations of the thermals are conserved (25). Pollutants emitted within the ABL quickly lose their identity and become part of the well mixed layer, extending up to heights of several kilometers. The time scale associated with the vertical mixing ranges from 30 minutes to several hours during afternoon conditions (26, 27).

In the late afternoon, the surface heating is reduced and the ground temperature becomes nearly equal to the overlying atmosphere. At this point, there is no upward heat flux to sustain the convective thermals and the boundary layer collapses, resulting in the bulk of the former daytime boundary layer becoming stable and stratified. The lack of vertical exchange confines the movement of air parcels to shallow layers roughly parallel to the ground surface. The pollutants that were part of the daytime mixing layer are now transported in horizontal layers with varying wind directions and speeds. Unlike the daytime conditions, the nocturnal transport within the former mixed layer will exhibit substantial veer and shear, the top layers heading in directions and at speeds that may be markedly different from the movement of the bottom layer (28, 29). It is this repeated process of vertical mixing followed by vertical wind veer which is essentially responsible for regional dispersion.

During the evening hours, the vertical extent of the nocturnal mixing height is only on the order of hundreds of meters. An important feature of the nocturnal surface based mixing layer is that it generally coexists with a cooled surface. Thus, the mixing process has to overcome the damping imposed by a stable stratified surface layer. These processes usually create a sharp impermeable interface between the former mixed layer and the surface based nocturnal mixing layer.

The interaction of a plume with the mixing layer is presented in Figure 3. As shown, emissions from a source in the nocturnal mixing layer are trapped in that shallow layer, resulting in high concentrations for a given emission rate. The pollutants confined to the surface layer are mixed to and exposed to the absorbing ground surfaces and to the chemical mix of that surface layer. Thus, dry deposition and potential chemical reactions significantly influence the fate of these pollutants. Emissions from a plume emitting above the mixing height or those left from the previous day's mixing are confined to a relatively shallow layer, de-coupled from the ground. Matter remaining above the surface layer overnight is not depleted by dry deposition, and chemical reactions occur with species mixed during the previous day and those emitted directly into the stratified atmosphere. During the morning hours, the mixing height begins to grow due to solar heating, diluting the pollutant concentrations within the mixing layer and entraining pollutants that were above the nocturnal layer.

Transport Simulation. The simulation of the regional transport in the Monte Carlo models is conducted by moving inertialess particles in the Eulerian frame according to the equation:


where x (=x,y,z) represents the particle position vector, and (= ) and u' (= u',v',w'), represent the time-averaged and fluctuating components of the flow field respectively. The application of equation (1) results in the transport and dispersion of the particles. The basic meteorological information required for this consists of the mean wind field and the fluctuating components u',v',w' at all times over the whole model domain.

The mean wind field is generated from meteorological models, such as the National Meteorological Center's synoptic scale model (21, 22) and Colorado State University's Regional Atmospheric Modeling System (RAMS) (23). The random component of the velocity field occurs at a resolution smaller than the meteorological model grids used to generate wind fields. Consequently, the fluctuating components are calculated from the available meteorological variables and treated as a sub-grid phenomena.

Vertical Diffusion. The above discussion, identified two key inputs for the simulation of vertical diffusion within the planetary boundary layer. First, an upper boundary or "lid" exists up to which convective mixing occurs, and second, the vertical daytime mixing occurs rather rapidly. In Monte Carlo simulations, the rapid daytime vertical mixing is implemented by changing the height of a particle at each 2-hour time step. Due to the rapid mixing, it is assumed that the particle height at each time step is independent of the height during the previous step, and is chosen according to a uniform probability distribution within the mixing layer. The exception is that material which has been previously transported above the local mixing height is subjected only to the mean wind vector.

A graphical illustration of the implementation of the above processes is presented by the three trajectories in Figure 3. The intense vertical mixing within the mixing layer is indicated by the crossing of the trajectory lines. The stable layers, during nocturnal hours, are illustrated in the lines of near constant height for each quantum.

Horizontal Diffusion. A common approach for the calculation of the horizontal diffusion is based upon the Prandtl mixing length model. This allows for the parameterization of the diffusion by an effective eddy diffusion coefficient K. The diffusion process is implemented as a random walk displacement of radius for each time step in the model. The value of the diffusion coefficient is a geographical function of the time of day and season. Values of K ranging between 102 m2/s in stable conditions and 106 m2/s during intense convective activity have been used in the past (14, 30). It is recognized that using an eddy diffusion coefficient to represent horizontal diffusion is a crude approximation. However, it will be shown below that horizontal diffusion has a small effect on regional dispersion.

Examples of Monte Carlo Simulations. The dispersion of a puff of emissions from a source in western Tennessee is illustrated in Figure 4. This figure was created by simulating the release of pollutants from a source using 25 particles. The position of the particles were plotted every two hours after their initial release at twelve in the afternoon and tracked for three days. The left hand picture contains the trajectories and height history of the particles when no horizontal diffusion was imposed, i.e. K = 0 m2/s. The right hand figure contains the trajectories and height history when a diffusion coefficient of K = 104 m2/s was employed for both day and night.

 Figure 4. Illustration of transport in the Monte Carlo simulation. The left figure displays the dispersion of a puff over three days of travel without horizontal eddy diffusion and the height history of each particle in the puff. The right figure is a similar puff, but an eddy diffusion coefficient of 10,000 m2/sec was used.

The initial six hours of transport from the source took place within the mixing layer. The wind velocity was relatively constant throughout the mixing layer, so the horizontal spreading of the plume was due to the turbulence. After approximately six hours of transport, the puff was in southern Illinois. At this point, the mixing layer had collapsed and the trajectories were transported in separate de-coupled layers with varying wind velocities producing horizontal spreading. This is seen in the left hand part of Figure 4 where the lateral spreading of the puff was due solely to wind veer. After approximately one day of travel, the airmass parcels were over the northwestern part of Ohio. At this point, it divided into a surface airmass heading towards New Jersey and Connecticut, and an upper air airmass heading towards Maine. The bifurcated plume continued to spread laterally, and after three days of transport the puff had a width of over 600 km encompassing most of New England and a height over three kilometers. It is evident that the large lateral spreading of the plume after three days of travel was primarily due to the wind veer rather than the horizontal eddy diffusion.

The simulation of airmass transport in the lowest few kilometers of the atmosphere is presented in Figure 5. This figure displays the position of particles on 01/16/1992 noon Greenwich mean time (GMT) which were released at some previous time from multiple sources uniformly distributed over the western US. In this sense, it is a flow visualization "experiment." The Eulerian to Lagrangian transformation was conducted using meteorological data generated by RAMS. A nested grid was employed where a grid of 60 km was used over the western US and a 24 km grid over Southern California and the Colorado Plateau (31).

A prominent characteristic of the western US is the complex terrain consisting of multiple large valleys and mountain peaks. Also, during January there is little solar input resulting in maximum daily mixing heights of approximately 600 m, as estimated from the RAMS data. The combination of complex terrain and poor vertical mixing resulted in the uneven spatial concentrations of particles displayed in Figure 5. The clumping of particles is due to the constant emission of particles into poorly ventilated areas. This is particularly noticeable in the Southwest where there are accumulations of particles in the San Joaquin Valley, Death Valley, and the Colorado River Valley. The poor vertical mixing also forces particles to travel around features, such as the Sierra Mountains. During the summer season, there is significantly more vertical mixing and horizontal dispersion resulting in more uniform particle patterns over the western US.

Back to Contents

Chemical Transformation and Removal Kinetics

In the Monte Carlo model, a semi-empirical approach is employed for the simulation of pollutant transformation and dry and wet deposition. The approach is deterministic in that rate equations are used for the calculation of kinetic processes applied to each quantum. At a quantum's release, it is assigned a species mass weight based upon the source's emission rate. Mass conservation is maintained at the quantum level by accounting for the mass that has been transformed and removed. The kinetic rate coefficients are dependent upon the chemical, meteorological, and geological environment of the quanta. Consequently, these coefficients vary with space and time creating variable transformation and removal rates.

The simulation of the kinetic processes is empirical, in the sense that the rate coefficient relationships are determined through a fitting process matching simulated to observed concentration and deposition fields. The main thrust of this approach is to make maximum use of existing meteorological, emission, receptor, etc. data in conjunction with a physical-chemical model to derive the best set of rate coefficient equations. This approach assumes that the correct pollutant transport and emission rates are known. It is hypothesized that given sufficient quality and quantity of data, it would be possible to identify a unique set of rate parameters, thereby simulating the physical/chemical processes.

SO2 - SO42- Simulation Over the Eastern US. Recently, a project was undertaken for the retuning of the Monte Carlo model for the simulation of the SO2 - SO42- system over the Eastern US during 1992. The following will present the resulting rate coefficient equations and some of the resultant simulated SO42- concentrations and wet deposition rates. The process by which the rate coefficient relationships were developed and their justification is given in Schichtel (32).

The rate of change of a quantum's SO2 mass can be expressed as:


where kt4, kd2, and kw2 are the oxidation, dry deposition, and wet deposition rate coefficients of SO2 respectively. Similarly, the rate of change of the SO42 mass can be expressed as:


where kd4 and kw4 are the dry and wet deposition rate coefficients of SO42- respectively. It was assumed that all SO2 and SO42- mass arose from the simulated sources only, thus, background SO2 and SO42- were neglected.

The tuning process was conducted using the meteorological data generated from the National Weather Service's Nested Grid Model (21, 22). The data were available on a stereographic grid with approximately 180 km size. Using the NGM data, a Lagrangian database was created by releasing particles from the center of each meteorological grid cell, east of Colorado. Approximately 350 sources were simulated in this manner. The three dimensional location of each particle, along with the specific humidity and temperature, were deposited in the Lagrangian database every two hours. Each particle was tracked for seven days, or until it was transported past the edge of the meteorological grid.

The NGM data contained a precipitation field. However, it was found that the precipitation field did not compare well with measured data. Deviations in the seasonal precipitation rates up to a factor of two were found depending on the season and location (32). The NGM precipitation field was replaced by hourly measured precipitation data from the US Control COOP Hourly Precipitation database (33). This was accomplished by averaging the precipitation rates over a two hour period, and spatially interpolating the data to the same grid used by the NGM meteorological data. The resulting database was then suitable for input into the kinetic module. A similar process was used for the creation of gridded total cloud cover data from hourly data in the National Solar Radiation Database (34).

Emissions were obtained from the 1985 NAPAP inventory (35). The original data were provided as point and area emissions on a 20 km grid. The data were transformed to the meteorological grid by summing all point and area emissions together for all emission grid cells which fell into the same meteorological grid cell. No attempt had been made to account for variations in 1992 emission from those in 1985.

The observed data utilized for the tuning process consisted of ambient fine sulfur aerosol concentrations from the IMPROVE (36) and the NESCAUM (37) monitoring networks, airport visibility data (38), ambient SO2 concentrations from the AIRS database (39), and total SO42- wet deposition data from the National Atmospheric Deposition Program/National Trends Network (40).

The final relations derived for the rate coefficients are:

kt2= 0.001 + SR (0.02 + 0.05 SH/SHmax + 0.1*PF) (4)

SR = (1-.65TS)GSR (5)

kd2 = Vd / H (6)

Vd = .22 + 2 * SR (7)

kd4 = /10 (8)

kw2 = Ws2/H * P = Sqr(40000 / (SO2CC + 2)) * P (9)

kw4 = Ws4/H * P = 250 * P (10)

where SR = ground level solar radiation [kW/m2]
TS = total cloud cover
GSR = ground level solar radiation in the absence of clouds [kW/m2]
SH = specific humidity [kg Water/kg Air]
PF = precipitation flag, i.e. PF = 1 if it rains PF = 0 otherwise
Vd = deposition velocity [cm/s]
H = mixing height [m]
Ws2, Ws4 = washout ratios for SO2 and SO42- respectively
SO2CC = SO2 column concentration [mg/m2]
P = precipitation rate [m/hr]

To illustrate how the meteorological variables affect the kinetic rate coefficients, a quantum was released from St. Louis on July 15 at 10:00 GMT, and tracked for four days. Figure 6 presents the quantum's trajectory, and six of the airmass history variables necessary for the calculation of the rate coefficients. The resulting rate coefficients and a corresponding sulfur budget are presented in Figure 7. As shown, the quantum was released, and remained, below the mixing height for the initial fourteen hours. Thus, it was exposed to the ground and a SO2 dry depositions rate of approximately 5%/hr occurred. This resulted in about 38% of the SO2 to be dry deposited. During this time, the solar input cycled from 0 to 0.6 back to 0 kW/m2 producing transformation rates from 0.1 %/hr to 3 %/hr. This caused approximately 12% of the SO2 to be transformed to SO42-. With the end of the solar input, the mixing layer collapsed and the quanta was above the mixing height. These conditions resulted in no deposition of the SO2 and little transformation. During the beginning of the second day, the ground level solar flux increased and the quantum encountered a light precipitating airmass. The precipitation caused total removal rates of the SO42 to be as high as 12%/hr removing 60% of the remaining SO42- mass over the next 12 hours. The precipitation also washed out a significant fraction of the SO42-. However, the enhanced transformation rates during precipitation, caused the overall SO42- mass to increase to 20% of the total sulfur. After this period the precipitation increased, and after two days of travel, nearly all of the SO2 and SO42- mass of the quantum had been dry and wet deposited.

 Figure 6. The variation of six airmass history variables for a quantum that was released from St. Louis on July 15, 1992 at 10:00 GMT.

 Figure 7. The variation of the kinetic rate coefficients along a quantum's trajectory and the corresponding quantum's sulfur budget. The quantum was released from St. Louis on July 15, 1992 at 10:00 GMT.

Simulated SO42- concentrations. Using the above rate equations, the Monte Carlo model was run for all of 1992, generating 24-hour concentration and deposition fields over the Eastern US. The surface concentrations were calculated for receptor volumes defined by the meteorological grid cell and the mixing height. Figure 8, presents the average simulated SO42- concentration fields during the periods January - March, quarter 1, and July -September, quarter 3. Also, plotted for comparison are "observed" SO42- isopleth maps generated from the fusion of IMPROVE and NESCAUM SO42- data and airport visibility data (41). During quarter 1, the observed SO42- is nearly uniform over the east with a concentration of about 3 mg/m3. The simulated concentrations have more spatial variability then the observed, but reproduce the general spatial pattern and magnitudes. The simulation deviates from the observed the greatest in a belt from eastern Ohio to the Chesapeake Bay, with simulated concentrations greater than 5 mg/m3 compared to about 3 mg/m3 for the observed SO42-.

The simulation during quarter 3, also reproduces the general spatial pattern and magnitudes of the observations. The highest observed concentrations are from southern Pennsylvania and north Virginia to western Kentucky and Tennessee ranging from 6 to greater than 7 mg/m3. The simulated SO42- has a similar high concentration band, however, it is pushed slightly to the north and west of the observed SO42- band. The concentrations for both the observed and simulated SO42- progressively decrease in all direction from the high concentration band producing similar concentrations over New England (~ 3 mg/m3), the South (~5 mg/m3), and into the Mid west (2 - 5 mg/m3).

 Figure 8. The simulated SO42- concentration fields during quarters 1 and 3 1992, and the corresponding "observed" SO42-isopleth maps generated from the fusion of IMPROVE and NESCAUM fine sulfur measurements and airport visibility data.

The daily averaged simulated SO42- are seasonally compared to observations in correlation plots (Figure 10). Due to the coarse meteorological grid used, it is unreasonable to expect the model to simulate concentrations at a single receptor. Therefore, these comparisons were conducted by averaging the model results and the NESCAUM/IMPROVE observations over a New England and East Central spatial domains as defined in Figure 9. The simulated SO42- in New England compares well with the observations, with squared correlation coefficients varying between 0.56 during quarter 3 and 0.83 during quarter 2. The data points are centered around the one to one line, except for quarter 1 where the model simulation tends to underestimate the observed SO42- concentrations.

The correlation plots in the East Central domain have squared correlation coefficients varying between 0.35 during quarter 1 to 0.61 during quarter 3 with r2 = 0.59 for the entire year. The simulation reproduced the average observed concentrations for quarters 1 and 4, but tended to underestimate them the rest of the year.

 Figure 9. The spatial domains used to average the simulated and observed SO42- concentrations and wet deposition rates. The symbol represents the IMPROVE and NESCAUM aerosol monitoring sites, the NADP wet deposition sites, and the center of receptor grid cells.

 Figure 10. Comparison of the daily averaged simulated SO42- concentrations, during 1992, to the NESCAUM/IMPROVE measurements over a New England and East Central domain, as defined in Figure 9.

Simulated Total SO42-Wet Deposition Rates. The average simulated and observed total SO42- wet deposition rates for quarters 1 and 3 are presented in Figure 11. The observed deposition rates during quarter 1 range between 0.6 and 1.8 g/m2/yr with the largest depositions occurring in a band from New York around the Carolinas and into the deep South. The simulation reproduces this pattern, but tends to underestimate the deposition rates along the coast and at the model domain boundaries. For example, the simulated deposition rates in Florida are about 0.3 g/m2/yr compared to 1 g/m2/yr for the observations.

 Figure 11. The simulated total SO42- wet deposition rates during quarters 1 and 3 1992, and the corresponding observed total SO42- wet deposition rates generated from the NADP observations.

During quarter 3, the model results properly simulate the high deposition rates south of Lake Erie and Ontario, where they exceed 4.2 g/m2/yr, and the decreasing deposition gradient extending from this region. However, like quarter 1, the simulation underestimates the deposition rates along the east coast, Florida, and the western part of the model domain between a factor of 1.5 and 2.

Comparisons of the simulated total SO42- wet deposition rates to observed rates in New England and East Central for weekly data are presented in Figure 12. This figure was created using the spatial domains defined in Figure 9. The seasonal squared correlation's for New England show high correlation ranging from 0.5 to 0.94 with an overall squared correlation of 0.79 for the year. The simulated deposition rates underestimated the observed by about 25% over the year, 1.71 g/m2/yr compared to 2.15 g/m2/yr for the observed rates. In the East Central domain, the seasonal squared correlation varied between 0.71 and 0.89 with a yearly correlation of 0.83. The yearly average deposition rates were 2 and 2.2 g/m2/yr for the simulated and observed data respectively.

  Figure 10. Comparison of the simulated weekly total SO42- wet deposition rates, during 1992, to NADP observations over a New England and East Central domain, as defined in Figure 9.

Back to Contents

Source Receptor Relationship

The source receptor relationship defines the attribution of a receptor's concentration to specific sources and/or source types. The Monte Carlo model was designed to be used as a diagnostic tool for the investigation and quantification of the source receptor relationship. This tool is based upon the fact, that in the Monte Carlo approach physical/chemical processes are simulated as discrete events. Thus, the influence of each process on the source receptor relationship can be analyzed separately.

It can be shown, that from the Lagrangian formulation of the dispersion equation (42, 43) the source receptor relationship can be defined as:


where i identifies the receptor volume
j identifies the source volume
k identifies the receptor time
l identifies the source release time (k - l= pollutant age)
the mass within the receptor volume i at time k due to emissions from source j released at time l
the mass released from source j at time l
the transfer matrix, the relative fraction of mass released from the source volume j at time l that impacts the receptor volume i at time k
Pt the transit probability
Pk the kinetic probability

This equation states that the pollutant mass at a receptor i and time k, due to emissions from source j released at time l, is proportional to the mass released at the source j at time l, where the proportionality constant is the element of the transfer matrix Tijkl. The physical interpretation of a transfer matrix element is that it represents the fraction of matter emitted from a source that reaches the receptor in the chemical form of interest. The total mass contribution to a receptor volume from a source can be found by simply summing over the source release time l.

Each element of the transfer matrix can be decomposed into a transit probability, Pt and a kinetic probability Pk. Pt defines the transport of the emissions, and can be thought of as the probability of a conservative species being transported from the source to the receptor. Pk defines the kinetic processes of the emitted pollutant mass. Pk for a primary species, such as SO2, is the probability that the pollutant mass will not be transformed or removed before being transported to the receptor, while for a secondary species, Pk is the probability that the pollutant mass will be transformed and not removed before being transported to the receptor. For example, Pk for SO42- is the probability that a quantum's initial SO2 mass is transformed to SO42- and the resulting SO42- mass is not removed before being transported to the receptor. For a conservative species, Pk is unity.

The source receptor relationship is presented graphically in Figure 13. The maps in this figure represent the three components of the source receptor relationship, mE, Pt, and Pk, and the resulting average source contributions to the SO2 and SO42- column concentrations at a receptor located in Massachusetts during quarter 3, 1992. The source grid is overlaid each map. The emission rate map shows that the highest values are over the Ohio River Valley where they can exceed 106 tons/yr. In the transit probability map, there is a decreasing probability of source impact with distance from the receptor. For example, sources in Michigan are about 7 times less likely to impact the Massachusetts receptor then the source located in Vermont. There is also a higher likelihood of emissions impacting the Massachusetts receptor from sources to the west and south than from the north. The SO2 kinetic probability map also displays an inverse relationship with distance. Sources in the vicinity of the receptor have Pk higher than 50% while in the Ohio River Valley the Pk is approximately 20%. The inverse relation with distance is expected since SO2 is a reactive primary species, and the residence time is roughly proportional to transport distance. The SO2 source attribution map shows that the concentration is dominated by contributions from nearby sources. The Massachusetts and Connecticut source grids supply over 35% of the SO2 mass. The large contributions from these sources are the result of their relatively high transit and kinetic probabilities. However, the high emissions in the Ohio River Valley lead to this area contributing almost 25% of the mass.

 Figure 13. The source attribution of the SO2 and SO42- column concentrations at a receptor located in Massachusetts. a) 1985 NAPAP Emission rates. b) Transit probabilities, Pt. c) SO2 kinetic probabilities, Pk. d) SO2 source attribution map. e) SO42- kinetic probabilities, Pk. f) SO42- source attribution map. The X on each map identifies the location of the Massachusetts receptor.

The kinetic probability map for SO42- presents a very different picture then seen for the SO2 . SO42- is a secondary species, and it takes time before a significant fraction of the SO2 in an airmass can be transformed to SO42-. This leads to kinetic probabilities that at first increase and then decrease with distance from the receptor. As shown, the Pk is about 2 to 5% for sources close to the receptor. This increases to almost 10% in western Pennsylvania and Northern Virginia, and again decreases to 2-5% in Illinois and Missouri. The SO42- source attribution map shows that the nearby sources contribute less than 10% of the SO42- mass to the receptor while the Ohio River Valley contributes nearly 30%.

Back to Contents

Summary and Discussion

A redesigned and re-calibrated Monte Carlo model for the simulation of regional scale transport, transformation, and removal of atmospheric pollutants was presented. The model was developed in a modular fashion, separating the calculation of the pollutant transport from the kinetics. The pre-computed transport data could then be used for the calculation of kinetic processes. This procedure was based upon the assumption that pollutant species do not influence meteorology to a significant degree. Transport calculations were conducted using 3-D gridded mean wind fields, and simulating horizontal and vertical diffusion. The kinetic processes were calculated through kinetic rate equations, where the rate coefficients were based upon space and time dependent environmental variables.

The model was used to simulate the SO2 - SO42- system over the Eastern US during 1992. Comparisons, of simulated seasonal SO2 concentrations to observations, showed that the model was able to reproduce the spatial patterns and magnitudes. Comparisons, performed with daily observations, showed that the model adequately simulated SO2 concentrations in New England with seasonal r2 = 0.56 - 0.83. In the East Central US, the seasonal r2 = 0.35 - 0.61. Comparisons of simulated total SO42- wet deposition rates to observed rates also showed strong agreement. Correlations of the weekly data in both New England and the East Central had seasonal r2 = 0.5 - 0.94.

A possible explanation, for the better correlations of the SO42- concentrations to observations in New England than the East Central domain, is that New England may be more dominated by regional scale phenomena than the East Central domain. There are few large sulfur emitting sources in the New England region. High and low concentrations of SO42- will tend to correspond to the transport of SO2 and SO42- from high and low emitting regions. The East Central domain, on the other hand, has large sources within and around the domain. Also, the Great Smoky Mountains are at the eastern edge of the domain. This is a complex situation, and the proper simulation of the SO42- concentrations is more likely to be dependent upon mesoscale transport and kinetics processes than in New England.

The Monte Carlo model was designed in such a manner as to allow for the investigation and quantification of the physical/chemical processes in the source receptor relationship. This was demonstrated by examining the role of source emissions, pollutant transport, and kinetics processes responsible for the source contributions of SO2 and SO42- to a Massachusetts receptor site during the quarter 3, 1992. It was found that the SO2 concentrations were highly influenced by nearby sources, such as those in Connecticut and Vermont. This was a result of the inverse relationship of the transit and kinetic probabilities with distance in the source receptor relationship. The high emissions in the Ohio River Valley also lead to significant contributions to the Massachusetts receptor concentration. The SO42- concentrations were dominated by emissions from more distant sources, particularly from the Ohio River Valley. This was due to the fact, that the concentrations of secondary species in a plume first increase and then decrease with time. Consequently, emissions from distant sources generally have less of a likelihood of impacting the receptor, but when they do, a higher fraction of the emissions is in the form of secondary pollutants than is the case with emissions from nearby sources.

Acknowledgment. Although the research described in this article has been funded in part by the U.S. Environmental Protection Agency under cooperative agreement CR 823756 and CR 818969, it has not been subjected to the Agency's peer and administrative review and therefore may not necessarily reflect the views of the Agency, and no official endorsement should be inferred. The information in this document has also been funded in part by Colorado State University Subcontract G-2887-4.

Back to Contents


1. Ottar, B. Ambio. 1976, 6, 203-206.

2. Hidy, G.M.; E.Y. Tong; P.K. Mueller. Design of the Sulfate Regional Experiment (SURE). Electric Power Research Institute, 1976; EC-125.

3. Pack, D.H.; G.J. Ferber; J.L. Heffter; K. Telegadas; J.K. Angell; W.H. Hoecker; L. Machta. Atmos. Environ. 1978, 12, 425-444.

4. Eliassen A. Atmos. Environ.. 1978, 12, 479-487.

5. Fisher, B.E.A. Atmos. Environ.. 1978, 12, 479-487.

6. Venkatram, A.; B. Ley; S. Wong. Atmos. Environ. 1982, 16, 249-257.

7. Draxler R.R.; Elliott W.P. Atmos. Environ. 1977, 11, 35-40.

8. Husar, R.B., D.E. Patterson, J.D. Husar, N.V. Gillani. Atmos Environ. 1978, 12, 549-568.

9. Husar, R.B.; D.E. Patterson; W.E. Wilson. A semi-empirical approach for selecting rate parameters for A Monte Carlo regional model. CAPITA Publication: Washington University St. Louis MO, 1986, 86r2.

10. Rolph, G. D.; Roland R Draxler; Rosa G. de Pena. Atmos. Environ. 1992, 26, 73-93.

11. Chang, J.S.; R. A. Brost; I.S.A. Isaksen; S. Madronich; P. Middleton; W. R. Stockwell; C.J. Walcek.. J. Geophys. Res. 1987, 92, 14681-14700.

12. Venkatram, A.; P.K. Karamchandani; P.K. Misra. Atmos. Environ. 1988, 22, 737-747.

13. Middleton P.; J. S. Chang; J.C. Del Conal; H. Geiss; J. M. Rosinkski. Atmos. Environ. 1988, 22, 1195-1208.

14. Patterson, D.E.; R.B. Husar; W.E. Wilson; L.F. Smith. J. Appl. Meteor. 1981, 20, 404-420.

15. Hall, C.D. Quart. J. Roy. Meteor. Soc. 1975, 101, 235-244.

16. Alsmiller, F.S.; Alsmiller R.G.; H.W. Bertini; C.L. Begovich. J. Appl. Meteor. 1979, 18, 17-26.

17. Boughton, B.A.; J.M. Delaurentis; W.E. Dunn. Boundary-Layer Meteorology. 1987, 40, 147-163.

18. Wiel, J.C. Atmos. Environ. 1994, 28, 3433-3448.

19. Luhar, A.K.; Rao, K.S. Atmos. Environ. 1994, 28, 3417-3431.

20. National Research Council. Science and Judgment in Risk Assessment. National Academy Press: Washington DC, 1994.

21. Hoke J.E.; N.A. Phillips; G.J. DiMego; J.J. Tuccillo; J.G. Sela. Wea. and Forecasting 1989, 4, 323-334.

22. Rolph G.D. National Climatic Data Center Report, 1992; TD-6140.

23. Pielke, R.A.; W.R. Cotton; R.L. Walko; C.J. Tremback; W.A. Lyons; L.D. Grasso; M.E. Nicholls; M.D. Moran; D.A. Ewlsey; T.J. Lee; J.H. Copeland. Meteorol. Atmos. Phys. 1992, 49, 69-91.

24. Lenschow, D.H.; Stephens, P.L. Boundary layer Meteor. 1980, 19, 509 - 532.

25. Willis, G.E.; Deardorff, J.W. J. Atmos. Sci. 1974, 31, 1297-1307.

26. Gillani, N.V. Atmos. Environ. 1978, 12, 569-588.

27. Gryning S.E., A.A.M. Holtslag, J.S. Irwin, and B. Sivertsen. Atmos. Environ. 1987, 21, 79-89.

28. Smith, F.B.; Hunt, R.D. Atmos. Environ. 1978, 12, 461-477.

29. Carras J.N.; Williams D.J. Atmos. Environ. 1980, 15, 2205-2217.

30. Munn, R.E.; Bolin B. Atmos. Environ. 1971, 5, 363-402.

31. Stocker, R.A.; M. Uliasz; R.A. Pielke. "Aerosol and Atmospheric Optics: Radiative Balance and Visual Air Quality, Volume A." In Proceedings of the International Specialty Conference, Snowbird, Utah, September 1994.

32. Schichtel, B.A. The retrieval of pollutant emission fields from ambient concentration precipitation chemistry data. Ph.D. diss., Washington University, 1996.

33. National Center for Atmospheric Research (NCAR), DS505.0 NCDC TD3240, US Controled Hourly precipitation, Daily 1948-Cont.

34. National Renewable Energy Laboratory (NREL). September 1992. User's Manual: National Solar Radiation Database (1961-1990), Volume 1.0. Distributed by National Climatic Data Center, Federal Building, Asheville, NC 28801.

35. Irving, Patricia M. Acidic Deposition: State of Science and Technology, Vol. I, Report 1, (Emissions). Government Printing Office: Washington, DC, 1991; 20402-9325.

36. Sisler J.F.; D. Huffman; D.A. Latimer; W.C. Malm; and M. Pitchford. Spatial and temporal patterns and the chemical composition of the haze in the Unites States: An analysis of data from the IMPROVE network, 1988-1991. CIRA, CSU, Fort Collins, CO, 1993; #ISSN No. 0737-5352-26.

37. Poirot, R. L.; P. J. Galvin; N. Gordon; S. Quan; A.V. Arsdale; R.G. Flocchini. "Annual and seasonal fine particle composition in the Northeast: second year results from the NESCAUM monitoring network." In the A&WMA conference, Vancouver, Canada, 1991; 91-49.1.

38. Husar R.B.; Elkins J.B.; Wilson W.E. "US visibility trends, 1960-1992, regional and national." In the 87th Annual Air & Waste Management Meeting, Cincinnati, OH, 1994.

39. US Environmental Protection Agency. AIRS User's Guide, Volume AQ1: AQS Data Dictionary, February, 1994, EPA-454/B-94-005.

40. National Atmospheric Deposition Program/National Trends Network. NADP/NTN Coordination Office, Natural Resource Ecology Laboratory, Colorado State University, Fort Collins, CO, 1993.

41. Husar, R.B.; Husar, J.D. Fine particle maps derived from regional PM2.5 and visibility data. Final Report to Radian Corporation, NC, 1995; EPA Contract Number 68-D2-0160, Work Assignment Number 2032.

42. Lamb, R.G.; Seinfeld, J.H. Environmental Science & Technology, 1973, 7, 253-261.

43. Husar, R.B. Possible remedies to "acid rain." CAPITA Publication: Washington University St. Louis MO, 1983; 83p4.

Submit your comments, feedback, questions, and ideas pertaining this page. Your input will be automatically added to the existing annotations. In order to add a new comment, you must be registered with the CAPITA People's Page.