A multimedia fate and chemical transport modeling system for pesticides: I. Model development and implementation

We have combined the US EPA MM5/MCIP/SMOKE/CMAQ modeling system with a dynamic soil model, the pesticide emission model (PEM), to create a multimedia chemical transport model capable of describing the important physical and chemical processes involving pesticides in the soil, in the atmosphere, and on the surface of vegetation. These processes include: agricultural practices (e.g. soil tilling and pesticide application mode); advection and diffusion of pesticides, moisture, and heat in the soil; partitioning of pesticides between soil organic carbon and interstitial water and air; emissions from the soil to the atmosphere; gas–particle partitioning and transport in the atmosphere; and atmospheric chemistry and dry and wet deposition of pesticides to terrestrial and water surfaces. The modeling system was tested by simulating toxaphene in a domain that covers most of North America for the period from 1 January 2000 to 31 December 2000. The results show obvious transport of the pesticide from the heavily contaminated soils in the southern United States and Mexico to water bodies including the Atlantic Ocean, the Gulf of Mexico and the Great Lakes, leading to significant dry and wet deposition into these ecosystems. The spatial distributions of dry and wet depositions differ because of their different physical mechanisms; the former follows the distribution of air concentrations whereas the latter is more biased to the North East due to the effect of precipitation.


Introduction
Pesticides are commonly used agrochemicals, but they have adverse effects on human health and the environment. Moreover, pesticides can be transported far away from their original application locations to contaminate important ecosystems such as the Great Lakes and the Arctic. Studies of the transport and fate of pesticides use a variety of approaches, but most current models have different limitations. For example, back trajectory models have no detailed soil simulations or atmospheric chemistry [1][2][3]; mass balance models [4] treat the heterogeneous atmosphere with uniform properties without considering realistic atmospheric conditions. A recent study coupled a three-layer fugacitybased soil-air exchange model with a simplified (12 layers from the surface to 7 km height) three-dimensional regionalscale dispersion model to address the long range atmospheric transport of pesticides [5,6] but did not consider physical and chemical interactions of pesticides with other atmospheric species in the atmosphere. This model underestimated annual dry toxaphene deposition to Lake Ontario [5] by a factor of more than 39 compared to observation-based data reported by Burniston et al [7]. Obviously, a more accurate modeling system is needed to simulate the fate and transport of pesticides and hence the deposition to sensitive ecosystems. Previous studies addressed environmental pesticides as an individual pollutant issue separately from other pollutant issues such as atmospheric aerosols, ozone and acid deposition; however, it is increasingly clear that atmospheric pollutants can interact with each other, and the interactions may affect the fate and transport of pollutants significantly. Since pesticides can exist in both gas and particle phases (i.e. on or in particulate matter) and can participate in atmospheric chemistry, their interactions with other natural and anthropogenic atmospheric species may affect their fate, budget, transport, chemistry, and deposition. Thus it is necessary to simulate pesticides along with other atmospheric species in the same modeling system so that the complex physical and chemical interactions among all the important species including pesticides can be accurately modeled. In the following, we will report the development and application of a single integrated model having capabilities to simulate all important processes in the soil, in the atmosphere, and on vegetation surface as well as agricultural practices such as toil tilling and pesticide application mode. To our knowledge, no other existing models can address all these processes. We based the development on toxaphene because of the information collected in the previous studies and also because toxaphene continues to be a widely concerned persistent, toxic, and bio-accumulative pesticide in many regions, and has been targeted for global elimination under the United Nations Environmental Program (UNEP) Stockholm Convention. The resulting model, however, can be used to predict the environmental behavior of a wide range of current use and banned pesticides (significant residues remaining in the soil due to past applications). The development of the modeling system is described in the present report, which is part 1 of a two-part series. Details of the model performance and a comparison with field measurements are reported in part 2 [8].

