During the 1970’s, air quality measurements had shown 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 design 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 (Ottar, 1976), in the US the Sulfate Regional Experiment (SURE) was conducted (Hidy, Tong, and Mueller, 1976), along with numerous plume studies.

Associated with these studies were the development of models for the simulation of regional scale dispersion. The initial models tended to be simple, using aggregated and parameterized meteorology and rate coefficients to simulate the physical/chemical processes. These included trajectory models (Pack, 1978; Eliassen, 1978), statistical models (Fisher, 1978; Venkatram et al., 1982), as well as simple diffusion and box models (Draxler, 1977; Husar et al., 1978; Patterson et al., 1981). 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 (Husar et al., 1986; Rolph and Draxler, 1990).

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) (Chang et al., 1987) and the Acid Deposition and Oxidant Model (ADOM) (Venkatram et al., 1988). These models attempt to simulate atmospheric processes, like that 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 (Middleton et al., 1988). This confines these models to be operated in a "black box" approach.

In this chapter, 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 (Patterson et al., 1981; Husar et al., 1986). The current 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.

In the Monte Carlo approach to air pollution simulation, individual quanta 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 of the simulation. The quanta are often referred to as particles; in this paper these two terns will be used interchangeably.

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 for 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 3-1). An airmass is depicted as 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 transition probabilities. 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.

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 upon receptor concentrations.

Monte Carlo simulations of atmospheric processes were developed in the 1970’s (Hall, 1975; Alsmiller et al., 1979). However, this technique lost favor in the 1980’s as other approaches were pursued. Due to the multiple benefits of this modeling approach, it has again increased in popularity during the 1990’s for mescoscale simulations (Boughton et al, 1987; Luhar and Rao, 1994; Wiel, 1994), and is believed to offer considerable improvements over conventional Gaussian-plume models (National Research Council, 1994). 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. (1981) and Husar et al. (1986).

The independence of transport from chemical processes allowed for the separate development of a meteorological module, for the calculation of pollutant dispersion, and a kinetic module for the calculation of pollutant concentrations and deposition rates. The flow of raw meteorological, emission, chemical, etc. data through these two modules is presented in Figure 3-2. The meteorological module processes the input data creating a Lagrangian airmass history database. This is accomplished by computing the transport of conservative species, simulated by quanta, 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, (Hoke et al., 1989; Rolph, 1992) and Colorado State University’s Regional Atmospheric Modeling System, RAMS (Pielke et al, 1992). The meteorological database is then converted into Lagrangian airmass histories by the Eulerian/Lagrangian transformer. Along specified time increments of the quantum’s trajectory, its 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.

The kinetic module combines a source oriented Lagrangian airmass history database with emissions rates, surface meteorological data, etc. 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 calculated 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 saved out. 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.

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 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 (AL). The ABL experiences diurnal (Figure 3-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 that 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 that develop above a surface that is considerably warmer than the overlying air. Over the upper part of the boundary layer, thermals may lose their buoyancy but retain sufficient momentum to distort the capping inversion, entrain overlying cooler air, and then fall back into the mixed layer (Lenschow, 1980).

Figure 3-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, 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 (Willis and Deardorff, 1974). 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 (Gillani, 1978; Gryning et al., 1987).

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 (Smith and Hunt, 1978; Carras and Williams, 1980). It is this repeated process of vertical mixing followed by vertical wind veer that is essentially responsible for regional dispersion.

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:

(3-1)

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 (3-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 (Hoke et al., 1989; Rolph, 1992) and Colorado State University’s Regional Atmospheric Modeling System (RAMS) (Pielke et al., 1992). The random component of the velocity field occurs at a resolution smaller than the meteorological model grids used to generate wind fields. These fluctuating components are sub grid phenomena and calculated from the available meteorological variables.

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 . 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 ranging between 10² m2/s in stable conditions and 10⁶ m2/s during intense convective activity have been used in the past (Patterson et al., 1981; Munn and Bolin, 1971). 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.

The dispersion of a puff of pollutants emitted from a source in western Tennessee is illustrated in Figure 3-4. This figure was created by simulating the release of pollutants from the source by 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. = 0 m²/sec. The right hand figure contains the trajectories and height history when a diffusion coefficient of = 10⁴ m²/sec was employed for both day and night.

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 3-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 of 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 3-5. This figure displays the position of particles on 01/16/1992 noon Greenwich mean time (GMT) that 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 (Stocker et al., 1994).

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

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 semi-empirical, in that it uses deterministic rate equations for the calculation of kinetic processes applied to each quantum. The simulation of the kinetic processes is empirical, in the sense that the rate coefficient are dependent upon the meteorological and chemical environment of the quanta, and the rate coefficient equations are determined through a fitting process matching simulated to observed concentration and deposition fields.

The semi-empirical approach for the generation of rate equations is not original. It was the underlying feature of the OECD modeling effort (Eliassen, 1978; Eliassen and Saltbones, 1983), and was also used in the original Monte Carlo Model (Husar et al., 1986). The methodology of the semi-empirical approach applied in this study is outlined below. This methodology was first put forth by Husar et al., (1986) for the tuning of the original Monte Carlo model.

The semi-empirical approach involves a loop consisting of four procedures (Figure 3-6). The first procedure is the collection of input data sets. These data sets provide the necessary meteorological, emission, and measured concentration data for the generation and verification of simulated pollutant species. The next step is to formulate the model, creating the rate equations dependent upon the input data that can, in principle, reproduce the observed field and lab observations. The third step is to determine any free constants in the rate equations. These constants are selected such that the model results best reproduce the observed data. The last step involves the analysis of the results, and their comparison to the observed data. In this process, patterns of deviations are sought that can help explain the residuals. This in turn, provides insights into revising the formulation of the rate equations, and potential problems or deficiencies with any of the first three steps. With the new insights, any of the previous three procedures can be returned to and revised. The process is complete when no new information can be discovered from the residuals.

In the execution of the approach, the following guidelines are followed: 1) The process is always begun with real world data that are used as a guide for all other steps. 2) Keep the model as "simple" as possible, such that mathematical complexity does not obscure the physical relationships, or interfere with the ability to perform what if analysis. 3) Discard any added formulations of the rate coefficient equations that do not improve the model fit to data or the understanding of the physical and chemical processes governing the kinetics.

