Abstract

An unsteady state Lagrangian diffusion model with consideration for the removal of sulphur dioxide through chemical transformation in the planetary boundary layer has been developed and evaluated in comparison with 2 existing models. The data used were from episodes of industrial emissions in a sub-urban area monitoring studies in South Africa. The effect of varying transformation rate was investigated using the three different overall rate constants for the advection of SO₂ in smoke plume between a source (power station smoke stack) and a receptor point (ground based sampling site) under the influence of the local meteorology. The concentration distribution in the along-wind direction revealed that SO₂ oxidation rate as second order with deposition as first order predicted the observed ambient concentration to an accuracy of about 91.9%.

---

INTRODUCTION

The concept of incorporating the reaction terms in atmospheric diffusion enables the forecasting of concentration distribution of chemically reactive species in the along-wind direction during air transport (Lusis and Phillips, 1977; Omstedt and Rodhe, 1978; Carmichael and Peters, 1984). The signature of any homogeneous chemical reaction is unique with specific atmospheric conditions, therefore as meteorological conditions varies the mixture concentrations and compositions varies thereby altering the reaction pattern. This has necessitated the investigation of the behaviour of SO₂ oxidation and deposition during diffusion with advection in the planetary boundary layer over Mpumalanga Highveld area of South Africa.

Over southern Africa, regional emissions of SO₂ is significantly high especially over the northern sub-region of the South African Highveld due to the high density of SO₂ emitting industries which include from fuel refining process industries, coal-fired power plants and open cast coal mines (Terblanche et al., 1993; Held et al., 1993, 1994). This sub-region which is dominated by grassland vegetation is characterised by mostly dry weather with less than annual rainfall with conditions ranging from unstable to neutral atmosphere (Held et al., 1996). These properties allows for a well-mixed daytime composition of windborne materials dispersion reaching the ground and throughout the mixing layer.

In this study, an unsteady state Lagrangian diffusion model was developed and evaluated to predict the concentration distribution of tropospheric SO₂ undergoing a homogeneous oxidation with dry deposition as the removal mechanism in a smoke plume during the advection (Fig. 1). The model was formulated to predict the downwind mass concentration of SO₂ after oxidation to sulphate similar to that applied in Lusis and Phillips (1977) with vertical height limited by the inversion layer (Pasquill, 1974). In earlier studies, atmospheric SO₂ oxidation rate applied in diffusion model was often estimated as first orders in the rate expression while in recent studies it is shown as a second order rate constant (Seinfeld and Pandis, 1998; Herrmann et al., 2000; Grgic and Bercic, 2001). The oxidation rate was developed utilising the reactant and product sulphur species sulphate (Igbafe, 2007; Levenspiel, 1999; Ronneau and Snappe-Jacob, 1978; Alkezweeny and Powell, 1977) measured with continuous SO₂ and sulphate analysers with short averaging time. The degree of accuracy of the diffusion model was examined with the comparative study between the output of the second order SO₂ oxidation model (Igbafe, 2007) and two existing first order models, an atmospheric simulated smog chamber study (Carmichael and Peters, 1984) and from published rate constant data (Möller, 1980).

MODEL DEVELOPMENT

The diffusion models were proposed based on the following assumptions:

•	All atmospheric diffusion of gases and submicron particles (such as sulphates) are equal and that the vertical distribution of reactant and product concentrations is of the Gaussian pattern throughout the mixing layer (Fig. 1)

•	That the rate of reaction is second order (Freiberg, 1975; Seinfeld and Pandis, 1998; Warneck, 1999; Levenspiel, 1999; Herrmann et al., 2000; Grgic and Bercic, 2001) and the deposition rate is first order (Alkezweeny and Powell, 1977; Ronneau and Snappe-Jacob, 1978; Sehmel, 1980)

•	Sulphates are formed from the transformation of emitted sulphur dioxide in the puff

•	Background concentrations have negligible effect on the emitted sulphur concentration and insignificant contributions on the total sulphur at source point

•	That concentration changes due to diffusion exceed that due to chemical reactions

•

That wind directions at ground level are considered to fairly good indications of the origin of polluted air masses. Boundary layer wind speeds were estimated from surface winds based on power exponent’s law described by Seinfeld and Pandis (1998), De Nevers (2000) and Tyson and Preston-Whyte (2000)

•	That from defined source points, the travelled times of air advection through the sampling site depends on the surface winds and mixing heights