Modeling system description
To create the combined model, we coupled the US EPA's third generation SMOKE/CMAQ air quality modeling system [9,10] with the dynamic soil model: pesticide emission model (PEM) [11,12]. The combined modeling system can simulate all important processes involving pesticides, including the physical partitioning of the pesticide between the gas and aerosol phases and its chemical removal by reaction with OH. The OH concentration fields as well as those of the three aerosol modes are created by the model using criteria emissions inventories and thus have realistic spatial and temporal distributions, thereby increasing the accuracy of the pesticide transport and removal processes.
The PEM model includes additional processes such as soil tilling that brings pesticide residues from the deeper soil to the surface. This enhances the rate of volatilization to the atmosphere and provides a more realistic description of the pesticide loss from the soil. PEM is driven by the same meteorology that drives the atmospheric model, which allows it to reproduce the strong diurnal and seasonal emission variations observed for pesticides [11,12]. This is an improvement on fugacity-based models used by previous studies [5,6] that lack this agricultural and meteorological sensitivity. Figure 1 shows the overall flow diagram of the modeling system. The components are: MM5 (fifth generation mesoscale meteorology model [13,14], which uses the Pleim-Xiu land surface model [15][16][17]); MCIP (meteorology-chemistry interface processor) [18]; PEM [11,12]; SMOKE (sparse matrix operator Kernel emissions modeling system) [19]; CMAQ (community multiscale air quality modeling system [9,10]); MPIP (meteorology-PEM interface processor) and PSIP (PEM-SMOKE interface processor). MPIP and PSIP were developed for this project.

Model framework
MM5 produces meteorology for all models. MCIP processes the MM5 output to provide meteorological input for both SMOKE and CMAQ and computes dry deposition velocities for individual gas phase species. For this study, MCIP was modified to calculate dry deposition velocities for pesticides. MPIP prepares hourly meteorology from the MM5 output for use by the PEM model. PEM simulates the processes involving the pesticide behavior in the soil and on the surface of vegetation and calculates pesticide emissions. These emissions are then prepared by PSIP for input into the SMOKE emission processor, which simulates criteria pollutant emissions and combines them with pesticide emissions and prepares them for input into the modified CMAQ model that will be described in detail in the following sections.
In the toxaphene study, we ignore the effect of the pesticide air concentration (C air ) on the toxaphene flux from the soil, because toxaphene is very strongly bound to organic soil carbon [20] and the resistance to its volatilization to the atmosphere lies primarily in the soil. Thus the air concentrations do not affect the flux significantly except at very low soil concentrations such as those resulting from deposition and re-emission to/from non-agricultural soil. In the toxaphene study, therefore, the feedback from CMAQ to PEM involving C air was not implemented.

Soil process and pesticide emission modeling
PEM can estimate emissions both from current applications and from soil residues due to past applications. PEM can address three different modes of pesticide application: direct spray, in-soil application and seed treatment. The PEM model simulates such soil processes as vertical advection and diffusion of moisture, heat and pesticides in both bare and vegetated soils. Soil tilling is simulated by mixing the soil layers in the top 10 cm, which enhances the emissions into the atmosphere by bringing fresh pesticide residues to the soil surface. A canopy model is also incorporated in PEM, to account for such processes as spray interception by the vegetation canopy, with subsequent volatilization and wash off during precipitation events. For the toxaphene study, however, only soil residues are used in calculating emissions, because there was no usage in the domain during the study period.