This approach for the generation of rate coefficient equations has several advantages. First, at any stage of the development there is a functional tool that can provide best guess estimates. Deviations between the current simulation and the measured data provide clues and possible new insights as to how the physico-chemical frame could be improved. Thereby, potentially generating new knowledge. The model is forced to fit a set of observed data, therefore it provides sophisticated interpolation and extrapolation of the tuning concentration data in space and time. Also, given adequate quantity and quality of data the rate coefficient equations may evolve from primarily empirical to a mainly physico-chemical capable of predications.

The disadvantages of this approach are that the resulting rate coefficient equations may not generate unique rate constants or even be correct. Also, the rate coefficient equations are most suitable to the conditions for which they were developed. For example, if the chemical mix of the atmosphere changes markedly, then the kinetic rate constants derived previously may be inapplicable.

The approach above was applied to the formulation of kinetic rate equations for the simulation of SO₂ and SO₄^2- concentration and deposition fields over the Eastern US during 1992. Below is a description of the four procedures involved in the tuning process outlined above, and the looping required to arrive at a set of suitable rate equations.

The first step in the semi-empirical approach for the generation of the rate equations is the collection of relevant data sets. The data requirements can be classified into three categories, meteorological, emissions, and ambient monitoring data.

NGM Meteorological Data. The primary meteorological data set used to drive the Monte Carlo model were obtained from the National Meteorological Center's Nested Grid Model (NGM). The NGM model is an operational forecasting model that provides a variety of objective weather analyses and forecast products on a number of different grid scales. It was developed at the National Meteorological Center (NMC) (Hoke et al., 1985; Hoke et al., 1989, Kanamitsu, 1989, Junker et al., 1989). The model is routinely operated and is used to support weather forecasting over the US. and the globe.

Precipitation. The proper simulation of the wet deposition of pollutants is highly dependent upon the precipitation field. One of the parameters of the NGM meteorological data is the precipitation rate. This precipitation field was used in the initial SO₂-SO₄^2- simulations where the simulated wet deposition rates were found to be very poor. This prompted a verification study of the NGM precipitation field. The study was accomplished by comparing the NGM precipitation fields to measured hourly data from the US Control COOP Hourly Precipitation database (NCAR, 1995). The COOP database consists of approximately 2900 stations located throughout all 50 states in the US.

Figures 3-8 and 3-9 present the seasonal precipitation fields from the NGM and measured data respectively. Average precipitation rates for a New England, Industrial Midwest, and Southeast region are listed in Table 3-2. During the winter seasons, Q1 and Q4, the NGM precipitation follows the spatial patterns of the measured data over the Eastern US. However, the Q1 precipitation rates are underestimated throughout the East up to a factor of 2. During Q2, the NGM data underestimate the precipitation rates in New England, ~0.6 m/yr compared to 1 m/yr, but in the south, the measured precipitation rates were slightly overestimated. The precipitation rates compare the least favorably during Q3, where the NGM data exhibits a significant decreasing gradient from the south (>3 m/yr) to the North (~ 1 m/yr). The measured data rates are rather uniform throughout the East at about 1.5 m/yr.

The varied spatial and temporal deviations between the NGM and measured precipitation rates precluded the use of the NGM data as the sole source of precipitation data for the derivation of rate coefficients and pollutant simulations. A new precipitation field, suitable for model input, was created by summing the measured precipitation rates over a two hour period, then spatially interpolating the data to the same grid used by the NGM meteorological data. The measured data covered only the continental US. The NGM precipitation data was used to fill in Canada, Mexico, and off the coast of the US.

Total Cloud Cover. The total cloud cover was used in the formulation of the SO₂ oxidation and dry deposition rates. The total cloud cover data were derived from the National Solar Radiation Database (1992). This database contains the synoptic surface reports from approximately 290 stations reporting synoptic weather over the United States every hour. Figure 3-10 presents the noon cloud cover on May 8, where it is seen that in the Midwest the sky was clear but it was cloudy over the East. The total cloud cover data was transformed to the same stereographic grid used by the NGM meteorological data. This was accomplished using a technique similar to that used for the precipitation data. The grid points outside of the US were set to 60%, the average cloud cover over the US during 1992.

Lagrangian Airmass History Database. The Lagrangian airmass history database was created from the NGM meteorological database. The modeled sources were placed in the center of each grid cell over the US, southern Canada, and Northern Mexico, essentially creating a uniform source field (Figure 3-11). Each modeled source then accounts for all emissions from point and area sources within the grid cell. The database contained the three dimensional location of each particle, along with the specific humidity and temperature, at two hours time increments. Each particle was tracked for seven days, or until it was transported past the edge of the meteorological grid.