Chemical sub-model: The oxidation of SO₂ has been described as a very slow process in the formation of sulphate and a result it is the rate-controlling step. This stage is followed by a faster reaction between the sulphur trioxide formed from SO₂ oxidation and water vapour (Calvert et al., 1978; Pienaar and Helas, 1996; Warneck, 1999; Herrmann et al., 2000). The first oxidation rate model applied in the diffusion model was based on the conservation of estimated mass for the total sulphur in air present as the sulphur dioxide and sulphate in plume from the coal-fired power plants downwind traversing an air monitoring station 20 km away (Igbafe, 2007). The kinetic model designated as model type-1 during advection is:

(1)

Where:

 

T (K)	=	The ambient temperature
SO2	=	The sulphur dioxide concentration at any given time t, in the mixing layer

 

The second rate expression designated as model type-2 was derived from smog chamber studies for tropospheric SO₂ oxidation (Carmichael and Peters, 1984) -involved 72 reactions with major species including NO, NO₂, SO₂, O₃, HNO₂, HNO₃, NO₃, H₂O₂, CO, CO₂, SO₃ and HSO₃-, is given by

(2)

The third rate model described by Möller (1980) and designated as model type-3 was derived from published rate constant data of 67 reactions of SO₂ oxidation with different reaction mechanisms to account for the varying atmospheric conditions. The mechanisms considered were photochemical oxidation, homogeneous gas phase oxidation with radicals, liquid phase oxidation and gas to particle conversion.

 

Fig. 1:	Diffusion model pattern for the disappearance SO2 to sulphate

 

(3)

Equation (1-3), are the expressions for the disappearance by reaction of SO₂ in the atmosphere required in the diffusion model.

Atmospheric diffusion during advection: The diffusion approach proposed in this paper is similar to that established in Lusis and Phillips (1977) where the transformation rate constant of the chemical sub-model is put into Eq. (4) expression. In reaction kinetics the rate of a chemical reaction is expressed either as a function of the disappearances of the reactants or the formation of the products. Since, the detailed chemical reactions involved in the formation of particulate sulphate are not completely discussed in this study, it is more convenient to express the transformation only in-terms of the disappearance of SO₂ by oxidation and deposition. Hence, formulated rate expressions of Eq. (1-3) are incorporated into diffusion with advection model. For stationary situations and homogeneous turbulence the common atmospheric diffusion formula for the mean concentration of a species emitted from a continuous, elevated point source is the Gaussian formula (Seinfeld and Pandis, 1998). Based on the mean wind speed in the along-wind direction coupled with the assumption that the mean crosswind speed and mean vertical wind speed equals zero, the equation for the contribution of a puff at a receptor point from a continuous point source for real time atmospheric application of partial absorption at z = 0 and the presence of an impermeable upper boundary with 0≤z≤H (i.e., inversion layer), is given by Eq. (4). In Eq. (4) the mean concentration (μg m^-3) at the receptor point (x, y, z) at time t expressed in terms of the emission rate, q, at the source point (x₀, y₀, z₀) at time t₀ given as:

(4)

Where,

(4a)

 

Where:

 

 

q (μg)	=	The mass of specie in puff
(m s-1)	=	The mean wind speed at emission height
(m2 s-1)	=	The mean eddy diffusivity in the along-wind direction
(m2 s-1)	=	The mean eddy diffusivity in the crosswind direction

 

(m2 s-1)	=	The mean eddy diffusivity in the vertical direction
x (m)	=	The downwind receptor point from puff centre
y (m)	=	The crosswind receptor point from puff centre
z (m)	=	The vertical receptor height from ground

H (m)	=	The inversion layer height above ground
g (m)	=	The vertical term

Given the following boundary conditions that equals zero as distances x and y approaches infinity and that:

Where, v_d is the deposition velocity that measures the degree of absorption of the Earth’s surface while K_zz is the vertical eddy diffusivity. At steady state is assumed

(Seinfeld and Pandis, 1998). Replacing the diffusivities with dispersion coefficients and simplifying Eq. (4) gives:

(5)

 

(5a)

 

Where, σ_x, σ_y and σ_z are the along-wind, crosswind and vertical dispersion coefficients, respectively. When the principle of conservation of mass for an unsteady state emission from a continuous point source is applied to atmospheric diffusion with transformation, the changes in emission source strength q with time is related to the reaction rate, η, according to Eq. (6):

(6)

Equation (6) is similar to that described by Lusis and Phillips (1977). Combining Eq. (1), (5) into (6) and (2), (5) into (6) as well as (3), (5) into (6) with both rate and deposition constants expressed in per seconds gives:

(7)

Where:

 

(7a)

 

(8)

(9)

Integrating Eq. (7-9) for x-, y-, z-coordinates over the travelled time (between the emission source and monitoring station) by the reacting fluid yields the following: Equation (7) reduces to:

(10)

Where,

(10a)

While, Eq. (8) becomes:

(11)

And Eq. (9) simplifies to:

(12)

In Eq. (10-12), the summation term was truncated at n = 15 because beyond this value of n, the exponential terms were approximately zero. To further simplify these equations, the smoke puff was assumed to be horizontally symmetric, where, σ_x is made equal to σ_y. In order to, apply these models to the area under study (Mpumalanga Highveld), with the boundary layer condition of the reference height, z, relative to the Monin-Obukhov length, L, of z/L<0, the deviation of wind velocities (dispersion coefficients) in the crosswind and vertical directions (Seinfeld and Pandis, 1998) are necessary. For crosswind direction, the dispersion coefficient is expressed as:

(13)

Where,

(13a)

While, for the vertical direction, it is expressed as:

(14)

Where,

(14a)

Where, f is the Lagrangian time scale that specifies the characteristics of the atmospheric boundary layer. The terms u^* and w_* are the friction velocity and convective velocity scale, respectively, while, h_ux is the mixed layer height and L is the Monin-Obukhov length. Simplifying Eq. (10-12) in terms of σ_v and σ_w gives:

(15)

Where,

(15a)

And

(15b)

While, rearranging Eq. (11) gives to:

(16)

And Eq. (12) becomes:

(17)

The terms u^*, w_*, h_us and L, dictate the characteristic behaviour of the tropospheric boundary layer are referred to as the boundary layer parameters. The procedure given below was used to determine the boundary layer parameters.

Determination of the boundary layer parameters: In this study, the expressions used for calculating the planetary boundary layer parameters for convective (daytime) conditions are described below. According to Oke (1987), the sensible heat of the overall solar radiation heat flux in the boundary layer is dominant by the sensible heat compared to the latent heat, therefore, the sensible heat flux, φ is determined from the heat balance formula simplified into:

(18)

Where:

 

Rn (W m-2)

=

The net radiation

 

Equation (18) assumes that the soil heat flux is 10% of the net radiation (Holtslag and Van Ulden, 1983) and a 0.8 for the Bowen ratio for grassland as given by Oke (1987). The frictional velocity u_* (m s^-3) and Monin-Obukhov length L(m), which are interrelated parameters were calculated with methods described in literatures (Panofsky and Dutton, 1984; Holtslag and Van Ulden, 1983; Van Ulden and Holtslag, 1985; Perry, 1992) expressed as:

(19)

The determination of the frictional velocity u_* is based on the following conditions defined by the ratio of reference height for wind speed, z, to the Monin-Obukhov length, L. Therefore, when

(20)

and

(21)

 

(21a)

 

(22)

Where, k is von Karman’s constant (k = 0.41). Ψ_m and μ are unit-less quantities. A reference height for wind speed of 10 m was used for the area since it was a grassland and the roughness height z₀ was 0.05 m for grass of between 0.25 and 1.0 m (Oke, 1987; Seinfeld and Pandis, 1998), u_z is the wind speed at reference height, z. u_* was determined by initially assuming a neutral condition with Ψ_m = 0 and obtaining an initial u_* using Eq. (19). With u_* determined, Monin-Obukhov length L^N was calculated from Eq. (23) (Holtslag and Van Ulden, 1983; Wyngaard, 1988) as:

(23)

The calculated L^N was then substituted as L into (19) and a new friction velocity u_*^N+1 was calculated which was used to calculate the new Monin-Obukhov length L^N+1 (Venkatram, 1980b). This iteration continues until |u_*^N+1-u_*^N|≤0.01 and |L^N+1-L^N|≤0.01 has been satisfied. φ is the sensible heat flux, T(K) is the temperature at height z. From Van Ulden and Holtslag (1985), the time dependent convective mixing heights for unstable h_us and neutral h_u conditions (since, it is almost impracticable to achieve stable conditions in the daytime) were calculated from the combination of a number of energy balance model rate equations discussed by Deardorff et al. (1980), Venkatram (1980a), Van Ulden and Holtslag (1985), Garratt (1992) and Cimorelli et al. (2004) utilizing morning potential temperature sonding before sunrise as well as the hourly changing surface heat fluxes. The mixed layer heights are given by:

(24)

(25)