Atmospheric chemistry and criteria emissions
Pesticides such as toxaphene can undergo four chemical processes in the atmosphere: reactions with O 3 , OH, NO 3 and photolysis. For most pesticides including toxaphene the reaction with OH dominates [21] and this was incorporated for this study into the CMAQ CB4 mechanism with linkages to both aqueous chemistry and aerosol processes. The CB4 mechanism is a lumped structure type that is the fourth in a series of carbon compound mechanism, and has detailed representation of organic compounds. The original chemical mechanism contains 36 species and 93 chemical reactions, of which 12 are photolytic. These lumped chemical reactions represent atmospheric interactions among compounds in the atmosphere, transforming reactants into intermediates and products, including OH radicals, which oxidize many pollutants such as pesticides in the atmosphere. The chemical reactions can also lead to formation of secondary aerosols, with which pesticides can be associated, from both inorganic and organic gaseous precursors. It is assumed that only gas phase toxaphene undergoes chemical reaction with hydroxyl radicals; particle phase toxaphene can become gas phase via partitioning mechanism (described in section 2.4). The rate constant used for the toxaphene reaction with OH is 2.65 × 10 −12 cm 3 molecule −1 s −1 , which is the average value of the rate constants for the reactions of the hexachloro and decachloro congeners of toxaphene [22].
In this study, 1999 US and 1995 Canadian criteria emission inventories were used to drive the CB4 chemical mechanism, evolving atmospheric chemistry that produces realistic oxidizing species and aerosols that can interact with pesticides. The US inventories (other than biogenic) were obtained from the US EPA [23]; US and Canadian land use data for calculating biogenic emissions as well as Canadian criteria inventories for other sources (i.e. area, mobile, and point) were obtained from the Ontario Ministry of the Environment.
These criteria emission inventories are processed by the SMOKE model to produce hourly, spatially gridded, and speciated criteria emissions for the specified domain. The SMOKE model is also used to merge the processed criteria emissions with pesticide emissions calculated by the PEM model to create input emission data for the modified CMAQ model.

Aerosols and gas-particle partitioning
In the present study, we have incorporated the K oa model [24][25][26] into the CMAQ system to simulate the partitioning of pollutants onto particulate matter (PM). The K oa model uses the octanol-air partition coefficient (K oa ) and fraction of organic matter in the aerosols ( f om ) to distribute semi-volatile material between gas and particle phases, using the customary parameterization: log K p = log K oa + log f om − 11.92. (1) Gas-particle partitioning is calculated according to the following equation: which expresses the fraction of pesticide on PM (φ) in terms of the mass density of total suspended particles (TSP, μg m −3 ) [24] and the gas/particle partitioning coefficient K p = C p /C g , (m 3 μg −1 ). C p and C g are, respectively, the pesticide concentrations in the particle phase (ng μg −1 ) and the gas phase (ng m −3 ). This model works with the aerosol module in CMAQ [27][28][29] to quantify the partitioning of pesticides to atmospheric particulate matter. The CMAQ aerosol module simulates PM in three log-normal modes: the Aitken ('I' mode) with diameters up to about 0.1 μm, the accumulation ('J' mode) with diameters between 0.1 and 2.5 μm and coarse particles having diameters between 2.5 and 10 μm. Recent field studies [30][31][32][33][34][35] suggest that coarse-mode aerosols can contain significant amounts of organic materials, so for the purpose of this study, we have assumed that the mass fraction of organic matter in the coarse particles is equivalent to that in fine particles (PM 2.5 ; particles smaller than 2.5 μm in diameter).

Dry deposition, wet deposition and transport
The deposition rate at which a pesticide is transferred to the surface of the earth is affected by many physical, chemical, and biogenic factors including properties of the pesticide, characteristics of the surface, and the strength of atmospheric turbulence. In CMAQ the dry toxaphene deposition rate is calculated by multiplying the dry deposition velocity (m s −1 ) with surface air concentrations (kg m −3 ). The basic equations to calculate dry gas deposition velocities can be found in Byun et al [18]. The dry deposition velocity in each grid cell is calculated using a series-resistance model in MCIP, which was modified for the purposes of this study by the inclusion of pesticides. Dry deposition of PM (and the associated pesticide) is also calculated using a deposition velocity approach. Size dependent deposition velocities for particles are determined by Brownian diffusion, gravitational settling, and turbulence, the details of which are in [27] and [29].
The methods applied in CMAQ for scavenging of pollutants by precipitation and wet deposition are described in [36] and [27]. The calculation of scavenging depends on whether the pollutant participates in aqueous chemistry. If the pollutant takes part in cloud water chemistry, the scavenging of this pollutant is dependent on Henry's law constant, dissociation constant, and cloud water pH; if the pollutant does not participate in cloud chemistry, the ending concentration and deposition amounts of the pollutant are estimated with the Henry law equilibrium equation. No aqueous phase cloud chemistry of pesticides is included in this study, but the uptake and concentration in cloud water are estimated using Henry's law. The temperature dependent Henry law constant for toxaphene used in this study is [37]: It is assumed that the accumulation-mode and coarsemode aerosols are completely absorbed by the cloud and rain water. Different from the other two modes, the Aitkenmode particles are considered as interstitial aerosols and are slowly absorbed into the cloud/rain water. Readers are referred to Binkowski [27] for detailed description of this process. Scavenging by snow and ice is not considered at this time.
The transport processes that primarily include vertical and horizontal advection, vertical and horizontal diffusion, and the mixing of pollutants by the parameterized subgrid-scale clouds [9] are all included in the simulation. The transport processes were formulated in conservation (i.e., flux) forms for the generalized coordinate system.