The 1985 NAPAP Annual emission inventory (Irving, 1991) for North American was used for the model tuning. The original data were provided as point and area emissions on a grid with resolution of 1 degree Longitude and 2/3 degree Latitude grid. The data were transformed to the meteorological grid by summing all point and area emissions together for all emission grid cells that fell into the same meteorological grid cell. No attempt had been made to account for variations in 1992 emission from those in 1985. However, it is estimated that the overall US SO₂ emissions decreased about 5% since 1985 (US Environmental Protection Agency, 1994).

The ambient concentration data constitute the most important data sets for the tuning of the model. The observed data utilized for the tuning process consisted of ambient fine sulfur aerosol concentrations from the IMPROVE (Sisler et al., 1993) and the NESCAUM (Poirot et al., 1991) monitoring networks, airport visibility data (Husar et al., 1994), and the SO₄^2- wet deposition data from the National Atmospheric Deposition Program/National Trends Network (NADP/NTN, 1993).

Sulfate Concentrations. The sulfate concentration data were obtained from the IMPROVE (Interagency Monitoring of Protected Visual Environments) and NESCAUM (Northeast States for Coordinated Air Use Management). The IMPROVE network has been operational since March 1988, as a cooperative effort between the US EPA, federal land management agencies, and state air agencies. The network has 55 sites located in National Parks, Wilderness areas, and rural area meeting the IMPROVE site protocol, 13 of which are east of the Mississippi. The network collects aerosol samples over a twenty four hour period twice weekly. The NESCAUM network consisted of 10 stations located in the northeastern US. These sites were located in rural areas found to be consistent with the IMPROVE site protocol (Poirot et al., 1991, Schichtel and Husar 1992). The network used the same sampling method as IMPROVE, but collected samples three times a week.

The 1992 seasonal fine sulfur aerosol from the IMPROVE and NESCAUM networks is presented in Figure 3-13. With the exception of the locations in Northern Minnesota, the fine sulfur concentrations are larger in the summer (Q2 and Q3) than the winter (Q1 and Q4). All seasons display the largest concentrations for the locations in the Central East, Tennessee through Virginia, and decrease in all directions from that area. The region of highest concentrations coincides with the highest SO₂ emission rates (Figure 3-12).

In the East, there are large areas that do not contain IMPROVE or NESCAUM monitoring sites, particularly in the Great Lake region and the South. To overcome the lack of stations, airport visibility data were used as a surrogate for the spatial interpolation of the fine sulfur aerosol between the monitoring sites, while maintaining the actual concentrations at the monitoring sites. This technique was first employed by Husar and Husar (1995) for the derivation of seasonal fine particle maps over the continental US. The airport visibility data came from the same surface observation network as the total cloud cover data (National Solar Radiation Database, 1992). Figure 3-14 presents the final seasonal isopleth maps of sulfate generated from the fusion of the aerosol and visibility data.

Wet Deposition Data. The SO₄^2- wet deposition data were obtained from the National Atmospheric Deposition Program/National Trends Network (NADP/NTN, 1993). This data contains the contributions from both SO₂ and SO₄^2- washout of the atmosphere. This network has 257 monitoring sites located throughout the US with every state having at least two monitoring sites. The network collects precipitation samples over a weeks time continuously, and has been operating since 1978.

Seasonal SO₄^2- wet deposition fields are presented in Figure 3-15. As shown, the smallest deposition rates occur during the winter seasons Q1 and Q4 with deposition rates ranging between 1.25 and 2 g/m²/yr in the East. The largest deposition rates occur during Q3 where they are greater than 3.5 g/m²/yr in Pennsylvania and New York. This high deposition region is slightly to the north east of the largest emission rates (Figure 3-12). The high deposition region coincides with high precipitation rates (Figure 3-9) where they are greater than 1.5 m/yr, while the precipitation rates in the high emission region are between 0.5 and 1.5 m/yr.

The kinetic processes governing the atmospheric concentrations of SO₂ and SO₄^2- are the transformation of SO₂ to SO₄^2- and their removal by dry and wet deposition. The rate of change of the SO₂ mass is:

(3-2)

where , , and are the oxidation, dry deposition, and wet deposition rate coefficients of SO₂ respectively. Similarly, the rate of change of the SO₄^2- mass is:

(3-3)

where and are the dry and wet deposition rate coefficients of SO₄^2- respectively.

In the SO₂ to SO₄^2- rate equations, there are five rate coefficients that must be defined via the iterative tuning process. The development of the final set of rate coefficient equations followed a long arduous path involving many iterations. In each iteration, the simulation was compared to the observed data via comparisons of daily, weekly, and seasonal SO₄^2- concentration and deposition fields, site by site correlation plots, and regional averages and correlation plots. These iterations resulted in the gathering of new data, reformations of the rate equations, and refitting of free coefficients. The following will not recount all of the steps and reasons leading to the final equations. However, some of the intermediate steps resulting in the more important insights will be shown.

The initial rate coefficients used were those developed by Husar et al., (1986) where they tuned the CAPITA Monte Carlo model using data for the Eastern US during 1977 - 1978. These rate coefficients with units of 1/hr are:

k_t= 0.002 + 0.036 SR (3-4)

SR = (1-.7TS)GSR (3-5)

k_d2 = V_d / H (3-6)

V_d = 0.002 + 0.024 * SR (3-7)

k_d4 = k_d2/10 (3-8)

k_w2 = W_s2/H_S * P = 75 * P (3-9)