Where, f is the Coriolis parameter estimated as 1.21x10^-4 s^-1 for the Highveld area, t(h) is the time of day and u (m s^-3) is the mean wind velocity. Another boundary layer parameter required for evaluation of the dispersion coefficients is the convective velocity scale w_*. The convective velocity scale required as a result of turbulent eddies in the vertical dispersion profile was determined using the buoyant production of turbulent kinetic energy expression in Van Ulden and Holtslag (1985) with the calculated h_us depending on the L obtained as:

 

Table 1:	Seasonal mixing heights over specific temperature range

 

(26)

In this study, the calculated mean convective boundary layer height for both unstable and neutral conditions on a seasonal basis are shown in Table 1. With the terms u^*, w_*, h_us and L known, the deviation of wind velocities σ_v and σ_w can be determined. When the dispersion coefficients are known, the mass concentrations distribution of the diffusing mixture can be obtained.

Determination of dispersed SO₂ relative to the source strength: According to Tyson and Preston-Whyte (2000), studies over the entire southern Africa, have shown that surface inversion heights over the Highveld ranged between 400 and 600 m. The model in this study was simplified by assuming, a constant value of 500 m as the daytime surface inversion height. Simplifying Eq. (15-17) for the mass fraction of SO₂ distribution between two points within the boundary layer in the daytime generate the following expressions.

Integrating Eq. (15) [model type-1] with respect to vertical distribution and rearranging in terms of the displacement time gives into the form:

(27)

Where,

(27a)

 

(27b)

 

 

(27c)

 

 

(27d)

 

 

(27e)

 

 

(27f)

 

While, for the model type-2, Eq. (16) becomes:

(28)

And for the model type-3, Eq. (17) gives:

(29)

Equation (27-29) represent the contribution of atmospheric stability on emission rates variations. They are developed for daytime emission rate change within the convective boundary layer for unstable condition at z <50 m. Integrating Eq. (27-29) numerically using Simpson’s formula defined in Eq. (30), from t₁ the start time to t₂ the time taken to reach the observation point yielded the mass fraction for SO₂ (Q_s/Q⁰_s)₂ present in a air parcel in puff at location 2 with a known plume age t₃ in relation to the fraction of SO₂ (Q_s/Q⁰_s)₁ remaining at an earlier plume location 1 and age t₁ within the day. If s_q is a downwind distance at any fixed point q from emission source between location 1 and 2, then t_q the travelling time is evaluated by:

(30)

(31)

(32)

t₁ and t₂ are the start and finish times respectively and κ is an even number of divisions between the start and finish times, while q are equal distant locations generated from theκ-divisions. The integral function Ψ_T(t) represents the sum of the individual integral function at a particular time t. The emission rate changes of SO₂ content downwind (Q_s/Q⁰_s)₂ over Elandsfontein area of the Mpumalanga Province were estimated with reference to the source emission rate (Q_s/Q⁰_s)₁ using Eq. (27-29). The advection model was tested with ambient air samples of sulphur dioxide at two distances from the emission source upwind the sampling site for 12 months.

EXPERIMENTAL

The oxidation of SO₂ has been described as a very slow process (Calvert et al., 1978; Seinfeld and Pandis, 1998) and is the rate-controlling step in the formation of sulphate (Pienaar and Helas, 1996) while, the reaction between the trioxide and water vapour is a faster reaction (Pienaar and Helas, 1996). It could then be assumed that the transformation into trioxide and sulphate are negligible at emission source. Hence, in addition to the model development field measurements were conducted for the deduction of SO₂ disappearances with advection.

The power station emission stacks are fairly close to each other about 200 m apart. The combined capacity of the two coal-fired power plants is 5400 MW with 10 emission stacks with 6 grouped three per stack-casing an the other 4 grouped 2 per stack-casing. Each stack has approximately 3 m inner-diameter by 250 m physical height.

Sampling methodology: Ambient air sampling episodes for evaluating the formulated diffusion model was conducted at Elandsfontein area of the Mpumalanga Highveld between September 2004 and August 2005. Elandsfontein is about 20 km away from two adjacent coal-fired power stations. These coal-fired power stations are the only significantly large sources of air pollution around the area. At Elandsfontein is located an air quality monitoring station free of adjacent sources of SO₂, this was to attribute the sole responsibility for all the episodes of SO₂ to the industries located near the sub urban town. It was equipped with instruments for ambient concentrations and meteorological measurements. The sampling site is an open area 1600 m above sea level with fairly flat topography and few undulating hills and sparse vegetation. Air masses and air trajectories allow for accurate indications of the origin of polluted (Ronneau and Snappe-Jacob, 1978). With defined source points, the travel time of transported pollutant from the source downwind of the sampling site is a function of the meteorology and could be evaluated based on the changing mixing heights with time of day. Based on the mean wind velocities of advecting air mass to the sampling site, the travelling time was estimated. It was assumed that for transformation purpose all SO₂ and sulphates emerged from the industrial emission source area since it was the only significant air pollution contributor within the area as well as data from the wind directional sector towards the power stations were most significant relative to other directions.