Implementation for the testing of the modeling system
The modeling system is tested by carrying out a study that investigates environmental behavior of toxaphene in North America. A large amount of toxaphene was used from the 1940s to the 1990s in North America, primarily in cottongrowing regions in the southern United States. Toxaphene was banned in the United States in 1982 and all use there ceased by 1986. It was not licensed for use in Canada and was deregistered in Mexico in 1993 [38,39]. The historical applications, however, left significant soil residues in the regions of heaviest use, and these residues continue to release toxaphene. Concern about the effect of this on distant watersheds, especially the Great Lakes, has led to several studies [3,40,41] that identified long range transport of toxaphene as a potential problem for several ecosystems.
The relatively localized application of toxaphene provides a geographically well-defined source in which the quantities used and the times of application were studied. This information along with available field measurements of toxaphene in North America provides an ideal means to test and evaluate the ability of the multimedia model we have developed to predict concentrations and depositions.

Model domain and study period
To address the issue of long range transport, the model domain (see figure 2)

Boundary and initial conditions
Boundary conditions for toxaphene were determined as follows: MacLeod et al [4] show that the minimum values of measured median air concentrations of toxaphene in the Arctic Archipelago and Yukon River-Aleutian regions are about 5 pg m −3 ; this value is considered to be the background air concentration and is used as the boundary condition for the northern, eastern, and western boundaries of the gridded domain.
The southern boundary condition was based on the field measurements reported in Hoh and Hites [41] and the study of Scholtz and Bidleman [20]. The Hoh and Hites measurements used for this purpose were conducted at a location in Louisiana, which is only 17 km from the Gulf of Mexico in which toxaphene had never been used. Samples of gas phase toxaphene collected at this site every 12 days from February 2002 to September 2003 yielded an arithmetic average toxaphene concentration of 61 ± 7 pg m −3 . Based on these data and other measurements in the southern United States, Scholtz and Bidleman [20] estimated an average background toxaphene concentration of 67 pg m −3 in this region and this concentration was used as the southern boundary condition for our study.
The initial and boundary conditions for criteria pollutant species were based on various measurements and were set as default values in the original CMAQ, and the initial conditions for toxaphene throughout the domain were set to zero.

Toxaphene residues
There is no current toxaphene usage in the study domain, but there are significant residues in the soil in both the United States and Mexico from historical applications, and emissions from these soil residues are the major source of atmospheric toxaphene in North America. Toxaphene residue data for the United States were obtained from Environment Canada [39,42]. Most toxaphene in the United States was limited to cotton-growing states in the south east, where toxaphene soil residue concentrations can be as high as 270 kg km −2 (about 350 × 10 3 kg/grid square). Gridded toxaphene residues for Mexico were calculated from the Mexican toxaphene usage estimated by MacLeod et al [4] and an assumed soil half-life of four years [3]. The residues in Mexico were allocated to the model domain using the gridded percentage of cropland [42] as a surrogate for spatial distribution. Figure 2 presents the merged US and Mexican residues in the study domain, showing that the highest toxaphene soil residue concentrations in some locations in Mexico are of similar magnitude to the concentrations in the southern United States. The soil residues in Canada are negligible since toxaphene was not licensed for use in this nation.