k_w4 = W_s4/H_S * P = 200 * P (3-10)

where SR ground level solar radiation [kW/m²]
TS total cloud cover
GSR ground level solar radiation in the absence of clouds [kW/m²]
V_d deposition velocity [cm/s]
H = mixing height [m]
W_s2, W_s4 washout ratios for SO₂ and SO₄^2- respectively
P precipitation rate [m/hr]
H_S Scale Height

It was assumed that all SO₂ and SO₄^2- mass arose from the simulated sources only, thus, background SO₂ and SO₄^2- were neglected. The surface concentrations were calculated for receptor volumes defined by the meteorological grid cell and the mixing height.

The rational for using the solar radiation to parameterize the oxidation of SO₂ to SO₄^2- was based upon the work of Husar et al., (1978) and Wilson (1981). Examining oxidation rate coefficients estimated from plume studies, they found that much of the variation between the estimated rate coefficients could be reduced by normalizing the data by the ground level solar radiation. In the initial formulation of the transformation rate coefficient, k_t, it was modulated by the solar radiation between a "background" oxidation rate of k_tb = 0.2 %/hr, and a maximum oxidation rate of k_tmax = 3.6 %/hr was used. The transformation rate coefficient equaled the background rate during the night and winter, while it equaled the maximum rate during clear mid afternoon summer days.

Although solar radiation is an important part of the oxidation of SO₂, in these equations it is merely a scaling parameter that provides diurnal, seasonal, and latitudinal variation of k_t. The ground solar radiation is calculated from the extraterrestrial solar radiation at the top of the atmosphere using standard procedures (Kondratyev, 1969).

The dry deposition of SO₂ is determined by atmospheric mixing, which delivers material to the surface, and by the absorptive properties of the surface. A convenient formulation for the dry deposition constant is k_d = V_d/H. H is the depth of the mixing layer, and V_d is the deposition velocity. The SO₂ dry deposition velocity varies spatially, diurnally and seasonally. Results from a number of experimental studies compiled in Garland (1978) and in a review of acid deposition by the National Research Council (1983) lists values of the SO₂ V_d between 0.1 - 2.2 cm/s depending on season, time of day, and surface conditions. The maximum V_d occurs midday during the summer, and the minimum values occur during the night and winter. To account for these variations, Husar et al., (1986) parameterized V_d using the solar radiation. They used a background V_d of 0.2 cm/s, and the maximum V_d of 2.4 cm/s.

The dry deposition velocity for SO₄^2- is limited by the aerosol fluid dynamical properties, and the SO₄^2- V_d is usually smaller than the SO₂'s (National Research Council, 1983; Seinfeld, 1986). Field experiments of the V_d of submicrometer diameter particles compiled in National Research Council (1983) has the SO₄^2- V_d ranging between 0.04 - 1 cm/s. In Husar et al., (1986), as well as in this study, the SO₄^2- V_d was set to one tenth of SO₂'s. In the Monte Carlo simulations, dry deposition was only applied to quanta below the mixing layer. Except for slight modification of the maximum deposition velocity constant, this equation was not changed in the tuning process.

The wet deposition rates developed in Husar et al., (1986) were based upon a probability of precipitation parameter. In this study, the washout ratio was employed. The washout ratio is a common parameterization of wet removal rates in dispersion models (Barrie, 1978; Barrie, 1981; Eliassen and Saltbones, 1983; Seinfeld, 1986), and is defined as:

W = = (3-11)

The wet deposition rate coefficient is then defined as:

k_w = (3-12)

where P is the precipitation rate and H_S is the scale height.

The theoretical washout ratio for SO₂, W_S2, ranges between 2300 and 2 * 10⁶ at a temperature of 0 ^oc and a pH range of 3 to 6 (Barrie, 1981). The washout ratio of sulfate aerosol is about an order of magnitude higher than SO₂ on average (Garland, 1978), with values between 3*10⁵ and 8*10⁵ in eastern Canada (Barrie, 1981), and average values of 10⁶ in Europe (Whelpdale, 1981).

In the initial simulations, the washout ratios were constant with W_S2 = 75*10⁵ and W_S4= 200*10⁵. Also, the scale height was set to a constant value of H_S = 1 km. A scale of 1 km was found to be the typical value estimated from simulated scale heights over the Eastern US using the Monte Carlo model. The wet deposition rates were applied to all quanta in a precipitating grid cell.

Examination of the results from the simulation using the above rate coefficients, revealed a systematic spatial and temporal variation between the simulated and observed SO₄^2- precipitation weighted concentrations. This is presented in Table 3-3, which contains the average seasonal SO₄^2- concentrations for a New England, Industrial Midwest, and Southeast region. During every quarter, the ratio of the simulated to NADP concentrations in New England are 10 to 60% lower then in the Industrial Midwest and Southeast, and the ratios in the Industrial Midwest are 2 - 30% larger than in the Southeast. Also, ratios during the winter, Q1 and Q4 are about 30% larger than during the summer, Q2 and Q3. These results indicate that the washout ratio should be larger during the summer than the winter, and larger in New England than the Industrial Midwest and Southeast.

A similar spatial/ temporal pattern to the wet deposition data was also seen in the SO₂ surface and column concentrations (Table 3-4). The winter season has larger SO₂ concentrations than the summer, and New England has the smallest concentrations. Consequently, the SO₂ concentrations could be used as a scaling parameter of the washout ratios to remove the deviations in the wet deposition results.