Sampling procedure: Ambient air sampling was conducted for 12 months on a continuous basis, for ambient concentration of SO₂ and particulate sulphate. SO₂ concentrations were measured with a Thermo Electron Corporation Instrument Model 43 C SO₂ pulse fluorescence analyser while the particulate sulphate was monitored with the Rupprecht and Patashnick automated Series 8400 S ambient particulate sulphate monitor. The instruments data were stored on a ten-minutes averaging time. Background concentrations of sulphate and SO₂ were determined during episodes of minimum SO₂ concentration (between 0.5 and 1 ppb). The Model 43 C SO₂ trace level analyser measures SO₂ pulse fluorescence released from excited sampled SO₂ gas. A detailed description of the SO₂ pulse fluorescence principle is given in Dittenhoefer and De Pena (1978). The series 8400 S particulate sulphate monitor consists of a pulse generator and a pulse analyser. The pulse generator operates by sampling ambient air with particulates of less than or equal to 2.5 microns through a humidifier to a reaction chamber where it is flashed at about 600°C by a platinum flash strip within 0.01 sec to release a pulse of SO₂ which is transmitted to the pulse analyser (Stolzenburg and Hering, 2000). At the pulse analyser the SO₂ concentration is determined in parts per billion and the signal returned to the generator where it is converted to μg m^-3 for particulate sulphate. The response times for the analysers are one minute for the SO₂ analyser and ten minutes for the particulate sulphate monitor.

RESULTS AND DISCUSSION

Using the concentration closest to the source point as a starting point, at the concentrations farther-out locations can then be predicted. The emission rates of SO₂ downwind (Q_s/Q⁰_s)₂ over Elandsfontein area of the Mpumalanga Province were estimated using Eq. (27-29). These variations in mass fraction are shown in Table 2.

 

Table 2:	Comparison of predicted and observed SO2

 

The outcome revealed that model type-1 with an R² of 0.9198 predicted the measured more accurately than model type-2 with an R² of 0.606 and much more compared to model type-3 with an R² of 0.596. Based on the coefficient of determination the trend-line analysis of the three model predictions are given by Eq. (33-35). That is, with model type-1, the SO₂ mass fraction predicted as against the measured is:

(33)

Whereas, with model type-2, the SO₂ mass fractional distribution is expressed as:

(34)

And with model type-3, SO₂ may be predicted using:

(35)

Where, y_p and y_m represents the predicted and measured SO₂ mass fractions respectively that were present during atmospheric diffusion of sulphur-containing smoke puff under the influence of transformation.

Analyses of the concentration distribution of SO₂ during atmospheric diffusion with chemical transformation and removal of sulphur dioxide in smoke puffs in the PBL have been established. Chemical sub-models of both first and second-order reaction types were applied in incorporated within the time-dependent Lagrangian puff model in order to predict the extent of SO₂ disappearances resulting from transformation during advection. Atmospheric diffusion of SO2 under the influence of chemical transformation was better predicted when the reaction kinetics was of a second order than when considered as a first order reaction. Consequently for subsequent prediction under same tropospheric conditions of Table 3, the expression of Eq. (33) is most suitable for the determination of ambient SO₂ at any location when the source strength is known.

 

Table 3:	Meteorological conditions throughout sampling period

 

ACKNOWLEDGEMENT

The authors wish to express their appreciation to Mr. R. Burger Climatology Research Group, University of the Witwatersrand for his suggestions. This work was sponsored by Climatology Research Group, University of the Witwatersrand, Johannesburg, South Africa.

Appendix

How to cite this article:

A.I. Igbafe and S.J. Piketh . Variations in the Order of Reaction: An Implication on Atmospheric Diffusion with Transformation and Deposition of SO₂ in Smoke Plume.
DOI: https://doi.org/10.36478/erj.2009.4.13
URL: https://www.makhillpublications.co/view-article/1994-5396/erj.2009.4.13