Monte carlo simulation of regional aerosol: transport and kinetics

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.

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.

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

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:

(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 (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 . 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² m²/s in stable conditions and 10⁶ m²/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. = 0 m²/s. The right hand figure contains the trajectories and height history when a diffusion coefficient of = 10⁴ m²/s 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 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.