A means of using the SO₂ concentration as a scaling parameter for the SO₂ washout ratio was developed from approximations of the SO₂ solubility equations by Barrie (1981). The equilibrium solution chemistry of SO₂ is described by the following equations:

SO2(g) + H2O(l) « SO2^.H20(aq) (3-13)

K_H = [SO2^.H20] / [SO₂]_g

SO2^.H20(aq) « H⁺(aq) + HSO3- (aq) (3-14)

K₁ = [H⁺] [HSO₃^-](aq) / [SO2^.H20]

HSO3- (aq) « H⁺(aq) + SO32- (aq) (1-3) (3-15)

K₂ = [H⁺] [HSO₃^-](aq) / [HSO₃^-]

K_H is the Henry's law coefficient and K₁ and K₂ are the first and second dissociation coefficients. The square brackets represent molar concentrations. The SO₂ in solution is distributed among the three species SO₂^.H₂0_, HSO₃^- and SO₃^2- in a manner that is strongly dependent on the pH and weakly dependent on temperature. Between typical atmospheric pH’s of 3 to 6, over 90% of dissolved SO₂ exists as bisulfite, HSO₃^-. Therefore, Equation 3-15 can be neglected, and the washout ratio can be expressed as (Barrie, 1978; Barrie, 1981):

W_S2 =[HSO₃^-]_aq/[SO₂]_g = k₁k_H/[H⁺]_aq (3-16)

This equation is dependent upon the assumption that the SO₂ in the air is in equilibrium with the SO₂ in the rain. Hales (1978) has shown that this is a valid assumption for long range transport.

If it is assumed that SO₂ is the predominant acidic substance, then Equations 3-13 and 3-14 yields (Barrie, 1981):

[H⁺]_aq = [HSO₃^-]_aq = (k₁k_H [SO₂]_g)^1/2 (3-17)

and the washout ratio becomes:

W_S2 = (k₁k_H/[SO₂]_g)^1/2 (3-18)

This equation provides the sought after scaling of the washout ratio by SO₂ concentrations. It was implemented into the Monte Carlo model using the column concentration of SO₂ via the equation:

W_S2 = (A/([SO₂CC]_g/H_S + B))^1/2 (3-19)

where A is the optimum constant value of k₁k_H, and B is a constant representing background acidity not accounted for by the SO₂ concentration. The constants A and B were determined by fitting the simulated to observed data, assuming a constant scale height of 1 km. The optimum values were found to be A = 0.63 (mol/l_aq)²/(mol/l_g) and B = 2 mg/m². The rate coefficient for the SO₂ wet removal then became:

= W_s2/H_S * P = (40000 / ([SO₂CC]_g + 2))^1/2 * P (3-20)

where [SO₂CC]_g has units of mg. The SO₄^2- wet removal rate coefficient equation was not changed. However, W_S4/H_S was increased from 200 to 250.

Table 3-5 presents the seasonal average wet deposition rates for the same sub-regions of the Eastern US as in Table 3-3, using the new SO₂ and SO₄^2- wet removal rate coefficient. As can be seen, the seasonal variation in the ratio of the simulated and NADP concentrations in table 3-3 has been removed. Also, some of the spatial variation has been removed. The Industrial Midwest and Southeast both have ratios of ~1, but the simulation in New England still underestimates the NADP concentrations by about 30%.

The oxidation of SO₂ to SO₄^2- occurs in both gas and aqueous phase chemical reactions. Aqueous phase oxidation of SO₂ takes place in hydrometers, such as clouds and fog. If the hydrometer is a precipitating cloud, then the products of the reactions will be removed from the atmosphere. If the cloud or fog evaporates, then the dissolved sulfur dioxide along with the resulting oxidation products will be released to the atmosphere. It is estimated that for every precipitating cloud, six have formed and evaporated (Hobbs 1994). Aqueous phase oxidation mechanisms generally have high oxidation rates, often greater than 10 %/hr and can be as high as 100% an hour. Consequently, the cycling of cloud formation, SO₂ oxidation, and cloud evaporation can be a significant oxidation mechanism in the atmosphere (Hegg and Hobbs, 1981; Calvert et al., 1983; Calvert et al., 1985; Hegg and Hobbs, 1989; Hidy, 1994).

The transformation rate coefficient defined by Equation 3-4, does not separate the gas and aqueous phase oxidation mechanisms, but uses SR to estimate an overall oxidation rate. However, SR does not provide any indication of SO₂ interaction with hydrometers. In fact, in equation 3-4, the transformation rate decreases with increasing cloud cover, not accounting for the enhanced oxidation rates due to higher likely hood of cloud interaction. To remedy this, the transformation rate coefficient was reformulated to include the parameter specific humidity. Specific humidity is the mass of water per unit mass of moist air, m_w/m_a. It is a conservative property so long as neither evaporation nor condensation occurs in the air, and is independent of temperature, pressure and volume. It is believed that the specific humidity is a valid parameter for the indication of aqueous phase transformation, because both the dew point temperature and the cloudbase decrease with increasing specific humidity (Humphreys, 1964). This results in a higher likelihood of cloud interaction with a polluted airmass.

Using specific humidity, the transformation rate coefficient became:

= k_tb + SR (k_{tSR max} + k_{tSH max} SH/SHmax ) (3-21)

where k_tb is the back ground transformation rate
k_{tSR max} the maximum transformation rate due to SR alone
k_{tSH max} the maximum transformation rate due to the specific humidity.
SH the specific humidity
SHmax the maximum specific humidity in the Eastern US, SHmax = 0.023

