- Email: [email protected]
Development of a model for reactive emissions from industrial stacks
Environmental Modelling & Software 20 (2005) 1239–1250 www.elsevier.com/locate/envsoft
Development of a model for reactive emissions from industrial stacks Luis E. Olcese, Beatriz M. Toselli* Departamento de Fı´sico Quı´mica/INFIQC, Facultad de Ciencias Quı´micas, Universidad Nacional de Co´rdoba, 5000 Co´rdoba, Argentina Received 1 May 2003; received in revised form 17 June 2004; accepted 18 August 2004
Abstract We have developed a model, CAPAS, capable of estimating short-term concentrations of primary and secondary pollutants resulting from point source emissions. The model is designed to simulate the complex interaction of plume dispersion and non-linear chemistry. The main features of the model are: (1) the horizontal and vertical resolution within the plume, which offers a realistic treatment of the entrainment process, (2) its flexibility with regard to choices of chemical kinetic mechanism and the inclusion of a powerful numerical solver to handle stiff systems, and (3) its simplicity of operation and adaptability to solve new problems. The realism of the plume model simulations was tested by comparing model calculations of plume concentrations with data of SF6 tracer concentrations and ozone concentrations. The results of the model compare well with both sets of field measurements. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Dispersion models; Plume reactivity; Air pollutants; Non-linear chemistry; Stiff solvers
Software availability Program title: CAPAS Developer: L.E. Olcese Contact address: Departamento de Fı´ sico Quı´ mica/ INFIQC, Facultad de Ciencias Quı´ micas, Universidad Nacional de Co´rdoba, 5000 Co´rdoba, Argentina. E-mail: lolcese@ fcq.unc.edu.ar. Tel.: C54-351-433-4169; fax: C54-351-433-4188. First year available: 2000 Hardware requirements: PC Software requirements: FORTRAN 90 compiler * Corresponding author. Fax: C54 351 433 4188. E-mail address: [email protected] (B.M. Toselli). 1364-8152/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2004.08.008
Program language: Availability:
FORTRAN free of charge, upon request to the authors
Abbreviations and symbols massi,j,k:
mass stacki: cell length: sy,j: s starty,j: s endy,j: sz,j:
Mass of the specie i in the sub cell j in the cell k (If k is not present, it means that it is applied to any cell) Mass of the specie i emitted by the stack Length of every cell Pasquill sigma for horizontal dispersion at sub cell j Pasquill sigma for horizontal dispersion at the beginning of the sub cell j Pasquill sigma for horizontal dispersion at the end of the sub cell j Pasquill sigma for vertical dispersion at sub cell j
1240
s startz,j: s endz,j: absi: volj,k:
vol under surfacej: exchangej:
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
Pasquill sigma for vertical dispersion at the beginning of the sub cell j Pasquill sigma for vertical dispersion at the end of the sub cell j Absorption coefficient in the surface for specie i Volume of the sub cell j in the cell k (If k is not present, it means that it is applied to any cell) Volume of the sub cell j that is located under the surface Volume exchanged between the sub cells j and jC1
1. Introduction Air quality models can be classified into two broad groups: one that includes an extensive treatment of the chemical reactions, mainly oriented to predict formation of secondary pollutants; and another group oriented to simulate dispersion of primary pollutants in different situations (street canyons, industrial plumes, etc.). Usually, chemical industries, power plants, refineries, etc., use plume models like ISCST3 (US EPA, 1995a,b) to estimate the air pollution impact of their production processes. ISCST3 is perhaps the most widely used model for estimating near-field concentrations of non-reactive pollutants. It is a steady-state Gaussian plume model, i.e. dispersion parameters such as meteorological conditions and emission rate are assumed constant while the plume is in route from the source to the receptor. Some recent applications based on Gaussian models can be found in the works of Goyal and Sidharta (2004) and Canepa and Ratto (2003). The main drawback of these models is the inability to handle formation of secondary pollutants like ozone, or to handle the rate of disappearance of the emitted specie when the rate is different from a first order decay. This assumption is valid for CO or SO2, but is erroneous for species like NOX or VOC in urban atmospheres, or in polluted rural areas. In these cases, the formation of secondary pollutants becomes important, and the use of another type of models is necessary. Reactive plume models are often used to estimate the local or short- to medium-range (i.e. up to a few hundred km) impacts of power plants or smelters on air quality. Issues of interest typically include ozone and particulate matter concentrations above the ambient air quality standards of a region, visibility degradation and acid deposition. These models use detailed gas-phase mechanisms to describe the chemical transformations. The particular case of simulating the dispersion of reactive species emitted by an industrial stack has been treated only in few models, mainly the Reactive Plume Model (RPM) (Stewart and Liu, 1981). The RPM
version IV (Morris et al., 1992) is a standard regulatory model used for estimating pollutant concentrations in the atmosphere, resulting from the emissions from point sources such as industrial stacks. It uses either point source emission estimates or initial plume concentrations as inputs, and calculates downwind concentrations, as the plume expands. RPM-IV makes use of a diffusion/reaction equation similar in form to the atmospheric diffusion equation, but with dispersion parameters that evolve with time. This model considers a control volume that moves downwind along the plume trajectory, and solves the diffusion/reaction equation by considering the nonlinear photochemistry and diffusion of chemical species from one cell into another cell. Other existing models that are useful to simulate dispersion of reactive plumes are ROME (Gabruk et al., 1999) and ADMS3 (Carruthers et al., 1994), but this one is not freely available. In this work we present a new model (CAPAS), capable of simulating reactive plumes, easy to use by an operator with little training and easily adaptable to solve new problems or to use new chemical mechanisms. The model code is written in standard FORTRAN 90 and it has a modular structure, thus the replacement of a particular module (i.e., the numerical integrator of the ordinary differential equations) by another one is straightforward. CAPAS model shares some characteristics with the RPM model, but includes several new ones, as it will be shown below.
2. Model formulation This model is based on the Lagrangian approach to the turbulent diffusion (Seinfeld, 1986). We have chosen this approach, and not the Eulerian one, due to the simplicity of the resulting approximations. Briefly, the model consists of an array of cells, each one similar to those used in the box models; the cells are arranged in the direction of the x axis in a continuous way and the shape of each cell is an elliptic truncated cone. The size of a cell is given by four times the Pasquill sigmas in the y and z axes (both different in the beginning and in the end), and user definable in the x axis, typically 5–30 m depending on the wind speed. Each cell is later divided in 10–30 annular and equi-spaced sub cells, with the concentration of the species assumed to be homogeneous in each sub cell. A scheme of the cells and sub-cells can be viewed in Fig. 1. The boundary conditions are defined by an external sub cell with the concentrations of the species taken from the atmospheric concentrations. The steps to set up the model can be summarized as follow: 1. All the gases are injected in the atmosphere without any initial interaction. The gases exiting from the
1241
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
mass of each species is recalculated (Section 2.8). The fraction of the mass of the gases in the region of a sub cell in the previous cell that is occupied by a different sub cell in the new cell is added to the new sub cell. 9. Concentrations of the species are updated depending on the chemical mechanism (Section 2.9). 10. Steps 6–9 are repeated until the final distance is reached.
Several cells
A common problem of many air pollution models is that due to instabilities of the numerical solvers and to the rapid change in the species concentration because of variations in the volume of the cells, some intermediate transient concentrations are negative. We have found that a way to avoid such spurious results is to set up the model such that works internally with the mass of each species, instead of concentration.
Sub cells in a cell
2.1. Initial distribution
Lateral view
Front view
Fig. 1. Schematic of the physical representation of the plume used in the CAPAS model to simulate the chemistry and dispersion of the pollutants.
stack are assumed to be well represented by a Gaussian distribution and located in the first cell; the initial distribution of pollutant mass emitted by the industrial stack is given by the error function in two dimensions (Section 2.1). 2. The gases temperature is the same as the exit temperature of the stack. 3. The final height of the plume is calculated (Section 2.2). 4. The meteorological conditions are assumed constant with time, over the time period of transport from the source to the receptors. 5. The ordinary set of stiff differential equations representing the chemical mechanism is solved to obtain the new values of the species concentration in order to advance the solution. In the case that the chemical mechanism includes photolysis reactions the zenith angle dependence of these reactions is incorporated (Section 2.3). 6. After this first cell, new atmospheric gases are added to the outermost sub cell, and the new temperature is calculated (Sections 2.4 and 2.5). 7. The cell grows in size by the increase in the sy and sz values, and the sub cells grow proportionally (Section 2.6). 8. Depending on the volume of the sub cell (Section 2.7) and the absorption coefficient on the surface, the
A bi-dimensional Gaussian distribution of the species is assumed between the boundaries of the initial cell. The mass of each species in each sub cell is given by the equation: massi;j;1 Z mass stacki ! cell length sy;j
!
pffiffiffi sy;outermost j ! 2 !! sy;j pffiffiffi erf sy;j1 ! 2
! erf
ð1Þ
The outermost ring concentrations are the same as the atmospheric concentrations. The lack of interaction between the plume and the atmosphere in the first cell is justified considering the turbulent nature of the gases coming from the stack, and the short time elapsed between the exit of the gases from the stack and the formation of the plume. 2.2. Final height of the plume The final height of the plume is calculated according to the classical equations proposed by Briggs (1969, 1971, 1974), and revised by Turner (1997). The complete set of equations used is shown in Appendix A. 2.3. Photolysis rate constants When the chemical mechanism includes photolysis reactions, the photolysis rate constants are previously calculated using the TUV 4.1 program (Madronich, 1993), at ten values of zenith angle (0 –10 –20 –30 –40 –50 –60 –70 –78 –86 ). The calculated constant
1242
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
for each photochemical reaction at each angle are entered as input data, and depending on the zenith angle at each cell, the corresponding photolysis rate constant is calculated using a cubic spline interpolation. The computational time needed to calculate the photolysis rate constants is very small, making unnecessary to consider uniform photolysis rate constants along the entire plume. Photolysis rate constants are calculated assuming clear sky conditions, with no clouds and aerosols. However, different cloud cover conditions and aerosol loading can be easily incorporated into the TUV model and therefore included in the CAPAS model.
m1 Z
number of species X
massi;j;k massi;j;k
iZ1
exchangej1 volj
ð4Þ
m2 Z
number of species X
massi;j;kC1 massi;j;k
iZ1
exchangej voljC1
ð5Þ
tempj Z
m1 tempj Cm2 tempjC1 m1Cm2
ð6Þ
2.4. Mass exchange 2.6. Sigma y and sigma z Due to the increase in volume of the plume, when the gases pass from the cell k to the cell kC1, the mass present in each sub cell must be adjusted. The model calculates the volume exchanged between sub cells j and jC1 by calculating the volumes of the sub cells from the innermost until the jth one, minus the same volume but in the previous cell. This can be expressed as: exchangej Z
numberX of sub cells
numberX of sub cells
jZ1
jZ1
voljC1
volj
ð2Þ
Consequently, the mass of a given species in the sub cell j is: massi;j;k Z massi;j;k1 Cmassi;jC1;k1 massi;j;k1
exchangej voljC1;k1
exchangej1 volj;k1
ð3Þ
In the first cells, due to the rapid change in the sub cells volume, the volume exchanged can be greater than the volume of the cell, leading to negative volumes. To overcome this problem, the model selects a new time step and calculates the new volume. If the new calculated volume is greater than zero, is accepted and the calculation continues. If not, the model iteratively decreases the step size until the calculated volumes are positive. 2.5. Temperature variation The temperature of each sub cell is adjusted in every step of the calculation, averaging the temperature of the fraction of the total mass of the sub cell that is not exchanged and the temperature of the fraction that is injected to the sub cell, weighted by the mass. In each cell, the temperature of the outer sub-cell remains constant. For a given sub cell j of a cell k, its temperature is calculated according to:
Both sigmas are calculated according to the procedure described by Pasquill (1961) depending on the atmospheric stability. Urban or rural conditions are used depending on the land use. If different conditions are found along the plume trajectory, the predominant environment is used. 2.7. Sub cell volume Depending on the position of the sub cell relative to the ground, there are three different possibilities to calculate the volume of each sub cell: (a) The sub cell does not reach the surface in any point. (b) The sub cell does reach the surface since its initial x. (c) The sub cell does reach the surface in a point between its initial and final x. The mathematical expression of each volume is obtained by analytical integration. In case (a), the volume is calculated as an annular elliptical truncated cone. In case (b), the volume of the sub cell under zZ0 is subtracted from the calculated in case (a). Case (c) is treated as a combination of the two previous cases, dividing the sub cell in two sections; the first one between the initial distance x and the point where the sub cell does reach the ground, and the second one from there to the final distance x. The volumes of each section are calculated in the same way as cases (a) and (b). The mathematical expressions for each case are included in Appendix B. 2.8. Absorption coefficient The interaction of each species with the surface (cement, grass, sand, etc.) the products formed and the rate constants of the involved reactions (dry and wet deposition, heterogeneous reactions, etc.) are very complex phenomena, and depend on many variables
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
(soil type, temperature, humidity, etc.). As a simplification, three cases are possible for each species: total absorption, partial absorption, and total reflection. To simulate these three situations in this model, we used a coefficient abs, specific for each species. This factor can take values ranging from 1 (the species is completely absorbed in the surface) to 0 (completely reflected). In each sub cell, massi;j Z massi;j massi;j ! absi
vol under surfacej volj
ð7Þ
The experimental difficulty of determining the value of the absorption factor, abs is one of the major handicaps of this model. Nevertheless, the value for the coefficient abs can be fitted to the observed atmospheric concentrations of selected species using iterative methods.
1243
2.10. Concentration at a given receptor The final goal of this air pollution model is to estimate the concentration of selected pollutants at some defined places, called receptors. The concentration of a specie at a given receptor, defined by x, y, z coordinates, is obtained looking at the concentration in the sub cell that included the receptor. Depending on the downwind distance to the stack, the wind speed, and the length of the cell the model calculates inside which cell a given receptor is located. To determine in which sub cell is placed the receptor, the angle formed between the receptor and the center of the plume is calculated together with the distance of the border of each sub cell. If the receptor is located outside the limits of the plume, atmospheric concentrations are used.
2.9. Chemical mechanism 3. Validation One of the most widely used numerical procedures to solve stiff systems of ordinary differential equations (ODEs) is the Gear’s method (Gear, 1971) with automatic time step selection and error control. The method is however; time consuming because it requires inversion of large matrices. This is the method used by RPM-IV model. To improve the performance of the model, the numerical integration subroutine used by CAPAS is based on the Rosenbrock method with Kaps–Rentrop coefficients (Bader and Deufhlhard, 1983), modified to take into account the sparsity of the Jacobian matrix (Duff and Reid, 1984). This (semi-)implicit method for treating complex non-linear kinetic mechanisms has been selected after a detailed study of the different available solvers, and has been chosen due to its stability and velocity (Olcese and Toselli, 1998). We have tested the numerical integrator with different chemical mechanisms proposed by Hov et al. (1989) and Seinfeld (1986), and the CBM-IV mechanism (Gery et al., 1989) obtaining good stability in all of the mechanisms under different initial conditions. For a complete description of the numerical solver and the different tests performed the interested reader should consult the work by Olcese and Toselli (1998). The selection of a given mechanism and the numerical method to solve the stiff ODEs is crucial to the performance of the model. Roughly, 90% of the computational time is spent in the numerical integration of the chemical mechanism. To test the CAPAS model we have chosen for the simulations the three mechanisms mentioned above, but any other chemical mechanism could have been used, since the subroutine is a general one. A list of the reactions included in the different chemical mechanisms is presented in Appendix C.
Several preliminary tests have been done to ensure mass conservation of an inert species emitted by the stack or incorporated later on to the plume. All the tests have produced good results, observing only small differences in the total mass of a given species between the initial and ending point, due to accumulation of small numerical inaccuracies. As an example, when a stack emits 5 g s1 of SF6, with wind speed of 5 m s1, considering a cell length of 10 m, the sum of the SF6 masses over all the sub cells in the first cell was 10.0031 g, and in the cell number 10,000 (100 km) was 10.0002 g. Another test performed to verify the capability of the model is to simulate dispersion of non-reactive pollutants. We have done simulations comparing the CAPAS model and the ISCST3 (US EPA, 1995a,b) model. Different simulations were carried out covering a wide range of atmospheric conditions (stability, wind speed, and height of the inversion layer) and physical characteristics of the stack (height, diameter, temperature, and exit velocity of the gases). The numerical results are similar in all the cases for both models, showing that CAPAS is capable to predict not only the chemical reactions that occur in the atmosphere but also the dispersion of the pollutants. For example, with conditions described in Table 1, the average difference between both models in the receptors placed each 500 m Table 1 Ambient conditions and stack parameters for comparison with ISCST3 model Ambient temperature Wind speed Stability class Inversion layer height
300 K 2 m s1 C 5000 m
Stack height Stack diameter Gases exit speed Gases exit temperature Emission rate
10 m 2m 14 m s1 400 K 1 g s1
1244
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
Tracer concentration [µg m-3]
3.0 ISCST3 CAPAS
2.5
2.0
1.5
1.0
0.5
0.0
0
5000
10000
15000
Table 2 Results of comparison with Indianapolis dataset X distance [km]
Y distance [km]
Measured value [ppt]
Calculated value [ppt]
Percentage of difference
1.46 1.93 1.25 2.94 3.90 0.95 0.86 1.91 2.94 0.70
0.17 0.36 0.38 0.35 0.66 0.2 0.31 0.11 0.75 0.26
304 264 211 187 182 180 154 151 145 128
360 245 240 159 182 156 157 153 165 123
18 7 14 15 0 13 2 1 14 4
20000
Distance from stack [m] Fig. 2. Comparison between the ISCST3 and the CAPAS models to simulate the dispersion of non-reactive pollutants.
under the center of the plume is (11G9)%, with a maximum of 30% (Fig. 2).
4. Results and comparison with experimental data Comparison with experimental data can be divided in two parts: Comparison with Indianapolis experiment (Olesen, 1994) to verify the dispersion of non reactive pollutants and a comparison with Stewart and Liu (1981) work, to evaluate the chemical transformations inside the plume. First, the comparison was done against the EPRI Indianapolis field study that involved SF6 tracer releases from the Perry K power plant in Indianapolis, Indiana, USA (latitude 39.8N and longitude 86.2W). A total of 170 h of tracer data are available from September and October, 1985, and represent all stability classes and most wind speed ranges. Meteorological data were taken from the Standard National Weather Source (NWS) and local observations. Concentrations were observed on a network of about 160 ground-level monitors on arcs at distances ranging from 0.25 to 12.0 km from the stack. Although the measurements were done in an urban ambience, wind tunnel study have suggested that the buildings do not influence the plume, which tended to rise one hundred meters or more above the stack top most of the time (Hanna and Chang, 1991). As meteorological conditions were not constant over the integration time (1 h), we have chosen one hour with similar conditions than the previous and following ones. The comparison was done for day 05/10/85 at 10 AM. In Table 2 is listed the location of the ten receptors which have higher values of measured SF6 concentration, along with the concentrations calculated by the CAPAS model and the percentage of difference. The largest difference (the model overestimated the SF6 concentration by 18%) is found in the receptor with the
highest concentration. The average difference in these ten receptors was (9G7) %. Tests done for other days of the same dataset have shown similar results regarding to the average differences. General performance of the CAPAS model can be considered good, in terms of dispersion of inert species. The second type of tests performed had the purpose to analyze the capability of the model to simulate reactive plumes. The publicly available database of industrial emissions and air quality data is very small. Some data are presented in the pioneer work by Stewart and Liu (1981), in their development of the RPM. In this work, they show their model results compared to experimental data collected by an airplane flying transversally to the plume at several points of two power plants: the Widows Creek Power Plant and the Oak Creek Power Plant. The most important data of the set of experimental conditions that can be extracted from the work are presented in Table 3. The reported results are the O3 average concentrations measured in the plume during each flight and the calculated RPM values. The CAPAS model gives as output concentrations at an array of virtual receptors; therefore the plume average was calculated by taking an average of the concentration calculated at the receptors located inside the plume, placed at 100 m each. There are two factors that turn difficult the comparison between our model and the experimental data: one is related to the fact that we do not know the emission rate of the stack and the second one is that we have no information concerning the starting and ending point of
Table 3 Atmospheric conditions during the experiment for the Widows Creek Power Plant (values extracted from Stewart and Liu, 1981) Date: August 23, 1978 Time: 04:00–07:00 am Stability class: B Wind Speed: 3 to 5 m s1 Atmospheric VOC/NOx ratio: 30:1
1245
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
0.12
0.18
O3 concentration [ppm]
O3 concentration [ppm]
0.20
RPM CAPAS Plume average Plume average (CAPAS)
0.14
0.10 0.08 0.06 0.04
0.16
RPM CAPAS Experimental data
0.14 0.12 0.10 0.08 0.06 0.04
0.02 0.02 0.00
0
20
40
60
80
100
120
Downwind distance [km] Fig. 3. Measured (Widows Creek Power Plant plume) and simulated (CAPAS and RPM) ozone concentration as a function of the downwind distance. The measured and RPM results have been adapted from Stewart and Liu (1981).
the plume. Hence, the results presented here can account for general trends but not absolute values. In Figs. 3–5 we show the results obtained with the CAPAS model compared with the RPM model and the experimental data. Tests varying the absorption coefficient (see Section 2.8 above) have been done, and better reproducibility of the experimental data was obtained when the absorption coefficient was set equal to zero for all the species.
5. Conclusions We have presented here a new model to estimate reactive gaseous emissions from industrial stacks. A
O3 concentration [ppm]
0.14 0.12
RPM CAPAS Experimental data
0.10 0.08 0.06 0.04 0.02 0.00 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 Lateral distance at 32 km from source [km]
Fig. 4. Measured (Oak Creek Power Plant plume) and simulated (CAPAS and RPM) ozone concentration as a function of the lateral distance looking at 32 km from the stack. The measured and RPM results have been adapted from Stewart and Liu (1981).
0.00 -4
-3
-2
-1
0
1
2
3
4
Lateral distance at 125 km from source [km]
Fig. 5. Measured (Oak Creek Power Plant plume) and simulated (CAPAS and RPM) ozone concentration as a function of the lateral distance looking at 125 km from the stack. The measured and RPM results have been adapted from Stewart and Liu (1981).
complete validation of the model has not been done due to the lack of experimental data, but the results are promising. The model is available free of charge to anyone that wishes to predict the impact of reactive emissions and compare the performance of CAPAS against other models. As there are very few experimental data, no conclusion about which model (RPM-IV, ROME, ADMS3 or CAPAS) reproduces best the reality and all of them must be considered alternatives and reevaluated when complete datasets of measurements in different kind of reactive plumes are available. The assessment also should include consideration of practical aspects of the models including cost, computing requirements, compatibility with meteorological and topographical data, range of applications, and speed of operation and versatility of outputs. Although there is not a clear reason to choose this model over the commonly used RPM-IV, we can say that the main differences are the CAPAS powerful numerical solver to handle stiff systems, the dissimilar geometry of the cell in the plume and the lesser operator requirements to run the CAPAS model. There are differences between the outputs of the two models but these do not follow any consistent pattern. It is not practicable to make generalized recommendations on the use of one or the other model. The choice of a model should be considered for each specific application depending on the particular situation and on the merits of different models in that situation. An important limitation of CAPAS is its inability to treat complex terrain; only smooth terrain can be simulated. This disadvantage is common to many dispersion-only models, and will be corrected in a future version of CAPAS by dividing the sub cells into minor
1246
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
parcels. This modification could also produce better agreement with the experimental measurements, by treating differently the gases coming from another area or line sources (freeways, garbage depots, etc.) or from others point sources like different industries. The CAPAS model was tested for a range of conditions to ensure that model results were physically and chemically consistent. The model will be evaluated against existing data in future work as they become available. A LIDAR instrument would be of great importance to develop a three dimensional map of pollution distribution. Knowing all the physical characteristics of the stack as well as its emission rate, the meteorological conditions, the ambient concentrations of the principal species and the concentrations in the plume can be measured. In the case that the rate of emission from the stack and the meteorological conditions could be kept constant during a given period of time (a couple of hours), the whole plume can be scanned at several distances, giving the chance to make fine adjustments to the model.
The result of the integral is
Volume Z
psya sza d23 psyb sza d22 psya szb d22 C C Cpsyb szb d2 3 2 2 psya sza d13 psyb sza d12 psya szb d12 3 2 2 psyb szb d1 ða1:4Þ
A.2. Case (b) The volume is calculated according to the integral that represents the volume of the sub cell over a given height h Z Volume Z
psy x sz x C2hpðxÞ d1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sz ðxÞ 2 2 pðxÞ sy ðxÞ pðxÞ sy x pðxÞ sy x sz ðxÞarcsin dx sy ðxÞ
Acknowledgements
where
We thank CONICET, SeCyT (UNC), Fundacio´n Antorchas, and the Agencia Nacional de Promocio´n Cientı´ fica y Tecnolo´gica (Pre´stamo BID 1201/OC-AR N PICT 06-06358) for partial support of the work reported here.
pðxÞZsy x
d2
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi h 1 sz ðxÞ
ða1:5Þ
ða1:6Þ
The analytical result of the previous integral is 2
Appendix A
a1 Zðsza d2 Cszb Þ h2
To calculate the volume of an elliptic truncated cone between distances d1 and d2, we define first four quantities
a2 Zðsza dCszb Þ h2
ða1:8Þ
a3 ZðhCszb Csza d1 Þð hCszb Csza d1 Þ
ða1:9Þ
s endy;j s starty;j sya Z syb Zs endy;j sya d2 length s endz;j s startz;j szb Zs endz;j sza d2 sza Z length
ða1:1Þ
2
ða1:7Þ
ða1:2Þ
A.1. Case (a) The volume is calculated according to the integral that represents the volume of the sub cell Z
d2
VolumeZ d1
psy ðxÞsz ðxÞdx
ða1:3Þ
s2 Z1=6
ða1:10Þ
pffiffiffiffiffi s5 Z2psya s3za d23 3arcsin 1= szb Csza d2 a1 d22 syb s3za pffiffiffiffiffi C3pd22 syb s3za 3arcsin 1= szb Csza d2 a1 !d22 sya szb s2za C3pd22 sya szb s2za 2 pffiffiffiffiffi 2 a1 d2 syb szb sza 6arcsin 1= szb Csza pffiffiffiffiffi 2 C6pd2 syb szb sza C4hsya a1 sza d2
ða1:11Þ
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
pffiffiffiffiffi s4 Z s5 5hsya a1 szb 2sya h3 log sza d2 Cszb pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C ðhCsza d2 Cszb Þð hCsza d2 Cszb Þ pffiffiffiffiffi pffiffiffiffiffi C 9h a1 syb sza C3arctgðh= a1 Þszb 3 sya pffiffiffiffiffi C 6h2 arctgðh= a1 Þsyb sza pffiffiffiffiffi 6h2 arctgðh= a1 Þsya szb pffiffiffiffiffi 2arcsinð1=ðszb Csza d2 ÞÞ a1 sya s3za d32
Fb Zgvs d2s
Fm Zv2s d2s
DT 4Ts
ða2:3Þ
Based on Pasquill stability, different formulae are used: ða1:12Þ
s5 Z1=s2za
ða1:13Þ
s3 Zs4 s5
ða1:14Þ
s1 Zs2 s3
ða1:15Þ ða1:16Þ
pffiffiffiffiffi pffiffiffiffiffi s5 Z s6 5hsya a2 szb 2sya h3 log sza d1 Cszb C a3 Þ pffiffiffiffiffi pffiffiffiffiffi C9h a2 syb sza C3 arctg h= a2 s2zb syb sza pffiffiffiffiffi pffiffiffiffiffi arctg h= a2 s3zb sya C6h2 arctg h= a2 syb sza pffiffiffiffiffi 6h2 arctg h= a2 sya szb pffiffiffiffiffi 2 arcsin 1=ðszb Csza d1 Þ a2 sya s3za d31 ða1:17Þ s6 Z1=s2za
ða1:18Þ
s4 Zs5 s6
ða1:19Þ
s2 Zs3 s4
ða1:20Þ
B.1. Unstable – neutral atmosphere
ða1:21Þ
Appendix B First, the modified height of the stack (h#s ) is calculated depending on the wind speed (ws). i h h#s Zhs C2ds wuss 1:5 ws !1:5us ða2:1Þ h#s Zhs ws R1:5us The buoyancy (Fb) and momentum (Fm) fluxes must be calculated. DT is the difference between the gases exit temperature (Ts) and the ambience temperature (Ta)
vs1=3 2=3
ds vs 2=3 Fb R55 ðDTÞc Z0:00575 Ts 1=3 ds
ða2:4Þ
If the difference between Ts and DT is greater or equal than (DT )c, the plume is dominated by heat flux, else by mechanical turbulence. B.1.1. Heat flux The final height (he) of the plume and the distance xf are given by 3=4 F 5=8 Fb !55 he Zh#s C21:425 b xf Z49Fb us ða2:5Þ 3=5 F 2=8 Fb R55 he Zh#s C38:71 b xf Z119ZFb us B.1.2. Mechanical turbulence The final height (he) of the plume and the distance (xf) to which it takes that value are given by he Zh#s C3ds
vs 2=5 xf Z119Fb us
ða2:6Þ
B.2. Stable atmosphere The stability parameter s is defined as follows: sZg
VolumeZs1 Cs2
ða2:2Þ
Ta 4Ts
Fb !55 ðDTÞc Z0:0297Ts
s3 Z 1=6
1247
vq=vz Ta
ða2:7Þ
where vq/vz is the change of potential temperature with height of the atmosphere. The value is 0.020 K m1 for E stability and 0.035 K m1 for F stability. The same criterion used in unstable and neutral atmospheres is used here, but (DT )c is calculated following: pffiffi ða2:8Þ ðDTÞc Z0:019582Ts vs s B.2.1. Heat flux The final height (he) of the plume and the distance (xf) to which it takes that value are given by 1=3 Fb us he Zh#s C2:6 xf Z2:0715 pffiffi ða2:9Þ us s s
1248
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
B.2.2. Mechanical turbulence Two different equations are computed, and the one used is the one that gives the greater final height: 1=3 Fm vs # pffiffi he Zhs C1:5 he Zh#s C3ds ða2:10Þ us us s The distance to the final height is: pus xf Z pffiffi 2 s
ða2:11Þ
B.3. Smaller distances than xf
Seinfeld mechanism (Seinfeld, 1986) Reactants / Products
Rate constant
NO2 / NO C O O C O2 C M / O3 C M NO C O3 / NO2 C O2 RH C OH C O2 / RO2 C H2O RCHO C OH C O2 / RC(O)O2 C H2O RCHO C 2 O2 / RO2 C HO2 C CO HO2 C NO / NO2 C OH RO2 C NO C O2 / NO2 C RCHO C HO2 RC(O)O2 C NO C O2 / NO2 C RO2 C CO2 OH C NO2 / HNO3 RC(O)O2 C NO2 / PAN PAN / RC(O)O2 C NO2
8.88 ! 103b 6.00 ! 1034c 1.80 ! 1014a 1.04 ! 1025c 6.44 ! 1025c 5.25 ! 1033c 8.21 ! 1012a 3.10 ! 1025c 3.10 ! 1025c 1.09 ! 1011a 4.66 ! 1012a 3.57 ! 104b
a
The coefficient of mechanical turbulence penetration (bj), must be calculated. 1 us bj Z C 3 vs
c
[cm3 molec1 s1]. [s1]. [cm6 molec2 s1].
ða2:12Þ
B.3.1. Unstable – neutral atmosphere The effective height for a given distance (x) is:
he Zh#s C
b
3Fm x b2j u2s
!1=3 ða2:13Þ
B.3.2. Stable atmosphere The effective height for a given distance (x) is: " # pffiffi sinðx s =u Þ s he Zh#s C 3Fm ða2:14Þ pffiffi b2j us s
CBM-IV Mechanism (Gery et al., 1989) Reactants / Products
Rate constant
NO2 / NO C O O / O3 O3 C NO / NO2
Photolysis 1b 8.383 ! 104 exp(1175/temp)b 2.643 ! 103 exp(1370/ temp)a 1.375 ! 104a 2.303 ! 102 exp(687/temp)a 3.233 ! 102 exp(602/temp)a 1.760 ! 102 exp(2450/ temp)a Photolysis 2b Photolysis 3b 1.147 exp(390/temp)b 3.260a 2.344 ! 103 exp(940/temp)a 2.100 ! 101 exp(580/temp)a Photolysis 4b 1.909 ! 104 exp(250/temp)a 3.660 ! 101 exp(1230/ temp)a 7.849 ! 102 exp(256/temp)a 1.900 ! 106a 2.110 ! 1016 exp(10897/ temp)b 2.600 ! 105 exp(530/temp)a 1.600 ! 1011c 6.554 ! 102 exp(806/temp)a Photolysis 5b 9.770 ! 103a 1.500 ! 105a 1.537 ! 103 exp(713/temp)a 7.600 exp(1000/temp)a 5.482 ! 103 exp(240/temp)a 1.640 ! 102 exp(749/temp)a 2.876 ! 1015 exp(10121/emp)b 1.909 ! 103 exp(380/temp)a 8.739 ! 101 exp(1150/temp)a 7.690 ! 1010 exp(5800/temp)c Photolysis 6b
O C NO2 / NO O C NO2 / NO3 O C NO / NO2 O3 C NO2 / NO3
In both cases, the value of the effective height of the plume must be compared with the previously calculated, and the bigger one must be used as the effective height.
O3 / O O3 / O1D O1D / O H2O C O1D / 2 OH O3 C OH / HO2 O3 C HO2 / OH NO3 / NO2 C O NO3 C NO / 2 NO2 NO3 C NO2 / NO C NO2
Appendix C
NO3 C NO2 / N2O5 N2O5 C H2O / 2 HNO3 N2O5 / NO3 C NO2
Høv mechanism (Høv et al., 1989) Reactants / Products
Rate constant
HC C OH / 4 RO2 C 2 ALD ALD / 2 HO2 C CO RO2 C NO / NO2 C ALD C HO2 NO2 / NO C O3 NO C O3 / NO2 C O2 O3 / O2 C O1D O1D C H2O / 2 OH NO2 C OH / HNO3 CO C OH / CO2 C HO2
6 ! 1012a 7.8 ! 105 exp(0.87/cos q)b 8 ! 1012a 1 ! 102 exp(0.39/cos q)b 1.6 ! 1014a 1.9 ! 104 exp(1.9/cos q)b 2.3 ! 1010a 1 ! 1011a 2.9 ! 1013a
a b
[cm3 molec1 s1]. [s1].
NO C NO / 2 NO2 NO C NO2 C H2O / 2 HONO OH C NO / HONO HONO / OH C NO OH C HONO / NO2 HONO C HONO / NO C NO2 OH C NO2 / HNO3 OH C HNO3 / NO3 HO2 C NO / OH C NO2 HO2 C NO2 / PNA PNA / HO2 C NO2 OH C PNA / NO2 HO2 C HO2 / H2O2 HO2 C HO2 C H2O / H2O2 H2O2 / 2 OH
1249
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250 Appendix (continued)
Appendix (continued)
Reactants / Products
Rate constant
Reactants / Products
Rate constant
OH C H2O2 / HO2 OH C CO / HO2 FORM C OH / HO2 C CO FORM / 2 HO2 C CO FORM / CO FORM C O / OH C HO2 C CO FORM C NO3 / HNO3 C HO2 C CO ALD2 C O / C2O3 C OH ALD2 C OH / C2O3 ALD2 C NO3 / C2O3 C HNO3 ALD2 / XO2 C 2 HO2 C CO C FORM C2O3 C NO / NO2 C XO2 C FORM C HO2 C2O3 C NO2 / PAN
4.720 ! 103 exp(187/temp)a 3.220 ! 102a 1.500 ! 104a Photolysis 7b Photolysis 8b 4.302 ! 104 exp(1550/ temp)a 9.300 ! 101a
OH C CRES / 0.4 CRO C 0.6 XO2 C 0.6 HO2 C 0.3 OPEN NO3 C CRES / CRO C HNO3 CRO C NO2 / OPEN / C2O3 C HO2 C CO OPEN C OH / XO2 C 2 CO C 2 HO2 C C2O3 C FORM OPEN C O3 / 0.03 ALD2 C 0.62 C2O3 C 0.7 FORM C 0.03 XO2 C 0.69 CO C 0.08 OH C 0.76 HO2 C 0.2 GLY OH C XYL / 0.7 HO2 C 0.5 XO2 C 0.2 CRES C 0.8 MGLY C 1.1 PAR C 0.3 TO2 OH C MGLY / XO2 C C2O3 MGLY / C2O3 C HO2 C CO O C ISOP / 0.6 HO2 C 0.8 ALD2 C 0.55 OLE C 0.5 XO2 C 0.5 CO C 0.45 ETH C 0.90 PAR OH C ISOP / XO2 C FORM C 0.67 HO2 C 0.13 XO2N C ETH C 0.2 ALD2 C 0.4 MGLY C 0.2 C2O3 O3 C ISOP / FORM C 0.4 ALD2 C 0.55 ETH C 0.2 MGLY C 0.1 PAR C 0.06 CO C 0.44 HO2 C 0.1 OH NO3 C ISOP / ! O2N XO2 C NO / NO2 XO2 C XO2 / XO2N C NO / NR / NR APIN / 7 PAR C 0.75 ALD2 C 0.75 OLE UNKN / 7.5 PAR C 0.875 ALD2 C 0.375 OLE NO3 / NO
6.100 ! 104a
PAN / C2O3 C NO2 C2O3 C C2O3 / 2 ! O2 C 2 FORM C 2 HO2 C2O3 C HO2 / 0.79 FORM C 0.79 XO2 C 0.79 HO2 C 0.79 OH OH / XO2 C FORM C HO2 PAR C OH / 0.87 XO2 C 0.13 XO2N C 0.11 HO2 C 0.11 ALD2 C 0.76 ROR – 0.11 PAR ROR / 1.1 ALD2 C 0.96 XO2 C 0.94 HO2 – 2.10 PAR C 0.04 XO2N C 0.02 ROR ROR / HO2 ROR C NO2 / O C OLE / 0.63 ALD2 C 0.38 HO2 C 0.28 ! O2 C 0.2 OH C 0.3 CO C 0.2 FORM C 0.02 XO2N C 0.22 PAR OH C OLE / FORM C ALD2 C XO2 C HO2 – 1 PAR O3 C OLE / 0.5 ALD2 C 0.74 FORM C 0.33 CO C 0.44 HO2 C 0.22 XO2 C 0.1 OH – 1 PAR NO3 C OLE / 0.91 XO2 C 0.09 XO2N C FORM C ALD2 – 1 PAR C NO2 O C ETH / FORM C 0.7 XO2 C CO C 1.7 HO2 C 0.3 OH OH C ETH / XO2 C 1.56 FORM C HO2 C 0.22 ALD2 O3 C ETH / FORM C 0.42 CO C 0.12 HO2 OH C TOL / 0.08 XO2 C 0.36 CRES C 0.44 HO2 C 0.56 TO2 TO2 C NO / 0.9 NO2 C 0.9 HO2 C 0.9 OPEN TO2 / CRES C HO2
1.739 ! 104 exp(986/temp)a 1.037 ! 104 exp(250/temp)a 3.700a Photolysis 9b 7.915 ! 103 exp(250/temp)a 1.180 ! 104 exp(5500/ temp)a 5.616 ! 1018 exp(14000/ temp)b 3.700 ! 103a 9.600 ! 103a
7.774 ! 103 exp(1815/ temp)b 1.203 ! 103a
6.250 ! 1016 exp(8000/ temp)b 9.545 ! 104b 2.200 ! 104a 1.756 ! 104 exp(324/temp)a
7.740 ! 103 exp(504/ temp)a 2.104 ! 101 exp(2105/ temp)a
a b c
3.250 ! 104a 2.000 ! 104 Photolysis 10b 4.400 ! 104a 8.030 ! 102 exp(500/ temp)a
2.453 ! 104 exp(116/temp)a
2.600 ! 104a Photolysis 11b 2.700 ! 104a
1.420 ! 105a
1.800 ! 102a
4.700 ! 102a 1.200 ! 104a 2.550 ! 101 exp(1300/temp)a 1.000 ! 103a 1.000b 1.000 ! 1010b 1.000 ! 1010b Photolysis 12b
[cm3 molec1 s1]. [s1]. [cm6 molec2 s1].
1.135 ! 101a
References
1.540 ! 104 exp(729/temp)a
Bader, G., Deufhlhard, P., 1983. Numerische Mathematik 47, 373. Briggs, G.A., 1969. Plume rise, U.S. Atomic Energy Commission Critical Review Series T/D 25075. Briggs, G.A., 1971. Some recent analyses of plume rise observations. In: Englund, H.M., Beery, W.T. (Eds.), Proceedings of the Second International Clean Air Congress. Academic Press, New York, pp. 1029–1032. Briggs, G.A., 1974. Diffusion estimation for small emissions, Environmental research laboratories air resources atmospheric turbulence and diffusion laboratory 1973 annual report, USAEC Rep ATDL-106 Natl. Oceanic Atmos. Admin., Washington, D.C. Canepa, E., Ratto, C.F., 2003. Algorithms to simulate the transport of pollutant elements: a model validation exercise and sensitivity analysis. Environmental Modelling and Software 18, 365–372.
3.000 ! 103 exp(411/temp)a
1.856 ! 101 exp(2633/ temp)a 3.106 ! 103 exp(322/temp)a
1.200 ! 104a 2.500 ! 102b
1250
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
Carruthers, D.J., Holroyd, R.J., Hunt, J.C.R., Weng, W.-S., Robins, A.G., Thomson, D.J., Smith, F.B., 1994. UK-ADMS, a new approach to modelling dispersion in the earth’s atmospheric boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 52, 139–153. Duff, I.S., Reid, J.K., 1984. Subroutine to Solve a General Sparse System of Linear Equations. Rutherford Appleton Lab. Gabruk, R.S., Sykes, R.I., Seigneur, C., Pai, P., Gillespie, P., Bergstrom, R.W., Saxena, P., 1999. Evaluation of the reactive and optics model of emissions (ROME). Atmospheric Environment 33, 383–399. Gear, C.W., 1971. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, New Jersey, USA. Gery, M.W., Whitten, G.Z., Killus, J.P., Dodge, M.C., 1989. A photochemical kinetics mechanism for urban and regional scale computer modeling. Journal of Geophysical Research 94, 12945– 12956. Goyal, P., Sidhartha, 2004. Modeling and monitoring of suspended particulate matter from Badarpur thermal power station, Delhi. Environmental Modelling and Software 19, 383–390. Hanna, S.R., Chang, J.C., 1991. Modification of HPDM for Urban Conditions and Its Evaluation using the Indianapolis Data Set. Final Report Prepared for EPRI by EARTH TECH, 196 Baker Ave., Concord. Hov, Ø., Zlatev, Z., Berkowicz, R., Eliassen, A., Prahm, L.P., 1989. Comparison of numerical techniques for use in air pollution models with non-linear chemical reactions. Atmospheric Environment 23, 967–983. Madronich, S., 1993. UV radiation in the natural and perturbed atmosphere. In: Tevini, M. (Ed.), Environmental Effects of UV Radiation. Lewis Publishers, Boca Raton, pp. 17–69. Morris, E.R., Chang, E.C., Shepard, S.B., Ligocki, M.P., 1992. ‘‘User’s Guide to Version IV of the Reactive Plume Model (RPMIV)’’. (SYSAPP-92/037), Systems Applications International, San Rafael, California. Olcese, L.E., Toselli, B.M., 1998. Fast and reliable numerical methods to simulate complex chemical kinetic mechanisms. International Journal of Chemical Kinetics 30, 349–358. Olesen H.R., 1994. Operational Short-Range Atmospheric Dispersion Models for Environmental Impact Assessments in Europe. Pasquill, F., 1961. The estimation of the dispersion of windborne material. Meteorology Magazine 90, 33–49.
Seinfeld, J.H., 1986. Atmospheric Chemistry and Physics of Air Pollution. John Wiley & Sons, New York. Stewart, D.A., Liu, M.K., 1981. Development and application of a reactive plume model. Atmospheric Environment 15, 2377–2393. Turner, D.B., 1997. The long lifetime of the dispersion methods in Pasquill in U.S. regulatory air modeling. Journal of Applied Meteorology 36, 1016–1020. US EPA, 1995a. User’s Guide for the Industrial Source Complex (ISC3) Dispersion Models Volume I – User Instructions, EPA-454/ b-95-003a. US EPA, 1995b. User’s Guide for the Industrial Source Complex (ISC3) Dispersion Models Volume II – Description of Model Algorithms, EPA-454/b-95-003b. Beatriz M. Toselli Professor Toselli received her Bachelor Degree in Physical Chemistry from The Universidad Nacional of Co´rdoba, Argentina in 1981 and her Ph. D. in Chemistry from the same University in 1986. She has been a Post Doc at The Space Physics Research Laboratory. The University of Michigan, Ann Arbor, Michigan, U.S.A. from 1987–1990 and a Research Scientist in the Laser Photochemistry Group, Atomic Energy of Canada Limited, Chalk River Laboratories, Canada from 1990–1991. Currently she is an Associate Chemistry Professor at the Department of Physical Chemistry, Universidad Nacional de Co´rdoba and an Independent Research Scientist from the National Research Council of Argentina (CONICET). She is the director of the Air Pollution Group of the department and her research interest are on energy transfer, laser spectroscopy chemistry and air pollution and radiative processes in the atmosphere.
Luis E. Olcese Dr. Olcese received his Bachelor Degree in Chemistry from The Universidad Nacional of Co´rdoba, Argentina in 1994 and his Ph. D. in Chemistry from the same University in 2000. Currently he is a Teaching Assistant at the Department of Physical Chemistry, Universidad Nacional de Co´rdoba and an Assistant Research Scientist from the National Research Council of Argentina (CONICET). He is one of the researchers of the Air Pollution Group and his research interests are atmospheric chemistry, mathematical models and air pollution models development.
Development of a model for reactive emissions from industrial stacks Luis E. Olcese, Beatriz M. Toselli* Departamento de Fı´sico Quı´mica/INFIQC, Facultad de Ciencias Quı´micas, Universidad Nacional de Co´rdoba, 5000 Co´rdoba, Argentina Received 1 May 2003; received in revised form 17 June 2004; accepted 18 August 2004
Abstract We have developed a model, CAPAS, capable of estimating short-term concentrations of primary and secondary pollutants resulting from point source emissions. The model is designed to simulate the complex interaction of plume dispersion and non-linear chemistry. The main features of the model are: (1) the horizontal and vertical resolution within the plume, which offers a realistic treatment of the entrainment process, (2) its flexibility with regard to choices of chemical kinetic mechanism and the inclusion of a powerful numerical solver to handle stiff systems, and (3) its simplicity of operation and adaptability to solve new problems. The realism of the plume model simulations was tested by comparing model calculations of plume concentrations with data of SF6 tracer concentrations and ozone concentrations. The results of the model compare well with both sets of field measurements. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Dispersion models; Plume reactivity; Air pollutants; Non-linear chemistry; Stiff solvers
Software availability Program title: CAPAS Developer: L.E. Olcese Contact address: Departamento de Fı´ sico Quı´ mica/ INFIQC, Facultad de Ciencias Quı´ micas, Universidad Nacional de Co´rdoba, 5000 Co´rdoba, Argentina. E-mail: lolcese@ fcq.unc.edu.ar. Tel.: C54-351-433-4169; fax: C54-351-433-4188. First year available: 2000 Hardware requirements: PC Software requirements: FORTRAN 90 compiler * Corresponding author. Fax: C54 351 433 4188. E-mail address: [email protected] (B.M. Toselli). 1364-8152/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2004.08.008
Program language: Availability:
FORTRAN free of charge, upon request to the authors
Abbreviations and symbols massi,j,k:
mass stacki: cell length: sy,j: s starty,j: s endy,j: sz,j:
Mass of the specie i in the sub cell j in the cell k (If k is not present, it means that it is applied to any cell) Mass of the specie i emitted by the stack Length of every cell Pasquill sigma for horizontal dispersion at sub cell j Pasquill sigma for horizontal dispersion at the beginning of the sub cell j Pasquill sigma for horizontal dispersion at the end of the sub cell j Pasquill sigma for vertical dispersion at sub cell j
1240
s startz,j: s endz,j: absi: volj,k:
vol under surfacej: exchangej:
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
Pasquill sigma for vertical dispersion at the beginning of the sub cell j Pasquill sigma for vertical dispersion at the end of the sub cell j Absorption coefficient in the surface for specie i Volume of the sub cell j in the cell k (If k is not present, it means that it is applied to any cell) Volume of the sub cell j that is located under the surface Volume exchanged between the sub cells j and jC1
1. Introduction Air quality models can be classified into two broad groups: one that includes an extensive treatment of the chemical reactions, mainly oriented to predict formation of secondary pollutants; and another group oriented to simulate dispersion of primary pollutants in different situations (street canyons, industrial plumes, etc.). Usually, chemical industries, power plants, refineries, etc., use plume models like ISCST3 (US EPA, 1995a,b) to estimate the air pollution impact of their production processes. ISCST3 is perhaps the most widely used model for estimating near-field concentrations of non-reactive pollutants. It is a steady-state Gaussian plume model, i.e. dispersion parameters such as meteorological conditions and emission rate are assumed constant while the plume is in route from the source to the receptor. Some recent applications based on Gaussian models can be found in the works of Goyal and Sidharta (2004) and Canepa and Ratto (2003). The main drawback of these models is the inability to handle formation of secondary pollutants like ozone, or to handle the rate of disappearance of the emitted specie when the rate is different from a first order decay. This assumption is valid for CO or SO2, but is erroneous for species like NOX or VOC in urban atmospheres, or in polluted rural areas. In these cases, the formation of secondary pollutants becomes important, and the use of another type of models is necessary. Reactive plume models are often used to estimate the local or short- to medium-range (i.e. up to a few hundred km) impacts of power plants or smelters on air quality. Issues of interest typically include ozone and particulate matter concentrations above the ambient air quality standards of a region, visibility degradation and acid deposition. These models use detailed gas-phase mechanisms to describe the chemical transformations. The particular case of simulating the dispersion of reactive species emitted by an industrial stack has been treated only in few models, mainly the Reactive Plume Model (RPM) (Stewart and Liu, 1981). The RPM
version IV (Morris et al., 1992) is a standard regulatory model used for estimating pollutant concentrations in the atmosphere, resulting from the emissions from point sources such as industrial stacks. It uses either point source emission estimates or initial plume concentrations as inputs, and calculates downwind concentrations, as the plume expands. RPM-IV makes use of a diffusion/reaction equation similar in form to the atmospheric diffusion equation, but with dispersion parameters that evolve with time. This model considers a control volume that moves downwind along the plume trajectory, and solves the diffusion/reaction equation by considering the nonlinear photochemistry and diffusion of chemical species from one cell into another cell. Other existing models that are useful to simulate dispersion of reactive plumes are ROME (Gabruk et al., 1999) and ADMS3 (Carruthers et al., 1994), but this one is not freely available. In this work we present a new model (CAPAS), capable of simulating reactive plumes, easy to use by an operator with little training and easily adaptable to solve new problems or to use new chemical mechanisms. The model code is written in standard FORTRAN 90 and it has a modular structure, thus the replacement of a particular module (i.e., the numerical integrator of the ordinary differential equations) by another one is straightforward. CAPAS model shares some characteristics with the RPM model, but includes several new ones, as it will be shown below.
2. Model formulation This model is based on the Lagrangian approach to the turbulent diffusion (Seinfeld, 1986). We have chosen this approach, and not the Eulerian one, due to the simplicity of the resulting approximations. Briefly, the model consists of an array of cells, each one similar to those used in the box models; the cells are arranged in the direction of the x axis in a continuous way and the shape of each cell is an elliptic truncated cone. The size of a cell is given by four times the Pasquill sigmas in the y and z axes (both different in the beginning and in the end), and user definable in the x axis, typically 5–30 m depending on the wind speed. Each cell is later divided in 10–30 annular and equi-spaced sub cells, with the concentration of the species assumed to be homogeneous in each sub cell. A scheme of the cells and sub-cells can be viewed in Fig. 1. The boundary conditions are defined by an external sub cell with the concentrations of the species taken from the atmospheric concentrations. The steps to set up the model can be summarized as follow: 1. All the gases are injected in the atmosphere without any initial interaction. The gases exiting from the
1241
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
mass of each species is recalculated (Section 2.8). The fraction of the mass of the gases in the region of a sub cell in the previous cell that is occupied by a different sub cell in the new cell is added to the new sub cell. 9. Concentrations of the species are updated depending on the chemical mechanism (Section 2.9). 10. Steps 6–9 are repeated until the final distance is reached.
Several cells
A common problem of many air pollution models is that due to instabilities of the numerical solvers and to the rapid change in the species concentration because of variations in the volume of the cells, some intermediate transient concentrations are negative. We have found that a way to avoid such spurious results is to set up the model such that works internally with the mass of each species, instead of concentration.
Sub cells in a cell
2.1. Initial distribution
Lateral view
Front view
Fig. 1. Schematic of the physical representation of the plume used in the CAPAS model to simulate the chemistry and dispersion of the pollutants.
stack are assumed to be well represented by a Gaussian distribution and located in the first cell; the initial distribution of pollutant mass emitted by the industrial stack is given by the error function in two dimensions (Section 2.1). 2. The gases temperature is the same as the exit temperature of the stack. 3. The final height of the plume is calculated (Section 2.2). 4. The meteorological conditions are assumed constant with time, over the time period of transport from the source to the receptors. 5. The ordinary set of stiff differential equations representing the chemical mechanism is solved to obtain the new values of the species concentration in order to advance the solution. In the case that the chemical mechanism includes photolysis reactions the zenith angle dependence of these reactions is incorporated (Section 2.3). 6. After this first cell, new atmospheric gases are added to the outermost sub cell, and the new temperature is calculated (Sections 2.4 and 2.5). 7. The cell grows in size by the increase in the sy and sz values, and the sub cells grow proportionally (Section 2.6). 8. Depending on the volume of the sub cell (Section 2.7) and the absorption coefficient on the surface, the
A bi-dimensional Gaussian distribution of the species is assumed between the boundaries of the initial cell. The mass of each species in each sub cell is given by the equation: massi;j;1 Z mass stacki ! cell length sy;j
!
pffiffiffi sy;outermost j ! 2 !! sy;j pffiffiffi erf sy;j1 ! 2
! erf
ð1Þ
The outermost ring concentrations are the same as the atmospheric concentrations. The lack of interaction between the plume and the atmosphere in the first cell is justified considering the turbulent nature of the gases coming from the stack, and the short time elapsed between the exit of the gases from the stack and the formation of the plume. 2.2. Final height of the plume The final height of the plume is calculated according to the classical equations proposed by Briggs (1969, 1971, 1974), and revised by Turner (1997). The complete set of equations used is shown in Appendix A. 2.3. Photolysis rate constants When the chemical mechanism includes photolysis reactions, the photolysis rate constants are previously calculated using the TUV 4.1 program (Madronich, 1993), at ten values of zenith angle (0 –10 –20 –30 –40 –50 –60 –70 –78 –86 ). The calculated constant
1242
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
for each photochemical reaction at each angle are entered as input data, and depending on the zenith angle at each cell, the corresponding photolysis rate constant is calculated using a cubic spline interpolation. The computational time needed to calculate the photolysis rate constants is very small, making unnecessary to consider uniform photolysis rate constants along the entire plume. Photolysis rate constants are calculated assuming clear sky conditions, with no clouds and aerosols. However, different cloud cover conditions and aerosol loading can be easily incorporated into the TUV model and therefore included in the CAPAS model.
m1 Z
number of species X
massi;j;k massi;j;k
iZ1
exchangej1 volj
ð4Þ
m2 Z
number of species X
massi;j;kC1 massi;j;k
iZ1
exchangej voljC1
ð5Þ
tempj Z
m1 tempj Cm2 tempjC1 m1Cm2
ð6Þ
2.4. Mass exchange 2.6. Sigma y and sigma z Due to the increase in volume of the plume, when the gases pass from the cell k to the cell kC1, the mass present in each sub cell must be adjusted. The model calculates the volume exchanged between sub cells j and jC1 by calculating the volumes of the sub cells from the innermost until the jth one, minus the same volume but in the previous cell. This can be expressed as: exchangej Z
numberX of sub cells
numberX of sub cells
jZ1
jZ1
voljC1
volj
ð2Þ
Consequently, the mass of a given species in the sub cell j is: massi;j;k Z massi;j;k1 Cmassi;jC1;k1 massi;j;k1
exchangej voljC1;k1
exchangej1 volj;k1
ð3Þ
In the first cells, due to the rapid change in the sub cells volume, the volume exchanged can be greater than the volume of the cell, leading to negative volumes. To overcome this problem, the model selects a new time step and calculates the new volume. If the new calculated volume is greater than zero, is accepted and the calculation continues. If not, the model iteratively decreases the step size until the calculated volumes are positive. 2.5. Temperature variation The temperature of each sub cell is adjusted in every step of the calculation, averaging the temperature of the fraction of the total mass of the sub cell that is not exchanged and the temperature of the fraction that is injected to the sub cell, weighted by the mass. In each cell, the temperature of the outer sub-cell remains constant. For a given sub cell j of a cell k, its temperature is calculated according to:
Both sigmas are calculated according to the procedure described by Pasquill (1961) depending on the atmospheric stability. Urban or rural conditions are used depending on the land use. If different conditions are found along the plume trajectory, the predominant environment is used. 2.7. Sub cell volume Depending on the position of the sub cell relative to the ground, there are three different possibilities to calculate the volume of each sub cell: (a) The sub cell does not reach the surface in any point. (b) The sub cell does reach the surface since its initial x. (c) The sub cell does reach the surface in a point between its initial and final x. The mathematical expression of each volume is obtained by analytical integration. In case (a), the volume is calculated as an annular elliptical truncated cone. In case (b), the volume of the sub cell under zZ0 is subtracted from the calculated in case (a). Case (c) is treated as a combination of the two previous cases, dividing the sub cell in two sections; the first one between the initial distance x and the point where the sub cell does reach the ground, and the second one from there to the final distance x. The volumes of each section are calculated in the same way as cases (a) and (b). The mathematical expressions for each case are included in Appendix B. 2.8. Absorption coefficient The interaction of each species with the surface (cement, grass, sand, etc.) the products formed and the rate constants of the involved reactions (dry and wet deposition, heterogeneous reactions, etc.) are very complex phenomena, and depend on many variables
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
(soil type, temperature, humidity, etc.). As a simplification, three cases are possible for each species: total absorption, partial absorption, and total reflection. To simulate these three situations in this model, we used a coefficient abs, specific for each species. This factor can take values ranging from 1 (the species is completely absorbed in the surface) to 0 (completely reflected). In each sub cell, massi;j Z massi;j massi;j ! absi
vol under surfacej volj
ð7Þ
The experimental difficulty of determining the value of the absorption factor, abs is one of the major handicaps of this model. Nevertheless, the value for the coefficient abs can be fitted to the observed atmospheric concentrations of selected species using iterative methods.
1243
2.10. Concentration at a given receptor The final goal of this air pollution model is to estimate the concentration of selected pollutants at some defined places, called receptors. The concentration of a specie at a given receptor, defined by x, y, z coordinates, is obtained looking at the concentration in the sub cell that included the receptor. Depending on the downwind distance to the stack, the wind speed, and the length of the cell the model calculates inside which cell a given receptor is located. To determine in which sub cell is placed the receptor, the angle formed between the receptor and the center of the plume is calculated together with the distance of the border of each sub cell. If the receptor is located outside the limits of the plume, atmospheric concentrations are used.
2.9. Chemical mechanism 3. Validation One of the most widely used numerical procedures to solve stiff systems of ordinary differential equations (ODEs) is the Gear’s method (Gear, 1971) with automatic time step selection and error control. The method is however; time consuming because it requires inversion of large matrices. This is the method used by RPM-IV model. To improve the performance of the model, the numerical integration subroutine used by CAPAS is based on the Rosenbrock method with Kaps–Rentrop coefficients (Bader and Deufhlhard, 1983), modified to take into account the sparsity of the Jacobian matrix (Duff and Reid, 1984). This (semi-)implicit method for treating complex non-linear kinetic mechanisms has been selected after a detailed study of the different available solvers, and has been chosen due to its stability and velocity (Olcese and Toselli, 1998). We have tested the numerical integrator with different chemical mechanisms proposed by Hov et al. (1989) and Seinfeld (1986), and the CBM-IV mechanism (Gery et al., 1989) obtaining good stability in all of the mechanisms under different initial conditions. For a complete description of the numerical solver and the different tests performed the interested reader should consult the work by Olcese and Toselli (1998). The selection of a given mechanism and the numerical method to solve the stiff ODEs is crucial to the performance of the model. Roughly, 90% of the computational time is spent in the numerical integration of the chemical mechanism. To test the CAPAS model we have chosen for the simulations the three mechanisms mentioned above, but any other chemical mechanism could have been used, since the subroutine is a general one. A list of the reactions included in the different chemical mechanisms is presented in Appendix C.
Several preliminary tests have been done to ensure mass conservation of an inert species emitted by the stack or incorporated later on to the plume. All the tests have produced good results, observing only small differences in the total mass of a given species between the initial and ending point, due to accumulation of small numerical inaccuracies. As an example, when a stack emits 5 g s1 of SF6, with wind speed of 5 m s1, considering a cell length of 10 m, the sum of the SF6 masses over all the sub cells in the first cell was 10.0031 g, and in the cell number 10,000 (100 km) was 10.0002 g. Another test performed to verify the capability of the model is to simulate dispersion of non-reactive pollutants. We have done simulations comparing the CAPAS model and the ISCST3 (US EPA, 1995a,b) model. Different simulations were carried out covering a wide range of atmospheric conditions (stability, wind speed, and height of the inversion layer) and physical characteristics of the stack (height, diameter, temperature, and exit velocity of the gases). The numerical results are similar in all the cases for both models, showing that CAPAS is capable to predict not only the chemical reactions that occur in the atmosphere but also the dispersion of the pollutants. For example, with conditions described in Table 1, the average difference between both models in the receptors placed each 500 m Table 1 Ambient conditions and stack parameters for comparison with ISCST3 model Ambient temperature Wind speed Stability class Inversion layer height
300 K 2 m s1 C 5000 m
Stack height Stack diameter Gases exit speed Gases exit temperature Emission rate
10 m 2m 14 m s1 400 K 1 g s1
1244
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
Tracer concentration [µg m-3]
3.0 ISCST3 CAPAS
2.5
2.0
1.5
1.0
0.5
0.0
0
5000
10000
15000
Table 2 Results of comparison with Indianapolis dataset X distance [km]
Y distance [km]
Measured value [ppt]
Calculated value [ppt]
Percentage of difference
1.46 1.93 1.25 2.94 3.90 0.95 0.86 1.91 2.94 0.70
0.17 0.36 0.38 0.35 0.66 0.2 0.31 0.11 0.75 0.26
304 264 211 187 182 180 154 151 145 128
360 245 240 159 182 156 157 153 165 123
18 7 14 15 0 13 2 1 14 4
20000
Distance from stack [m] Fig. 2. Comparison between the ISCST3 and the CAPAS models to simulate the dispersion of non-reactive pollutants.
under the center of the plume is (11G9)%, with a maximum of 30% (Fig. 2).
4. Results and comparison with experimental data Comparison with experimental data can be divided in two parts: Comparison with Indianapolis experiment (Olesen, 1994) to verify the dispersion of non reactive pollutants and a comparison with Stewart and Liu (1981) work, to evaluate the chemical transformations inside the plume. First, the comparison was done against the EPRI Indianapolis field study that involved SF6 tracer releases from the Perry K power plant in Indianapolis, Indiana, USA (latitude 39.8N and longitude 86.2W). A total of 170 h of tracer data are available from September and October, 1985, and represent all stability classes and most wind speed ranges. Meteorological data were taken from the Standard National Weather Source (NWS) and local observations. Concentrations were observed on a network of about 160 ground-level monitors on arcs at distances ranging from 0.25 to 12.0 km from the stack. Although the measurements were done in an urban ambience, wind tunnel study have suggested that the buildings do not influence the plume, which tended to rise one hundred meters or more above the stack top most of the time (Hanna and Chang, 1991). As meteorological conditions were not constant over the integration time (1 h), we have chosen one hour with similar conditions than the previous and following ones. The comparison was done for day 05/10/85 at 10 AM. In Table 2 is listed the location of the ten receptors which have higher values of measured SF6 concentration, along with the concentrations calculated by the CAPAS model and the percentage of difference. The largest difference (the model overestimated the SF6 concentration by 18%) is found in the receptor with the
highest concentration. The average difference in these ten receptors was (9G7) %. Tests done for other days of the same dataset have shown similar results regarding to the average differences. General performance of the CAPAS model can be considered good, in terms of dispersion of inert species. The second type of tests performed had the purpose to analyze the capability of the model to simulate reactive plumes. The publicly available database of industrial emissions and air quality data is very small. Some data are presented in the pioneer work by Stewart and Liu (1981), in their development of the RPM. In this work, they show their model results compared to experimental data collected by an airplane flying transversally to the plume at several points of two power plants: the Widows Creek Power Plant and the Oak Creek Power Plant. The most important data of the set of experimental conditions that can be extracted from the work are presented in Table 3. The reported results are the O3 average concentrations measured in the plume during each flight and the calculated RPM values. The CAPAS model gives as output concentrations at an array of virtual receptors; therefore the plume average was calculated by taking an average of the concentration calculated at the receptors located inside the plume, placed at 100 m each. There are two factors that turn difficult the comparison between our model and the experimental data: one is related to the fact that we do not know the emission rate of the stack and the second one is that we have no information concerning the starting and ending point of
Table 3 Atmospheric conditions during the experiment for the Widows Creek Power Plant (values extracted from Stewart and Liu, 1981) Date: August 23, 1978 Time: 04:00–07:00 am Stability class: B Wind Speed: 3 to 5 m s1 Atmospheric VOC/NOx ratio: 30:1
1245
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
0.12
0.18
O3 concentration [ppm]
O3 concentration [ppm]
0.20
RPM CAPAS Plume average Plume average (CAPAS)
0.14
0.10 0.08 0.06 0.04
0.16
RPM CAPAS Experimental data
0.14 0.12 0.10 0.08 0.06 0.04
0.02 0.02 0.00
0
20
40
60
80
100
120
Downwind distance [km] Fig. 3. Measured (Widows Creek Power Plant plume) and simulated (CAPAS and RPM) ozone concentration as a function of the downwind distance. The measured and RPM results have been adapted from Stewart and Liu (1981).
the plume. Hence, the results presented here can account for general trends but not absolute values. In Figs. 3–5 we show the results obtained with the CAPAS model compared with the RPM model and the experimental data. Tests varying the absorption coefficient (see Section 2.8 above) have been done, and better reproducibility of the experimental data was obtained when the absorption coefficient was set equal to zero for all the species.
5. Conclusions We have presented here a new model to estimate reactive gaseous emissions from industrial stacks. A
O3 concentration [ppm]
0.14 0.12
RPM CAPAS Experimental data
0.10 0.08 0.06 0.04 0.02 0.00 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 Lateral distance at 32 km from source [km]
Fig. 4. Measured (Oak Creek Power Plant plume) and simulated (CAPAS and RPM) ozone concentration as a function of the lateral distance looking at 32 km from the stack. The measured and RPM results have been adapted from Stewart and Liu (1981).
0.00 -4
-3
-2
-1
0
1
2
3
4
Lateral distance at 125 km from source [km]
Fig. 5. Measured (Oak Creek Power Plant plume) and simulated (CAPAS and RPM) ozone concentration as a function of the lateral distance looking at 125 km from the stack. The measured and RPM results have been adapted from Stewart and Liu (1981).
complete validation of the model has not been done due to the lack of experimental data, but the results are promising. The model is available free of charge to anyone that wishes to predict the impact of reactive emissions and compare the performance of CAPAS against other models. As there are very few experimental data, no conclusion about which model (RPM-IV, ROME, ADMS3 or CAPAS) reproduces best the reality and all of them must be considered alternatives and reevaluated when complete datasets of measurements in different kind of reactive plumes are available. The assessment also should include consideration of practical aspects of the models including cost, computing requirements, compatibility with meteorological and topographical data, range of applications, and speed of operation and versatility of outputs. Although there is not a clear reason to choose this model over the commonly used RPM-IV, we can say that the main differences are the CAPAS powerful numerical solver to handle stiff systems, the dissimilar geometry of the cell in the plume and the lesser operator requirements to run the CAPAS model. There are differences between the outputs of the two models but these do not follow any consistent pattern. It is not practicable to make generalized recommendations on the use of one or the other model. The choice of a model should be considered for each specific application depending on the particular situation and on the merits of different models in that situation. An important limitation of CAPAS is its inability to treat complex terrain; only smooth terrain can be simulated. This disadvantage is common to many dispersion-only models, and will be corrected in a future version of CAPAS by dividing the sub cells into minor
1246
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
parcels. This modification could also produce better agreement with the experimental measurements, by treating differently the gases coming from another area or line sources (freeways, garbage depots, etc.) or from others point sources like different industries. The CAPAS model was tested for a range of conditions to ensure that model results were physically and chemically consistent. The model will be evaluated against existing data in future work as they become available. A LIDAR instrument would be of great importance to develop a three dimensional map of pollution distribution. Knowing all the physical characteristics of the stack as well as its emission rate, the meteorological conditions, the ambient concentrations of the principal species and the concentrations in the plume can be measured. In the case that the rate of emission from the stack and the meteorological conditions could be kept constant during a given period of time (a couple of hours), the whole plume can be scanned at several distances, giving the chance to make fine adjustments to the model.
The result of the integral is
Volume Z
psya sza d23 psyb sza d22 psya szb d22 C C Cpsyb szb d2 3 2 2 psya sza d13 psyb sza d12 psya szb d12 3 2 2 psyb szb d1 ða1:4Þ
A.2. Case (b) The volume is calculated according to the integral that represents the volume of the sub cell over a given height h Z Volume Z
psy x sz x C2hpðxÞ d1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sz ðxÞ 2 2 pðxÞ sy ðxÞ pðxÞ sy x pðxÞ sy x sz ðxÞarcsin dx sy ðxÞ
Acknowledgements
where
We thank CONICET, SeCyT (UNC), Fundacio´n Antorchas, and the Agencia Nacional de Promocio´n Cientı´ fica y Tecnolo´gica (Pre´stamo BID 1201/OC-AR N PICT 06-06358) for partial support of the work reported here.
pðxÞZsy x
d2
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi h 1 sz ðxÞ
ða1:5Þ
ða1:6Þ
The analytical result of the previous integral is 2
Appendix A
a1 Zðsza d2 Cszb Þ h2
To calculate the volume of an elliptic truncated cone between distances d1 and d2, we define first four quantities
a2 Zðsza dCszb Þ h2
ða1:8Þ
a3 ZðhCszb Csza d1 Þð hCszb Csza d1 Þ
ða1:9Þ
s endy;j s starty;j sya Z syb Zs endy;j sya d2 length s endz;j s startz;j szb Zs endz;j sza d2 sza Z length
ða1:1Þ
2
ða1:7Þ
ða1:2Þ
A.1. Case (a) The volume is calculated according to the integral that represents the volume of the sub cell Z
d2
VolumeZ d1
psy ðxÞsz ðxÞdx
ða1:3Þ
s2 Z1=6
ða1:10Þ
pffiffiffiffiffi s5 Z2psya s3za d23 3arcsin 1= szb Csza d2 a1 d22 syb s3za pffiffiffiffiffi C3pd22 syb s3za 3arcsin 1= szb Csza d2 a1 !d22 sya szb s2za C3pd22 sya szb s2za 2 pffiffiffiffiffi 2 a1 d2 syb szb sza 6arcsin 1= szb Csza pffiffiffiffiffi 2 C6pd2 syb szb sza C4hsya a1 sza d2
ða1:11Þ
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
pffiffiffiffiffi s4 Z s5 5hsya a1 szb 2sya h3 log sza d2 Cszb pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C ðhCsza d2 Cszb Þð hCsza d2 Cszb Þ pffiffiffiffiffi pffiffiffiffiffi C 9h a1 syb sza C3arctgðh= a1 Þszb 3 sya pffiffiffiffiffi C 6h2 arctgðh= a1 Þsyb sza pffiffiffiffiffi 6h2 arctgðh= a1 Þsya szb pffiffiffiffiffi 2arcsinð1=ðszb Csza d2 ÞÞ a1 sya s3za d32
Fb Zgvs d2s
Fm Zv2s d2s
DT 4Ts
ða2:3Þ
Based on Pasquill stability, different formulae are used: ða1:12Þ
s5 Z1=s2za
ða1:13Þ
s3 Zs4 s5
ða1:14Þ
s1 Zs2 s3
ða1:15Þ ða1:16Þ
pffiffiffiffiffi pffiffiffiffiffi s5 Z s6 5hsya a2 szb 2sya h3 log sza d1 Cszb C a3 Þ pffiffiffiffiffi pffiffiffiffiffi C9h a2 syb sza C3 arctg h= a2 s2zb syb sza pffiffiffiffiffi pffiffiffiffiffi arctg h= a2 s3zb sya C6h2 arctg h= a2 syb sza pffiffiffiffiffi 6h2 arctg h= a2 sya szb pffiffiffiffiffi 2 arcsin 1=ðszb Csza d1 Þ a2 sya s3za d31 ða1:17Þ s6 Z1=s2za
ða1:18Þ
s4 Zs5 s6
ða1:19Þ
s2 Zs3 s4
ða1:20Þ
B.1. Unstable – neutral atmosphere
ða1:21Þ
Appendix B First, the modified height of the stack (h#s ) is calculated depending on the wind speed (ws). i h h#s Zhs C2ds wuss 1:5 ws !1:5us ða2:1Þ h#s Zhs ws R1:5us The buoyancy (Fb) and momentum (Fm) fluxes must be calculated. DT is the difference between the gases exit temperature (Ts) and the ambience temperature (Ta)
vs1=3 2=3
ds vs 2=3 Fb R55 ðDTÞc Z0:00575 Ts 1=3 ds
ða2:4Þ
If the difference between Ts and DT is greater or equal than (DT )c, the plume is dominated by heat flux, else by mechanical turbulence. B.1.1. Heat flux The final height (he) of the plume and the distance xf are given by 3=4 F 5=8 Fb !55 he Zh#s C21:425 b xf Z49Fb us ða2:5Þ 3=5 F 2=8 Fb R55 he Zh#s C38:71 b xf Z119ZFb us B.1.2. Mechanical turbulence The final height (he) of the plume and the distance (xf) to which it takes that value are given by he Zh#s C3ds
vs 2=5 xf Z119Fb us
ða2:6Þ
B.2. Stable atmosphere The stability parameter s is defined as follows: sZg
VolumeZs1 Cs2
ða2:2Þ
Ta 4Ts
Fb !55 ðDTÞc Z0:0297Ts
s3 Z 1=6
1247
vq=vz Ta
ða2:7Þ
where vq/vz is the change of potential temperature with height of the atmosphere. The value is 0.020 K m1 for E stability and 0.035 K m1 for F stability. The same criterion used in unstable and neutral atmospheres is used here, but (DT )c is calculated following: pffiffi ða2:8Þ ðDTÞc Z0:019582Ts vs s B.2.1. Heat flux The final height (he) of the plume and the distance (xf) to which it takes that value are given by 1=3 Fb us he Zh#s C2:6 xf Z2:0715 pffiffi ða2:9Þ us s s
1248
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
B.2.2. Mechanical turbulence Two different equations are computed, and the one used is the one that gives the greater final height: 1=3 Fm vs # pffiffi he Zhs C1:5 he Zh#s C3ds ða2:10Þ us us s The distance to the final height is: pus xf Z pffiffi 2 s
ða2:11Þ
B.3. Smaller distances than xf
Seinfeld mechanism (Seinfeld, 1986) Reactants / Products
Rate constant
NO2 / NO C O O C O2 C M / O3 C M NO C O3 / NO2 C O2 RH C OH C O2 / RO2 C H2O RCHO C OH C O2 / RC(O)O2 C H2O RCHO C 2 O2 / RO2 C HO2 C CO HO2 C NO / NO2 C OH RO2 C NO C O2 / NO2 C RCHO C HO2 RC(O)O2 C NO C O2 / NO2 C RO2 C CO2 OH C NO2 / HNO3 RC(O)O2 C NO2 / PAN PAN / RC(O)O2 C NO2
8.88 ! 103b 6.00 ! 1034c 1.80 ! 1014a 1.04 ! 1025c 6.44 ! 1025c 5.25 ! 1033c 8.21 ! 1012a 3.10 ! 1025c 3.10 ! 1025c 1.09 ! 1011a 4.66 ! 1012a 3.57 ! 104b
a
The coefficient of mechanical turbulence penetration (bj), must be calculated. 1 us bj Z C 3 vs
c
[cm3 molec1 s1]. [s1]. [cm6 molec2 s1].
ða2:12Þ
B.3.1. Unstable – neutral atmosphere The effective height for a given distance (x) is:
he Zh#s C
b
3Fm x b2j u2s
!1=3 ða2:13Þ
B.3.2. Stable atmosphere The effective height for a given distance (x) is: " # pffiffi sinðx s =u Þ s he Zh#s C 3Fm ða2:14Þ pffiffi b2j us s
CBM-IV Mechanism (Gery et al., 1989) Reactants / Products
Rate constant
NO2 / NO C O O / O3 O3 C NO / NO2
Photolysis 1b 8.383 ! 104 exp(1175/temp)b 2.643 ! 103 exp(1370/ temp)a 1.375 ! 104a 2.303 ! 102 exp(687/temp)a 3.233 ! 102 exp(602/temp)a 1.760 ! 102 exp(2450/ temp)a Photolysis 2b Photolysis 3b 1.147 exp(390/temp)b 3.260a 2.344 ! 103 exp(940/temp)a 2.100 ! 101 exp(580/temp)a Photolysis 4b 1.909 ! 104 exp(250/temp)a 3.660 ! 101 exp(1230/ temp)a 7.849 ! 102 exp(256/temp)a 1.900 ! 106a 2.110 ! 1016 exp(10897/ temp)b 2.600 ! 105 exp(530/temp)a 1.600 ! 1011c 6.554 ! 102 exp(806/temp)a Photolysis 5b 9.770 ! 103a 1.500 ! 105a 1.537 ! 103 exp(713/temp)a 7.600 exp(1000/temp)a 5.482 ! 103 exp(240/temp)a 1.640 ! 102 exp(749/temp)a 2.876 ! 1015 exp(10121/emp)b 1.909 ! 103 exp(380/temp)a 8.739 ! 101 exp(1150/temp)a 7.690 ! 1010 exp(5800/temp)c Photolysis 6b
O C NO2 / NO O C NO2 / NO3 O C NO / NO2 O3 C NO2 / NO3
In both cases, the value of the effective height of the plume must be compared with the previously calculated, and the bigger one must be used as the effective height.
O3 / O O3 / O1D O1D / O H2O C O1D / 2 OH O3 C OH / HO2 O3 C HO2 / OH NO3 / NO2 C O NO3 C NO / 2 NO2 NO3 C NO2 / NO C NO2
Appendix C
NO3 C NO2 / N2O5 N2O5 C H2O / 2 HNO3 N2O5 / NO3 C NO2
Høv mechanism (Høv et al., 1989) Reactants / Products
Rate constant
HC C OH / 4 RO2 C 2 ALD ALD / 2 HO2 C CO RO2 C NO / NO2 C ALD C HO2 NO2 / NO C O3 NO C O3 / NO2 C O2 O3 / O2 C O1D O1D C H2O / 2 OH NO2 C OH / HNO3 CO C OH / CO2 C HO2
6 ! 1012a 7.8 ! 105 exp(0.87/cos q)b 8 ! 1012a 1 ! 102 exp(0.39/cos q)b 1.6 ! 1014a 1.9 ! 104 exp(1.9/cos q)b 2.3 ! 1010a 1 ! 1011a 2.9 ! 1013a
a b
[cm3 molec1 s1]. [s1].
NO C NO / 2 NO2 NO C NO2 C H2O / 2 HONO OH C NO / HONO HONO / OH C NO OH C HONO / NO2 HONO C HONO / NO C NO2 OH C NO2 / HNO3 OH C HNO3 / NO3 HO2 C NO / OH C NO2 HO2 C NO2 / PNA PNA / HO2 C NO2 OH C PNA / NO2 HO2 C HO2 / H2O2 HO2 C HO2 C H2O / H2O2 H2O2 / 2 OH
1249
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250 Appendix (continued)
Appendix (continued)
Reactants / Products
Rate constant
Reactants / Products
Rate constant
OH C H2O2 / HO2 OH C CO / HO2 FORM C OH / HO2 C CO FORM / 2 HO2 C CO FORM / CO FORM C O / OH C HO2 C CO FORM C NO3 / HNO3 C HO2 C CO ALD2 C O / C2O3 C OH ALD2 C OH / C2O3 ALD2 C NO3 / C2O3 C HNO3 ALD2 / XO2 C 2 HO2 C CO C FORM C2O3 C NO / NO2 C XO2 C FORM C HO2 C2O3 C NO2 / PAN
4.720 ! 103 exp(187/temp)a 3.220 ! 102a 1.500 ! 104a Photolysis 7b Photolysis 8b 4.302 ! 104 exp(1550/ temp)a 9.300 ! 101a
OH C CRES / 0.4 CRO C 0.6 XO2 C 0.6 HO2 C 0.3 OPEN NO3 C CRES / CRO C HNO3 CRO C NO2 / OPEN / C2O3 C HO2 C CO OPEN C OH / XO2 C 2 CO C 2 HO2 C C2O3 C FORM OPEN C O3 / 0.03 ALD2 C 0.62 C2O3 C 0.7 FORM C 0.03 XO2 C 0.69 CO C 0.08 OH C 0.76 HO2 C 0.2 GLY OH C XYL / 0.7 HO2 C 0.5 XO2 C 0.2 CRES C 0.8 MGLY C 1.1 PAR C 0.3 TO2 OH C MGLY / XO2 C C2O3 MGLY / C2O3 C HO2 C CO O C ISOP / 0.6 HO2 C 0.8 ALD2 C 0.55 OLE C 0.5 XO2 C 0.5 CO C 0.45 ETH C 0.90 PAR OH C ISOP / XO2 C FORM C 0.67 HO2 C 0.13 XO2N C ETH C 0.2 ALD2 C 0.4 MGLY C 0.2 C2O3 O3 C ISOP / FORM C 0.4 ALD2 C 0.55 ETH C 0.2 MGLY C 0.1 PAR C 0.06 CO C 0.44 HO2 C 0.1 OH NO3 C ISOP / ! O2N XO2 C NO / NO2 XO2 C XO2 / XO2N C NO / NR / NR APIN / 7 PAR C 0.75 ALD2 C 0.75 OLE UNKN / 7.5 PAR C 0.875 ALD2 C 0.375 OLE NO3 / NO
6.100 ! 104a
PAN / C2O3 C NO2 C2O3 C C2O3 / 2 ! O2 C 2 FORM C 2 HO2 C2O3 C HO2 / 0.79 FORM C 0.79 XO2 C 0.79 HO2 C 0.79 OH OH / XO2 C FORM C HO2 PAR C OH / 0.87 XO2 C 0.13 XO2N C 0.11 HO2 C 0.11 ALD2 C 0.76 ROR – 0.11 PAR ROR / 1.1 ALD2 C 0.96 XO2 C 0.94 HO2 – 2.10 PAR C 0.04 XO2N C 0.02 ROR ROR / HO2 ROR C NO2 / O C OLE / 0.63 ALD2 C 0.38 HO2 C 0.28 ! O2 C 0.2 OH C 0.3 CO C 0.2 FORM C 0.02 XO2N C 0.22 PAR OH C OLE / FORM C ALD2 C XO2 C HO2 – 1 PAR O3 C OLE / 0.5 ALD2 C 0.74 FORM C 0.33 CO C 0.44 HO2 C 0.22 XO2 C 0.1 OH – 1 PAR NO3 C OLE / 0.91 XO2 C 0.09 XO2N C FORM C ALD2 – 1 PAR C NO2 O C ETH / FORM C 0.7 XO2 C CO C 1.7 HO2 C 0.3 OH OH C ETH / XO2 C 1.56 FORM C HO2 C 0.22 ALD2 O3 C ETH / FORM C 0.42 CO C 0.12 HO2 OH C TOL / 0.08 XO2 C 0.36 CRES C 0.44 HO2 C 0.56 TO2 TO2 C NO / 0.9 NO2 C 0.9 HO2 C 0.9 OPEN TO2 / CRES C HO2
1.739 ! 104 exp(986/temp)a 1.037 ! 104 exp(250/temp)a 3.700a Photolysis 9b 7.915 ! 103 exp(250/temp)a 1.180 ! 104 exp(5500/ temp)a 5.616 ! 1018 exp(14000/ temp)b 3.700 ! 103a 9.600 ! 103a
7.774 ! 103 exp(1815/ temp)b 1.203 ! 103a
6.250 ! 1016 exp(8000/ temp)b 9.545 ! 104b 2.200 ! 104a 1.756 ! 104 exp(324/temp)a
7.740 ! 103 exp(504/ temp)a 2.104 ! 101 exp(2105/ temp)a
a b c
3.250 ! 104a 2.000 ! 104 Photolysis 10b 4.400 ! 104a 8.030 ! 102 exp(500/ temp)a
2.453 ! 104 exp(116/temp)a
2.600 ! 104a Photolysis 11b 2.700 ! 104a
1.420 ! 105a
1.800 ! 102a
4.700 ! 102a 1.200 ! 104a 2.550 ! 101 exp(1300/temp)a 1.000 ! 103a 1.000b 1.000 ! 1010b 1.000 ! 1010b Photolysis 12b
[cm3 molec1 s1]. [s1]. [cm6 molec2 s1].
1.135 ! 101a
References
1.540 ! 104 exp(729/temp)a
Bader, G., Deufhlhard, P., 1983. Numerische Mathematik 47, 373. Briggs, G.A., 1969. Plume rise, U.S. Atomic Energy Commission Critical Review Series T/D 25075. Briggs, G.A., 1971. Some recent analyses of plume rise observations. In: Englund, H.M., Beery, W.T. (Eds.), Proceedings of the Second International Clean Air Congress. Academic Press, New York, pp. 1029–1032. Briggs, G.A., 1974. Diffusion estimation for small emissions, Environmental research laboratories air resources atmospheric turbulence and diffusion laboratory 1973 annual report, USAEC Rep ATDL-106 Natl. Oceanic Atmos. Admin., Washington, D.C. Canepa, E., Ratto, C.F., 2003. Algorithms to simulate the transport of pollutant elements: a model validation exercise and sensitivity analysis. Environmental Modelling and Software 18, 365–372.
3.000 ! 103 exp(411/temp)a
1.856 ! 101 exp(2633/ temp)a 3.106 ! 103 exp(322/temp)a
1.200 ! 104a 2.500 ! 102b
1250
L.E. Olcese, B.M. Toselli / Environmental Modelling & Software 20 (2005) 1239–1250
Carruthers, D.J., Holroyd, R.J., Hunt, J.C.R., Weng, W.-S., Robins, A.G., Thomson, D.J., Smith, F.B., 1994. UK-ADMS, a new approach to modelling dispersion in the earth’s atmospheric boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 52, 139–153. Duff, I.S., Reid, J.K., 1984. Subroutine to Solve a General Sparse System of Linear Equations. Rutherford Appleton Lab. Gabruk, R.S., Sykes, R.I., Seigneur, C., Pai, P., Gillespie, P., Bergstrom, R.W., Saxena, P., 1999. Evaluation of the reactive and optics model of emissions (ROME). Atmospheric Environment 33, 383–399. Gear, C.W., 1971. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, New Jersey, USA. Gery, M.W., Whitten, G.Z., Killus, J.P., Dodge, M.C., 1989. A photochemical kinetics mechanism for urban and regional scale computer modeling. Journal of Geophysical Research 94, 12945– 12956. Goyal, P., Sidhartha, 2004. Modeling and monitoring of suspended particulate matter from Badarpur thermal power station, Delhi. Environmental Modelling and Software 19, 383–390. Hanna, S.R., Chang, J.C., 1991. Modification of HPDM for Urban Conditions and Its Evaluation using the Indianapolis Data Set. Final Report Prepared for EPRI by EARTH TECH, 196 Baker Ave., Concord. Hov, Ø., Zlatev, Z., Berkowicz, R., Eliassen, A., Prahm, L.P., 1989. Comparison of numerical techniques for use in air pollution models with non-linear chemical reactions. Atmospheric Environment 23, 967–983. Madronich, S., 1993. UV radiation in the natural and perturbed atmosphere. In: Tevini, M. (Ed.), Environmental Effects of UV Radiation. Lewis Publishers, Boca Raton, pp. 17–69. Morris, E.R., Chang, E.C., Shepard, S.B., Ligocki, M.P., 1992. ‘‘User’s Guide to Version IV of the Reactive Plume Model (RPMIV)’’. (SYSAPP-92/037), Systems Applications International, San Rafael, California. Olcese, L.E., Toselli, B.M., 1998. Fast and reliable numerical methods to simulate complex chemical kinetic mechanisms. International Journal of Chemical Kinetics 30, 349–358. Olesen H.R., 1994. Operational Short-Range Atmospheric Dispersion Models for Environmental Impact Assessments in Europe. Pasquill, F., 1961. The estimation of the dispersion of windborne material. Meteorology Magazine 90, 33–49.
Seinfeld, J.H., 1986. Atmospheric Chemistry and Physics of Air Pollution. John Wiley & Sons, New York. Stewart, D.A., Liu, M.K., 1981. Development and application of a reactive plume model. Atmospheric Environment 15, 2377–2393. Turner, D.B., 1997. The long lifetime of the dispersion methods in Pasquill in U.S. regulatory air modeling. Journal of Applied Meteorology 36, 1016–1020. US EPA, 1995a. User’s Guide for the Industrial Source Complex (ISC3) Dispersion Models Volume I – User Instructions, EPA-454/ b-95-003a. US EPA, 1995b. User’s Guide for the Industrial Source Complex (ISC3) Dispersion Models Volume II – Description of Model Algorithms, EPA-454/b-95-003b. Beatriz M. Toselli Professor Toselli received her Bachelor Degree in Physical Chemistry from The Universidad Nacional of Co´rdoba, Argentina in 1981 and her Ph. D. in Chemistry from the same University in 1986. She has been a Post Doc at The Space Physics Research Laboratory. The University of Michigan, Ann Arbor, Michigan, U.S.A. from 1987–1990 and a Research Scientist in the Laser Photochemistry Group, Atomic Energy of Canada Limited, Chalk River Laboratories, Canada from 1990–1991. Currently she is an Associate Chemistry Professor at the Department of Physical Chemistry, Universidad Nacional de Co´rdoba and an Independent Research Scientist from the National Research Council of Argentina (CONICET). She is the director of the Air Pollution Group of the department and her research interest are on energy transfer, laser spectroscopy chemistry and air pollution and radiative processes in the atmosphere.
Luis E. Olcese Dr. Olcese received his Bachelor Degree in Chemistry from The Universidad Nacional of Co´rdoba, Argentina in 1994 and his Ph. D. in Chemistry from the same University in 2000. Currently he is a Teaching Assistant at the Department of Physical Chemistry, Universidad Nacional de Co´rdoba and an Assistant Research Scientist from the National Research Council of Argentina (CONICET). He is one of the researchers of the Air Pollution Group and his research interests are atmospheric chemistry, mathematical models and air pollution models development.