Physical and chemical properties
Toxaphene is a complex mixture of chlorinated bornanes having more than 32 000 congeners of which more than 600 have been identified by chromatography [43]. However, it is not practical to simulate so many individual congeners. Therefore, this study simulates the pesticide as a single species using the physical and chemical properties for technical toxaphene, which are listed in table 1.

Domain-wide results and discussion
We used the model to assess the continental-scale transport of toxaphene for a one-year period. In figure 3, we show the results for the modeled annual mean gas phase (figure 3(a)) and particle phase ( figure 3(b)) toxaphene concentrations for the year 2000. The similarity between the distributions of air concentration (figure 3) and residues (figure 2) confirms that North American soil residues are the major source for the atmospheric toxaphene concentrations in the region. The concentrations in different parts of the domain vary over about four orders of magnitude, reflecting the large differences in emissions from the soil, reaching 1.0 × 10 −7 ppmV in Alabama-the region that has the highest soil residues-and decreasing with distance from the most heavily contaminated soils.
While the concentrations are highest in the regions of heaviest application, there are nevertheless significant concentrations of toxaphene-especially in the gas phaseover the Gulf of Mexico and the Atlantic ocean. These are obviously due to atmospheric transport during the oneyear simulation. Such extensive transport of the semi-volatile material in the gas phase is another noteworthy result of this study. Measurable-and in most cases significant-gas phase concentrations are transported throughout the eastern half of the continent and for hundreds of kilometers offshore from the source. The agreement between model predictions and landbased field measurements that was described in the companion paper [8] suggests that this prediction of the concentrations over the coastal waters is reasonably accurate and is thus only about one to two orders of magnitude less than the concentration at the point of heaviest application. It further reveals that in a period of one year, toxaphene emissions can affect the entire Gulf of Mexico, parts of the Caribbean and the Continental Shelf up to about 40 • N.
There are some differences in the spatial distribution between gas (figure 3(a)) and particle ( figure 3(b)) phases; the toxaphene concentrations in particles seem to be more spread compared to those in the gas phase. This is because particle phase toxaphene is assumed not to undergo chemical reactions, but gas phase toxaphene takes part in atmospheric chemistry that is an important removal process. The differences in the spatial distribution confirm that the atmospheric transport and removal of the pesticide are influenced by gas-particle partitioning.
The dispersion of the pesticide through the atmosphere of course leads to its deposition in those locations to which it is transported. The model predictions of this deposition for the year 2000 are shown in figure 4 for the dry (a) and wet (b) deposition mechanisms. The spatial distribution of dry deposition follows that of the air concentration, so it also affects the same parts of the Gulf of Mexico and the Atlantic Ocean mentioned above. Figure 4(a) also shows there is significant dry deposition into the Great Lakes area. The model predicts that the dry deposition into the lower Great Lakesespecially lake Erie and parts of lake Ontario-can be of the order of 50-100 mg ha −1 yr −1 , which is only about one or two orders of magnitude less than that at the locations of application.
The spatial distribution of wet deposition (figure 4(b)) differs somewhat from that for dry deposition because of the effect of precipitation. It is more biased to the North East and thus affects the Great Lakes and the offshore Atlantic more strongly than the Gulf of Mexico and the Caribbean. The absolute amount deposited to the Great Lakes by wet deposition, however, is similar to that from dry deposition, which can be on the order of a few tens of milligram/hectare for the year 2000.
The present study considered only emissions resulting from the soil residues remaining from the original application of the pesticide, but the predicted deposition rates are significant, suggesting that future studies may also need to consider re-emissions if a more complete picture of the pesticide distribution is to be obtained. This study shows that water bodies such as the Great Lakes, the coastal areas of the Atlantic Ocean and the Gulf of Mexico are being contaminated by deposition of toxaphene. Total deposition to the lower Great Lakes was occurring at the rate of 50-100 mg ha −1 during the 2000 calendar year. The further spread to high latitudes is also implied by these results.