The specific humidity term in the equation is multiplied by the ground level solar radiation, because the aqueous phase SO₂ oxidation is dependent upon atmospheric concentrations of H₂O₂. Atmospheric concentrations of H₂O₂ are primarily the result of gas-phase reactions driven by solar radiation, and aqueous phase SO₂ oxidation generally decreases at night and in the winter relative to summer days (Calvert et al., 1985).

The specific humidity does not differ substantially for an airmass within or outside of a cloud. Therefore, the specific humidity may be an indicator of aqueous phase oxidation, but as it is used in Equation 3-21 it is a bulk parameter accounting for both gas and aqueous phase chemistry.

The addition of the specific humidity to the transformation rate coefficients improved the simulation of sulfate in all of the measured parameters. In Figure 3-16 and 3-17, the daily averaged simulated SO₄^2- are seasonally compared to observations in correlation plots. These plots were created by averaging the model results and the NESCAUM/IMPROVE observations over the New England and Central East spatial domains as defined in Figure 3-18. The transformation rates used in the simulated data in Figure 3-16 was = 0.001 + 0.036 SR, while = 0.001 + SR (0.03 + 0.03 SH/0.023) in Figure 3-17. All other rate coefficients were identical. As shown, the correlation coefficients between simulated and observed data increased every quarter with it increasing from r² = 0.31 to r² = 0.39 in the year correlation plot of Central East.

Further examination of the SO₄^2- concentration residuals revealed an underestimation of the sulfate concentrations during and after many precipitation events. It was speculated that this was a result of underestimating the aqueous phase oxidation that must have occurred due to cloud interaction. To account for this, the transformation coefficient was further modified by adding the term k_{tprec max} *SR*PF, where k_{tprec max} is the maximum transformation rate during precipitation events, and PF is a precipitation flag. k_{tprec max} was set to 10%/hr

Figure 3-18. The spatial domains used to average the simulated and observed SO₄^2- concentrations and wet deposition rates. The symbol n represents the IMPROVE and NESCAUM aerosol monitoring sites, the NADP wet deposition sites, and ¡ the center of simulated receptor grid cells.

The final relations derived for the rate coefficients with units 1/hr are:

= 0.001 + SR (0.02 + 0.05 SH/SHmax + 0.1*PF) (3-22)

SR = (1-.65TS)GSR (3-23)

= V_d / H (3-24)

V_d = .22 + 2 * SR (3-25)

= /10 (3-26)

= W_s2/ H_S * P = Sqr(40000 / (SO₂CC + 2)) * P (3-27)

= W_s4/H_S * P = 250 * P (3-28)

where SR = ground level solar radiation [kW/m2]
TS = total cloud cover
GSR = ground level solar radiation in the absence of clouds [kW/m²]
SH = specific humidity [kg Water/kg Air]
SHmax = maximum specific humidity [kg Water/kg Air]: SHmax = 0.023
PF = precipitation flag, i.e. PF = 1 if it rains PF = 0 otherwise
V_d = deposition velocity [cm/s]
H = mixing height [m]
W_s2, W_s4 = washout ratios for SO₂ and SO₄^2- respectively
SO₂CC = SO₂ column concentration [mg/m²]
P = precipitation rate [m/hr]
H_S = Scale Height [m]: H_S = 1000 m

Implicit in these rate equations are nonlinear kinetic rate processes, such as the transformation of SO₂ to SO₄^2-. This transformation is dependent upon the ambient concentrations of SO₂ as well as OH and H₂O₂ and other species. If these rate coefficient are applied to regions and times, other than the Eastern US during 1992, then the concentrations of the species that the rate coefficients are dependent upon will most likely change. If the changes are significant, then it may be necessary to retune the constants in the equations along with possible changes in the formulations of the rate coefficient.

The rate coefficients averaged over each month in Vermont, Ohio, and Alabama and the input meteorological parameters are presented in Figure 3-19. The averaging was performed by weighting each rate coefficient by the SO₂ or SO₄^2- mass of the quanta that it was applied to. The transformation rate coefficient displays a strong seasonal pattern increasing 4.5 times from an average rate of 0.2 - 0.3%/hr in January to 0.8 - 1.4%/hr in July (Figure 3-19A). The seasonality is driven by the solar radiation which peaks in June and the specific humidity which peaks later in the summer (Figure 3-19C). During the summer months, the precipitation and specific humidity components comprise more than 50% of the transformation coefficient . However, during the winter months the solar radiation and background components dominate the transformation rate, particularly in Vermont. The transformation rate also is spatially dependent, increasing from North to

South (Figure 3-19A). The contributions from the specific humidity and precipitation displays more of the spatial variation than from the solar radiation component. The summer transformation rate due to the specific humidity and precipitation terms was ~0.4%/hr accounting for 50% of the transformation rate in Vermont. This increased to ~1 %/hr accounting for 70% of the transformation rate in Alabama.

The removal rate coefficients show little seasonal variation at the three stations (Figure 3-19B). However, the wet removal rates have high month to month variation due to variable precipitation (Figure 3-19D). The lack of seasonality in the dry deposition rates is due to the offsetting seasonality of the dry deposition velocities and mixing heights (Figure 3-19E). This is because, the dry deposition velocities were parameterized by the solar radiation, and the mixing height is driven by the solar radiation (see section 3.2.1). The SO₂ dry deposition rates are higher then the wet deposition rates, 1.7 - 3.2 %/hr compared to 0.2 - 1.3 %/hr, with the highest rates occurring in Alabama. The primary removal of SO₄^2- is via wet deposition, where the rates vary from about 1 to 5%/hr.

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. Below, the simulated seasonal SO₄^2- concentration and wet deposition fields are presented and compared to the observed data. The wet deposition data contains contributions from SO₂ and SO₄^2-. Also, correlation plots of the simulated and observed concentration and wet deposition rates are presented for the New England and Central East regions defined in Figure 3-18.

Simulated SO₄^2- concentrations. Figure 3-20, presents the seasonal simulated SO₄^2- concentration fields with Q1 corresponding to January - March. During the winter, Q1 and Q4, the simulated SO₄^2- is fairly uniform over the East with typical concentrations between 2.5 and 4.5 mg/m³. The highest concentrations are found in a belt from eastern Ohio to the Chesapeake Bay, where they can be greater than 5 mg/m³. The winter simulation reproduces the general spatial patterns of the observed data in Figure 14. However, the highest observed concentrations are around Alabama and Mississippi (>4 mg/m³), as opposed to Virginia for the simulation.

The simulations during the summer, Q2 and Q3, have more variability than the winter seasons, with typical concentrations between 4 and 7 mg/m³. Concentrations above 5.5 mg/m³ are found in a belt from the Chesapeake bay to Tennessee, with the highest values in the Ohio River Valley, ~7.5 mg/m³. The summer simulation reproduces the general spatial pattern and magnitudes of the observations in Figure 3-14. However, the highest observed concentrations are to the south and east of the simulation, ranging from southern Pennsylvania and northern Virginia to western Kentucky and Tennessee. The concentrations for both the observed and simulated SO₄^2- progressively decrease in all directions from the high concentration band producing similar concentrations over New England (~ 3 mg/m³), the South (~5 mg/m³), and into the Midwest (2 - 5 mg/m³).

The daily averaged simulated SO₄^2- are seasonally compared to observations in correlation plots for the New England and Central East regions (Figure 3-21). The simulated SO₄^2- 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 SO₄^2- concentrations. The correlation plots in the Central East domain have r² between 0.35 during Q1 to 0.61 during Q3 with r² = 0.59 for the entire year. The simulation reproduced the average observed concentrations for Q1 and Q4, but tended to underestimate them the rest of the year.

Simulated SO₄^2- Wet Deposition Rates. The seasonal simulated SO₄^2- wet deposition rates are presented in Figure 3-22. These deposition rates contain contributions from both SO₂ and SO₄^2- wet deposition. The Q1 simulation has deposition rates ranging between 0.5 and 2 g/m²/yr with the largest depositions occurring in a band from New York around the Carolinas and into the deep South. The simulation matches the observed data well (Figure 3-15), but tends to underestimate the deposition rates along the coast. The Q2 simulated deposition field is very similar to Q1, but the deposition rates are greater, ranging between 1.2 and 2.8 g/m²/yr. The Q2 simulation reproduces the bulk of the observed spatial pattern. However, the deposition rates are underestimated around the Great Lakes. The largest disparity occurs in New Jersey and New York where the observed deposition rates are about 3.5 g/m²/yr compared to less than 2 g/m²/yr for the simulation.

During quarter 3, the model results properly simulate the high deposition rates south of Lake Erie and Ontario, where they exceed 4.2 g/m²/yr, and the decreasing deposition gradient extending from this region. However, like the other quarters, the simulation underestimates the observed deposition rates along the east coast, Florida, and the western part of the model domain between a factor of 1.5 and 2. The Q4 simulation and observed deposition rates are similar to Q3, but they are smaller, with the largest values less than 2.5 g/m²/yr in the high deposition region south of Lake Erie and Ontario.

Comparisons of the simulated total SO₄^2- wet deposition rates to observed rates in New England and Central East for weekly data are presented in Figure 3-23. 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/m²/yr for the observed rates. In the Central East 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/m²/yr for the simulated and observed data respectively.

Using the airmass history database, and the SO₂ - SO₄^2- kinetics derived in the previous section, the two components of the transfer matrix, transit probability, P_t and kinetic probability, P_k were calculated for SO₄^2- surface concentrations and the sum of the SO₂ and SO₄^2- wet deposition rates. This was accomplished using the procedure outlined in Section 2-5. The probabilities were calculated for each source grid cell modeled in the airmass history database (Figure 3-11). The receptors occupied the same grid cells used as the sources, and their volumes were defined by the grid cell area and mixing height. The resulting transfer matrix database consisted of 504 sources and receptors. The receptor time and age each had two hour increments, and the maximum age was 7 days. This database was then aggregated as needed following the methods discussed in section 2-3.

In Figures 3-24 - 3-31, Q2 receptor oriented transfer matrices are presented for SO₄^2- surface concentrations and the sum of the SO₂ and SO₄^2- wet deposition rates for seven receptor sites located throughout the contiguous US. The sum of the SO₂ and SO₄^2- wet deposition rates will be referred to as the SO₄^2- wet deposition rate. These receptor oriented transfer matrices display the average relative source contributions to the SO₄^2- concentrations and wet deposition rates at each receptor. Included in each plot is the total relative source contribution calculated by integrating the transfer matrix over the spatial domain. Assuming that all sources have the same emission rates, then a receptor with a total relative source contribution that is two times larger than another receptor's would also have twice as much mass at the receptor.

The transfer matrix can be used to identify the pathway of the airmasses impacting the receptor (Ashbaugh et al., 1985; Poirot and Wishinski, 1986; Vasconcelos, 1995). If the airflow to a receptor predominately comes from one direction, then the sources along this pathway will have a larger relative impact than other sources. Inherent to all transfer matrices is the fact that nearby sources have a larger relative impact then distance sources. This creates sharply decreasing gradients from the receptor site, as seen in the left hand plots of Figure 3-24A, and washes out the contribution of distant sources reducing the ability to identify airmass pathways. Poirot and Wishinski (1986) called this bias central tendency. Based upon geometric reasons, they determined that weighing each transfer matrix element by the distance from the source to the receptor removes the central tendency, as seen by the right hand plots of Figure 3-24A. With the central tendency removed it is possible to identify sources that have large contributions to the receptor for other reasons then their proximity to the receptor.

The transfer matrices show that there is a predominate airmass pathway to each receptor that allows a different set of sources to contribute to the receptor concentrations and deposition rates. The airmass pathway varies considerable between the surface concentrations and wet deposition rates. At the Eastern receptor sites, sources impacting the SO₄^2- concentrations are over a broad area and tend to lie to the northwest and southeast of the receptor site. The sources to the northwest have higher relative contribution then other sources. At the two receptors in the West, Idaho and New Mexico (Figures 3-26B and 3-27), sources strictly from northwest contribute to the receptor concentrations. The sources contributing to the receptor’s wet deposited mass, Figures 3-28 - 3-31, are generally to the south and southwest of the receptors except in Idaho. This implies that the precipitating airmasses almost strictly come from the south and southwest while dry airmasses arrive from more expansive regions. This pattern is typical for all quarters particularly for receptor located east of the Rocky mountains.

The role of the kinetic rates in the transfer matrices is also evident in these figures. Examination of the transit probability, P_t, component for the Minnesota and Illinois receptors reveals that an airmass is more likely to come from the south than indicated in Figure 3-24. But as shown by the wet deposition transfer matrices, Figure 3-28, airmass from the south encounter more precipitation events than airmass from the north. Precipitation efficiently removes the SO₂ and SO₄^2- mass causing the southern sources to have smaller contributions to the receptor’s SO₄^2- concentrations then the northern sources.

The role of the kinetic rates is also evident in the total relative source contributions. The total relative source contribution varies little for the receptors east of the Rocky Mountains at approximately 0.0032 (Figures 3-24 - 3-26), but it is smaller for the New Mexico and Idaho receptors at 0.0016 (Figures 3-26B- 3-27). The rather constant Eastern total relative source contribution is due to the transformation and removal kinetics not vary substantially with space in the East (Figures 3-19). In the West, the transformation rates are lower than the East. This results in smaller relative source contributions compared to those in the East. However, the two western sites are close to the edge of the modeling domain, reducing the number of sources that can impact the receptors. This artificially reduces their total relative source contributions.

The total relative source contributions vary substantially for the wet deposition rates. In Oklahoma it is 0.009 (Figures 3-30A) compared to 0.0027 in Illinois (Figures 3-28B), with the transfer matrix value greater than 0.000016 1/km² for most of the source impact domain to Oklahoma. The different relative source contributions is primarily due to the varying precipitation rates and frequencies at the receptors. At the Oklahoma receptor the precipitation rate is 1.9 m/yr compared to 0.4 m/yr in Illinois (Figure 3-9), and it rained 26% of the time in Oklahoma compared to 15% in Illinois (Figure 32).

1

This chapter presented the CAPITA Monte Carlo model. A Lagrangian model for the simulation of the transport, transformation, and removal of atmospheric pollutants. The model was developed in a modular fashion, separating the calculation of the pollutant transport from the kinetics. Transport calculations were conducted using 3-D gridded mean wind fields, and simulating horizontal and vertical diffusion through transition probabilities. The kinetic processes were calculated through kinetic rate equations, where the rate coefficients were based upon space and time dependent environmental variables and determined through a tuning process.

The model was tuned for the development of SO₂ - SO₄^2- kinetic rates in the Eastern US. The best rate coefficients had transformation rates that exhibit substantial seasonal variation and a mild north-south spatial gradient. Monthly averages of k_t varied between 0.2 - 0.3 %/hr in January and ~0.8-1.4 %/hr in July, with the larger transformation rates occurring in the Southeast. It was found that the best SO₂ washout ratio was inversely dependent upon the square root of the SO₂ concentration. This dependence is based on the inverse solubility of SO₂ with solution acidity. This non linearity resulted in the washout ratio in Vermont being 35% smaller than in Ohio where the simulated SO₂ concentrations were among the highest. The final rate coefficients were able to simulated the spatial patterns and magnitudes of the SO₄^2- concentrations and wet deposition rates. Correlations comparing simulated to observed SO₄^2- concentrations in a New England region for each season were r² = 0.56 - 0.83, and r² = 0.35 - 0.61 in a Central East region. Seasonal correlations between simulated and observed SO₄^2- wet deposition rates in both New England and Central East were between r² = 0.5 - 0.94.

Transfer matrices were generated for SO₄^2- surface concentrations and wet deposition rates. Examination of Q2 transfer matrices over the US revealed distinct variations between the concentration and wet deposition transfer matrices. The transfer matrices for each receptor exhibited a preferential flow with precipitating airmasses coming almost exclusively from south and southwest of the receptor, while the dry airmasses came from a broad region biased to the northwest of the receptors. Total relative source contributions were relatively constant for SO₄^2- concentrations, but varied substantially for the wet deposition transfer matrices due to varied spatial precipitation rates.