Environ Eng Res > Volume 27(6); 2022 > Article
Eslami, Fatehifar, and Kaynejad: A 3D mathematical evaluation of the emission into the air of reactive BTEX compounds: A new approach for mechanism reduction

### Abstract

BTEX species are abundant volatile organic compounds that are classified as main pollutants by several environmental protection agencies. In this study, a new 3D air pollution dispersion model with the capability of taking into account the BTEX’s reactions was developed and evaluated. The Toxchem model is employed to estimate the amount of emitted BTEX from several area sources, followed by simulating the transport of the reactive species through a 3D model. Also, reduced mechanisms were developed, based on Master Chemical Mechanism (MCM), for the simulation of the atmospheric chemistry of BTEX. The application of the mechanism reduction method yielded a mechanism of 43 species and 45 reactions. Based on findings, the deviation of the reduced mechanism from the whole mechanism was < 4% throughout the simulation. In addition, the result showed that, during a 12 h period of simulation, the effect of the atmospheric chemical reaction on reducing the final concentration of benzene, toluene, ethylbenzene, m-xylene, o-xylene, and p-xylene was about 0.3, 1, 1.4, 5, 3.1, and 3.2%, respectively. Lastly, it was indicated that the estimated emission rates by Toxchem and simulated concentrations by the dispersion model were in good agreement with the reported experimental values.

### 1. Introduction

BTEX group, including benzene, toluene, ethylbenzene, and xylene, is considered as the most typical component of volatile organic compounds (VOCs) pollution in the air. Currently, it is well-documented that the toxicity, mutagenicity, and carcinogenicity of BTEX compounds lead to a major public health problem even at low concentrations [1]. Short-term allergic diseases such as headache, dyspnea, fatigue, nausea, itchy skin, mouth ulcers, sore throats, mental confusion, and dizziness, along with long-term exposure problems including leukemia, central nervous system problems, dysfunction of speech, and kidney and liver tissue, are considered as just a few problems which threaten human health by BTEX [14]. In addition to their short- and long-term health troubles, BTEX compounds have created some growing global concerns regarding their detrimental influences on the environment. They adversely affect aquatic, terrestrial, and aerial ecosystems, intensify the global greenhouse effect and form tropospheric ozone and other secondary pollutants [1, 4, 5].
To evaluate the fate of a pollutant in the environment, the first step is to estimate its emission rate from each polluting unit. A review of the different methods available to estimate VOC emissions indicated that Toxchem [68], BASTE [9], and WATER9 [1012] have been widely employed in previous researches. The development of these emission models is mainly based on the mass balance of several compounds for each unit operation facility in which volatilization, sorption, and biodegradation are considered as the main removal mechanisms [6]. For simplicity, in these general fate models, the pollutant concentration inside the units is assumed to be homogeneous and liquid levels are considered constant [13]. Besides, the models need some parameters such as influent pollutant concentration and operational parameters of units to predict the emission rate [6].
Several researchers focused their attention on estimating the emission rates of VOCs from different sources, especially from wastewater treatment plants, through mathematical models. For example, Benis et al. [10] evaluated the fate of acrylonitrile and styrene around an industrial wastewater treatment plant using the WATER9 model. Behnami et al. [6] assessed the fate of benzene, toluene, and styrene in a full-scale petrochemical wastewater treatment plant using the Toxchem model. They revealed that in the treatment process, the equalization basin is the major source of emission due to its higher concentration of VOCs, larger surface area, and turbulence effects. Shakerkhatibi et al. [7] investigated the fate of styrene and acrylonitrile emitting from the wastewater treatment unit of a petrochemical plant using the Toxchem model. They indicated that the volatilization of styrene and acrylonitrile is the main mechanism for their removal from the wastewater treatment unit. Zwain et al. used Toxchem to simulate the fate and removal mechanisms of phenol in the full-scale moving bed biofilm reactor sewage treatment plant. They showed that Toxchem is an efficient model to predict the fate of phenol [8].
The second step to reach the goal of assessing the fate of pollutants in the air is the estimation of dispersion patterns. Several advanced air quality dispersion models including AERMOD [10, 14, 15], ISCST3 [16], and CALPUFF [17, 18] were used to estimate dispersion patterns. The models not only provide an estimation of ambient chemical concentration in the area of interest but also help to evaluate health risk assessment analysis. In this way, the results of the estimation of pollutant emission rate from the units (first step) along with meteoroidal data, characteristics of sources, and wastewater properties are used as input data for dispersion modeling [6, 10].
A considerable number of studies were carried out in this case. For instance, Behnami et al. [6] employed AERMOD to predict the ground-level concentrations of benzene, toluene, and styrene releasing from 5 emission sources existing in a wastewater treatment plant. Benis et al. [10] used the AERMOD model to simulate the atmospheric behaviors of emitted acrylonitrile and styrene from the wastewater pretreatment plant of Tabriz Petrochemical Company. In another study, Abdul-Wahab used the original ISC (Industrial Source Complex) model to estimate the concentration of pollutants emitted from the Mina Al-Fahal refinery in the Sultanate of Oman [19]. Despite the popularity of mentioned models, some numerical methods are also applied for simulating pollutants dispersion in the air. For example, Fatehifar et al. [20] used a multiple cell model to evaluate the releasing of pollutants from Tabriz oil refinery stacks by considering a variable location of stacks and wind blowing angle. Xiang et al. [21] proposed statistical models combined with a linear regression model and a non-linear logistic one to model the PM2.5 concentration. Assuming fixed atmospheric parameters, Regland et al. [22] introduced a two-dimensional approach for modeling pollutant emissions from surface emission sources. Nourbakhsh et al. [23] simulated the dispersion patterns of gaseous pollutants (CO, SO2, and NO2) emitting from stacks and flares of gas refineries, and observed an acceptable agreement between their CFD model and experimental data. In another study, Issakhov and Baitureyeva developed a 2D and 3D CFD model to numerically investigate the rising pollution from thermal power plants. In addition, they calculated the pollution concentration level at various distances from the emission source [24].
Furthermore, some scholars have focused their attention on the reactive dispersion of VOCs and NOx as the primary pollutants that lead to the production of secondary pollutants. In this way, Baker et al. [25] evaluated the transmission of reactive pollutants on a street. Considering only three reactions for chemical transmissions, they simulated the emission of NOx and O3 by using a large-eddy simulation model and concluded that chemical reaction plays an important role in transmitting pollutants into the street. Similarly, Moradpour et al. [26] assessed the effects of green roofs on reactive pollutant dispersion within an urban street canyon by using a CFD model coupled with NO-NO2-O3 photochemistry and energy balance models. Additionally, Kitabayashi et al. [27] designed a model based on the Gaussian model to estimate NOx concentrations by considering the atmospheric chemical reactions. In another study, Tetzlaff et al. [28] modeled the emission of reactive pollutants into the atmosphere from an industrial zone in southern Germany using a large-scale networking model. Baik et al. [29] used the Navier-Stokes equations to model the emission of NOx and O3 pollutants in urban streets. They found that the chemical transmission of pollutants for O3 was comparable to that of the dispersion and diffusion phenomena, while the magnitude of the chemical reaction term for NOx is negligible to that of the advection or turbulent diffusion term. Besides, several reduced chemical mechanisms have been presented in the literature to describe the atmospheric chemistry of pollutants. For example, Joelsson et al. [30] presented a reduced mechanism, based on the MCM v3.3.1, for modeling atmospheric chemistry. Their reduction method was a semi-stochastic method based on the heuristic Ant Colony Optimization concept. They showed that their reduced mechanism successfully predicts the concentration of ozone, nitrogen oxides, and other important compounds in simulations of several cases up to five simulated days. As another example, Bright et al. [31] developed a reduced chemical scheme comprising 51 chemical species and 136 reactions, based upon a subset of the MCM v3.1 for the simulation of street canyon atmospheric chemical processing. At the first step, the reduction was achieved by removing night time only chemistry. Then, further reduction was attained by eliminating parent compounds and any unique daughter products, which had little effect on the key chemical intermediates under street canyon conditions.
The main objectives of the present study were (1) estimating the amount of emitted BTEX from several area sources existing in a wastewater treatment plant (WWTP), (2) developing a large scale, high resolution, and complementary mathematical model for BTEX dispersion into the air, and (3) introducing a precise mechanism reduction approach for eliminating non-effective reactions and compounds. To achieve the goals, the emission rates of BTEX compounds from several area sources of WWTP in a petrochemical company were estimated by Toxchem (version 4.1) model. Then, the concentration of BTEX compounds was calculated in a wide area of the case study by using a general 3D mathematical model including molecular diffusion, advection, and atmospheric reactions, as well as a multi-box model. A splitting method as a powerful numerical technique was also used along with an optimized and reduced form of the Master Chemical Mechanism (MCM).

### 2.1. Case Study

A petrochemical company located in Iran was considered as the case study. The company mainly manufactures special chemical and polymer materials such as acrylonitrile, butadiene, styrene, polyethylene, and polystyrene. The WWTP of the petrochemical company is designed to treat the considerable amount of produced wastewater with a mean flow rate of 200 m3/h on non-rainy days and a maximum of 230 m3/h on rainy days. The release of extremely volatile pollutants is inevitable since wastewater treatment is usually carried out in outdoor basins. The volatile compounds enter the air through evaporation and air striping processes. Fig. S1 illustrates the schematic flow diagram of the wastewater treatment plant and existing area sources. Moreover, Table S1 displays the characteristics of the entering wastewater such as BTEX concentration, flow rate (Qin), total suspended solids (TSS), total chemical oxygen demand (TCOD), and biochemical oxygen demand (BOD5).

### 2.2. Method for Concentration Measurement

Analysis of BTEX compounds in the wastewater was conducted using an Agilent technologies gas–chromatography system (model: 7890A GC system). Analysis was carried out using a capillary column of 30 m length and 0.32 mm ID gas chromatograph (cp-sil-5) with a flame ionization detector (FID). The temperature of the column was 40°C for the first 1 min and then reached 100°C with the ramp of 20°C/min. The injector and detector temperature were kept fixed at 230°C and 300°C, respectively [32].
The concentration of BTEX in air samples was measured using GC (Agilent-6890) based on the standard method of EPA-TO-17 with an FID detector. The (Agilent-6890) GC can provide fast analysis of gas samples. In this method, the separation and identification of compounds are based on boiling point and polarity [33].

### 2.3. Estimating the Emission Rate of BTEX Compounds to the Atmosphere

Generally, models are considered as convenient tools for determining the volatile organic compound emission from WWTP. Furthermore, they reasonably estimate emission rates without imposing costs and difficulties of direct measurement. Principally, these models are based on the mass balance for each basin in a WWTP. Assuming that the volume of the basin is constant, the mass balance equation is written as Eq. (1) [10, 34]:
##### (1)
$Vl (dCldt)=QCl,0-QCl+RV+RS+Rad+Rab+Rb$
where Cl,0 and Cl are the input and output concentrations of the soluble pollutants inside the basin expressed in mg/m3, respectively, Vl indicates the basin volume (m3), and Q shows the volumetric flow rate of wastewater into the basin (m3/h). RV, RS, Rad, Rab, and Rb are the removal rate of these compounds by evaporation, air stripping, adsorption, liquid absorption, and biological removal (mg/h), respectively. Toxchem is a widely used capable model for predicting contaminants fate and hazardous air pollutants emission within/from wastewater treatment plants [7, 35]. In this study, Toxchem, as a computer-based software including both steady-state and dynamic models, is used to estimate the amount of BTEX emission from a petrochemical wastewater treatment plant [36].
Of all the mechanisms of removal, volatilization, sorption onto solids, and biodegradation are of paramount importance. For a system open to the atmosphere, the rate of volatilization from its surface is given by Eq. (2):
##### (2)
$rv=KvClfnonV$
where rv is rate of volatilization (mg/h), Kv is volatilization mass transfer coefficient (m/h), Cl is concentration of volatile compound in the water (mg/m3), fnon is pH dependent fraction of non-dissociated compound, and V is volume of process vessel (m3).
At equilibrium, in a liquid/solid system, the concentration of a contaminant on a solid sorbent can be expressed by the following linear equation (Eq. (3)):
##### (3)
$Cx=KpCl$
where Cx is concentration of contaminant in solid phase (μg/g), Kp is sorption partition coefficient (L/g), and Cl is concentration of contaminant in liquid phase (μg/L).
In Toxchem model, pollutant biodegradation is assumed to be represented by the Monod expression, as presented in Eq. (4).
##### (4)
$rb=μmY(ClKs+Cl)XV$
where rb is biodegradation rate (mg/h), μm is maximum specific growth rate (1/h), Y is maximum yield coefficient (g/g), X is volatile fraction of total suspended solids (mg/L), Ks is half saturation constant (mg/L), Cl is contaminant liquid phase concentration (mg/L), and V is vessel volume (L).
To estimate emission rates, using Toxchem software, some information such as the dimensions of each treating unit, meteorological data, the concentration of all compounds in the wastewater entering each treatment unit, and equipment specifications including liquid depth, waterfall height, surface area, and weir length, were used as the input data [3739]. Table S2 indicates the dimensions of all units. Toxchem forecast and simulate the distribution of BTEX between air, liquid, and solid phases by solving all the mass balance equations for all of the unit processes in the treatment system, simultaneously [36]. The emission of BTEX compounds into the air from each basin was as input data for dispersion modeling.

### 2.4. Mathematical Dispersion Model and Boundary Conditions

In the present study, a Multi-Box model is used to predict the concentration of BTEX in the atmosphere around the studied petrochemical company. The Multi-Box model used in this study is based on the model of Fatehifar et al. [40], which is developed in terms of applying reaction terms in the main equation of pollutant transport, considering the multiple area sources for pollutant emissions, and using a high-resolution space domain. To achieve the goal, first of all, the domain of the study is spatially divided into some sub-boxes. The selected space domain has a length of 10 km, width of 2 km, and height of 800 m, which is discretized into 2,000 × 400 × 160 boxes (or grid points) with a volume of 5 × 5 × 5 m3. Fig. S2 displays a schematic view of the gridded space domain, the origin of the coordinate system, wind direction, and WWTP location in the petrochemical company.
In the next step, the mass balance equation was applied to each box. Given the basic equation of mass conservation, four major mechanisms, including (a) mass transfer mainly through wind velocity, (b) molecular diffusion mechanism, (c) removal by deposition, and (d) production and consumption by chemical reactions, are at work in a pollution dispersion model. The mathematical description of the mentioned physical and chemical processes in air pollution results in creating a system of partial differential equations, according to Eq. (5) [41, 42]:
##### (5)
$∂Cgj∂t=-∂(uCgj)∂x-∂(vCgj)∂y-∂(wCgj)∂z+∂∂x(Kx∂Cgj∂x)+∂∂y(Ky∂Cgj∂y)+∂∂z(Kz∂Cgj∂z)+Ej-(k1,j+k2,j)Cgj+Pj(Cg1,Cg2,....,Cgq), j=1,2,.....,q$
where the first three terms on the right of this equation represent the advection of component j in three main directions of x, y, and z. The next three terms describe the diffusion in the three corresponding directions. Ej shows the amount of emission from the pollutant sources, and the next term (k1,j+k2,j) indicates the mathematical equations for dry and wet deposition in the studied area. Finally, the last term (Pj) addresses the pollutants’ reaction. u, v, and w are wind velocity components, and defines the concentrations of the chemical species in the atmosphere. The diffusion coefficients are illustrated by Kx, Ky, Kz and the number of pollutants is shown by the parameter q.
Furthermore, the initial and boundary conditions for all components participating in the reaction, except the primary components of the air such as N2, O2, OH, O3, are as follows:
##### (6)
$Initial condition (t=0) Cgj=0 j≠N2,O2,OH,O3,…$
##### (7)
$B.C≠1 for x=0 Cgj=0 j≠N2,O2,OH,O3,…$
##### (8)
$B.C≠3 for y=0 Cgj=0 j≠N2,O2,OH,O3,…$
##### (9)
$B.C≠4 for y=W Cgj=0 j≠N2,O2,OH,O3,…$
##### (10)
$B.C≠5 for z=0 ∂Cgj∂z=0 j≠N2,O2,OH,O3,…$
##### (11)
$B.C≠6 for z=H ∂Cgj∂z=0 j≠N2,O2,OH,O3,…$
In the above equations, W and H are the width and height of the studied region, respectively. The boundary conditions are obtained based on the following concepts and assumptions:
1. The wind entering the area is pollution-free (Eq. (7)).

2. Pollution concentration at y = 0 and y = W remains zero; i.e. they are considered as the boundary conditions at infinity (Eqs. (8) and (9)).

3. The pollutants are reflected from the earth (Eq. (10)).

4. There is no temperature inversion, and the pollutants can move with no limitation in the vertical direction up to the mixing height (Eq. (11)).

5. Temperature is assumed to be constant (isothermal condition).

### 2.5. Numerical Solution

Choosing an appropriate numerical solution depends on the relevant physical and chemical processes governing the problem under study. Despite the great interest of many researchers to find an optimal numerical solution for the PDE of Eq. (5), little success has been achieved in this regard [43]. Therefore, the equation must be split into small components based on the pollutant transfer process, or sub-models, for solving easily. Splitting allows working with a few simple ODEs instead of one hard PDE. In addition, the method makes it possible to apply different numerical methods for solving each ODE equation. Thus, using an optimal method for the numerical solution of each of these equations significantly increases the accuracy of the dispersion model. The interesting features of the splitting method result in low computation time and a better understanding of the importance of each process [41, 44]. Five sub-models achieved by the splitting procedure are based on Eqs. (12)(16) as follows:
##### (12)
$∂Cgj,1∂t=-∂(uCgj,1)∂x-∂(vCgj,1)∂y$
##### (13)
$∂Cgj,2∂t=+∂∂x(Kx∂Cgj,2∂x)+∂∂y(Ky∂Cgj,2∂y)$
##### (14)
$∂Cgj,3∂t=-(k1,j+k2,j)Cgj,3(x,y,z0,t)$
##### (15)
$∂Cgj,4∂t=+Ej(x,y,z0,t)+Pj(Cg1,4,Cg2,4,....,Cgq,4), j=1,2,.....,q$
##### (16)
$∂Cgj,5∂t=-∂(wCgj,5)∂z+∂∂z(Kz∂Cgj,2∂z)$
The above-mentioned five sub-models represent horizontal advection, horizontal diffusion, dry and wet deposition, emission, chemical reaction terms, and vertical exchange, respectively. After calculating the five ODEs separately, Eqs. (12)(16) should be combined in each time-step to obtain the concentration of each component and calculate separately again in the next time interval. In this case, the data obtained at each point are considered as the input data for the next point. It is worth noting that it is difficult to directly apply Eqs. (12)(16) to actual processes in the atmosphere, despite the advantages mentioned above. Thus, the following assumptions are considered:
1. Each of the area sources is considered as several point sources.

2. The wind velocity is assumed to be only in the x-direction and as a function of height.

3. The system has no deposition.

### 2.6. Reaction Description and Method of Mechanism Reduction

Air pollutants are often modeled in passive or inert conditions. In other words, it is assumed that no reaction takes place during the chemical transmission. However, they are often chemically reactive and play a key role in some cases [25, 26]. The kinetic model of a reaction system is usually based on ordinary differential equations, along with a specified concentration at t = 0 called the initial condition (Eq. (17)).
##### (17)
$dCdt=f (k,C),C0=C (t=0)$
where C indicates the concentration of the components, k shows rate constants, and C0 is the initial concentration of components.
A set of possible reactions should be considered for describing atmospheric chemical reactions using differential equations. However, the selection of this reaction group should be optimized, and a wide range of the studied pollutant should be considered. Due to the fact that the chemical transmission of VOCs in the atmosphere can be explained well by oxidation processes, several reaction mechanisms are proposed for the photochemical oxidation of BTEX compounds in the air. The Master Chemical Mechanism (MCM), as a near-explicit reaction set, is a powerful database for inorganic and photolysis reactions that elaborates the detailed degradation of more than 100 emitted VOCs, along with the production of ozone, and other secondary pollutants [5, 45].
A subset of the MCM, including all possible reactions for the atmospheric oxidation of BTEX compounds, is employed in the present study. Here, the reaction rate is strongly dependent on the presence of sunlight. In this way, the interaction between sunlight and light-receiving molecules leads to the production of some free radicals from certain molecules. OH is considered as one of the radicals formed, which plays a key role in understanding the tropospheric chemistry of volatile pollutants. Hydroxyl radicals are produced in the troposphere through the three sets of reactions, including (a) the photolysis of O3, (b) the photolysis of nitrous acid HNO2, and (c) indirectly from the photolysis of formaldehyde HCHO [1, 46]. These reactions are illustrated as follows:
##### (18)
$O3+hν→O2+O1D$
##### (19)
$O1D+H2O→2OH$
##### (20)
$HONO+hν→OH+NO$
##### (21)
$HCHO+hν→H+HCO$
##### (22)
$H+O2→HO2$
##### (23)
$HCO+O2→CO+HO2$
##### (24)
$HO2+NO→OH+NO2$
It is well documented that the hydroxyl radical is a well-known oxidizing agent for VOCs, however there is a competition between VOCs and NOx to react with the OH radical in the troposphere. The hydroxyl radical reacts with the VOC when the VOC concentration is high. On the other hand, the NOx wins the competition when the NOx ratio is higher. Finally, OH reacts with both proportionally when the ratio of NOx/VOC is a certain value [47].
The existence of various generating and consuming reactions for OH radicals causes a comprehensive set of important atmospheric reactions for describing BTEX fate. The complete chemical mechanism extracted from MCM contains 119 species and 339 reactions. However, some of the extracted chemical reactions, and accordingly existing components, have a minor effect on the final concentrations of BTEX, which imposes various computational difficulties, along with their time-consuming process. For example, 119 × 108 ODEs must be handled at every time step for a space domain including 1,000 × 1,000 × 100 grid-points. Therefore, employing a helpful method for eliminating the ineffective reactions is necessary. Reducing the number of reactions not only decreases the computational time but also helps to easily describe the degradation of BTEX compounds, and identify effective components in the final BTEX concentrations.
To describe the phenomena by fewer reactions, the reduction process is performed by evaluating the effect of each reaction on the final concentration of BTEX. In this way, a one-dimensional model over a 12 h period and for fewer grid points is used to simulate the reactive dispersion of pollutants instead of a three-dimensional one. The steps used for reducing the procedure and extracting the effective reactions are as follows:
1. The initial concentration of components is assigned.

2. One-dimensional reactive dispersion simulation is conducted considering a system of equations involving 339 ODEs for a given time.

3. The final concentration of the BTEX is calculated (j = benzene, toluene, ethylbenzene, and xylene).

4. The first reaction (i = 1) is selected.

5. The ith reaction is eliminated from ODEs, and a new 1D simulation is conducted Then, the final concentration of the BTEX is recalculated.

6. The value of the objective function (O.F), which represents the percentage of change in the concentration of the BTEX before and after eliminating ith reaction, is calculated using Eq. (25).

##### (25)
$O.F=1n×∑j=1n|Cgj2-Cgj1|Cgj1×100$
where n represents the number of BTEX compounds.
7. If O.F > 1%, the ith reaction is selected as an effective reaction. Otherwise, it is discarded. It should be noted that the minimum of O.F can be altered optionally depending on the accuracy of the calculation.

8. The Next reaction is selected (i = i+1), and steps (5) to (7) are repeated until checking

the effect of all of the reactions. Fig. 1 displays the above-mentioned procedure.
At the final step, the three-dimensional simulation is performed using fewer reactions and lower species extracted from the algorithm.

### 3.1. Fate of BTEX

Fig. 2 illustrates the fate of each BTEX compound entering the wastewater treatment plant, resulted from the Toxchem model. As shown, a great part of the entering BTEX compounds is emitted into the atmosphere due to the high volatility and operational problems, and a little part is processed with oil separation or sludge treatment unit. The high emission percentage of BTEX into the air, especially for benzene, indicates that the volatilization is the main removal mechanism for VOCs, which is in line with the result of similar researches [7, 10]. It also shows that the WWTP has not sufficient capability to treat BTEX compounds. However, ethylbenzene is treated a little better than other compounds. The difference between the volatilization of different BTEX compounds may be interpreted using their intrinsic characteristics. The vapor pressure and boiling point temperature of BTEX compounds are shown in Table S3. It can be obviously concluded that the higher emission rate of benzene and toluene is related to their low boiling temperature and high vapor pressure. In other words, the lower the boiling temperature and the higher the vapor pressure, the easier it is for the species to evaporate. However, this argument about Ethylbenzene is not true. Although the boiling point of Ethylbenzene is lower than that of the xylene group, and its vapor pressure is slightly higher, it has less evaporation. In this case, it seems that the effect of other factors [6] such as turbulency and wastewater temperature should also be considered.
As shown in Fig. 3, Toxchem indicates that a large portion of the volatile organic compounds is emitted from the Equalization basin, which is in excellent accordance with the result of previous literature [6, 10]. In addition to the Equalization basin, the API unit plays an important role in air pollution in the wastewater treatment plant. It seems that being fully exposed to airflow, the presence of mixing and turbulence in the basin, high temperature of wastewater, high concentration of BTEX, and long residence time are probably the major causes of the great entry rate of BTEX compounds into the atmosphere through the Equalization basin and API. Similarly, the high concentration of BTEX compounds, along with the turbulence of the flow, causes a high level of emission in the trash screen. However, BTEX emission is negligible for the DAF, Clarifier, and Biological treatment units. It indicates that a lot of BTEX gets into the air before reaching these units, and a minor concentration of BTEX compounds is treated in these units. As a result, reducing the temperature and turbulence of wastewater and increasing the flow rate up to the designed value before entering the units can reduce the emission of volatile organic compounds [6].
Table S4 indicates the mass flow rate of the emission into the air of BTEX compounds, along with their corresponding input values. As shown, out of 1,180 kg of total entering BTEX, the huge amount of 868 kg (73.6%) of BTEX pollutants is released from WWTP into the atmosphere every day, among which benzene is a major contributor to the pollution with about 672 kg. The large volume of the emitted pollutants, in the long term, undermines the health of the staff and the surrounding environment. Therefore, the necessary measures should be taken accordingly. Since the process or equipment modification of petrochemical plants, in most cases, is not feasible [6], it seems that applying air stripping, steam stripping, and covering technologies, as the most applicable methods are the best choice for controlling VOCs emissions from the wastewater treatment unit. As most of the BTEX compounds are emitted from the API and Equalization basin (according to Fig. 3), covering them is a priority.
Furthermore, the estimated BTEX concentrations are compared to the measured values for API, Equalization, and Clarifier (plant outlet) treatment units to evaluate the capability and accuracy of Toxchem in calculating pollutants concentration (Table 1). The 4–8% difference demonstrates that the Toxchem has a high potential for estimating BTEX compounds [20]. Surprisingly, the model accuracy is greater when the concentration of a specified pollutant is higher, and vice versa. Thus, it can be concluded that the Toxchem has an excellent prediction capacity for higher concentrations compared to lower ones.

### 3.2. 1-D Reactive Transmission Analysis

Fig. 4 displays the effect of one-dimensional reactive transmission time on the conversion of BTEX compounds. Here, the conversion is defined as the ratio of the reacted material (at all locations) to the total emitted amounts at a specified time. As shown, the conversion of BTEX compounds increases by increasing the reaction time. It can be stated that the rise of conversion is, chiefly, for a couple of reasons. Firstly, as time goes on, the BTEX compounds react with more hydroxyl radical at each location. in the second place, the reactants migrate to new locations where OH radicals have not yet been used [47]. It is also observed that among all of the compounds, M-xylene is more affected by reaction and can be converted up to 3.1% over a 12 h period. Additionally, the negligible variation of benzene means that it can be classified among non-reactive compounds in the atmosphere Differences among the conversion rate of different VOCs can be interpreted using their corresponding rate constants (shown in Table S5). As can be seen, a reaction with a higher rate constant has a higher conversion rate [47].
In addition, non-effective reactions are recognized and eliminated from calculations by applying the mechanism reduction methodology, over a 12 h period. Accordingly, with a relative difference of 3–4%, a mechanism of 43 species and 45 reactions (13% of all reactions) has almost the same performance of the whole mechanism containing 339 reactions of 119 chemical species for the concentration predictions of BTEX. The obtained result is comparable to other previous studies in terms of compactness of reduced mechanism and accuracy of concentration prediction [30, 31]. Table S5 demonstrates all reactions existing in the reduced mechanism, along with their rate constants. It should be noted that some of the species in Table S5 are not actual species but are reactive intermediate species that allow VOCs to be broken down into smaller fragments. In addition, the index assigned to each species is the total number of C–H and C–C bonds it contains. For example, UDCARB14 is a representative carbonyl intermediate species with 14 C-C and C-H bonds. RA19AO2 is also a representative peroxy radical species with 19 C-C and C-H bonds [48, 49].
As shown in Table S5, the first twelve rows indicate the major reactions for BTEX consumption in which OH radical reacts with BTEX leading to the production of degraded species. Simultaneously, OH radicals are regenerated by the reaction of O1D, HONO, HO2, and NO through reactions 17, 21, and 22. However, as a serious competitor for BTEX, the existing reactants in reactions 14, 25, and 26 react with OH and limit BTEX consumption. Furthermore, some reactions such as 13, 18, 23, and 24 are considered as the competitors for OH because they consume O1D and NO components, which are regarded as some of the main sources of OH production. The adverse or favorable effect of other reactions on hydroxyl production can be similarly explained. The selected reactions and corresponding species are employed to participate in the 3D mathematical model and calculate the spatial concentration distribution.

### 3.3. 3-D model results

Fig. 5 illustrates the ground-level simulation result of the dispersion pattern of reactive BTEX compounds, in the cubic domain under study, over various time periods. The concentration profile of benzene at the starting moment of emission is shown in Fig. 5(a). As expected, the initial concentration profile indicates both the location of the treating units within the WWTP and the intensity of emission of each source. Based on the findings, the ground-level concentration of benzene increases up to 48 mg/m3 over the WWTP for 10 min after the emission (Fig. 5(b)). When the time increases, benzene travels farther distances, pollutes more spaces, and reaches the maximum value of 53 mg/m3 (Fig. 5(c)). Finally, for 12 h simulation time, no significant changes are observed for the distances less than 5 km, and consequently, pollution transfer occurs under the steady-state condition up to 5 km for the simulation times of more than 12 h (Fig. 5(d)). It should be noted that the high concentration of benzene over the WWTP is a serious threat to the life of the worker or people nearby. However, its concentration is less than 1 mg/m3 for distances between 1,000 and 2,000 m far from the emission sources. In addition, its value is lower than 0.05 mg/m3 for a distance greater than 3,000 m.
Fig. 6(a)–(c) illustrate similar patterns and their simulation result for 12 h after emission for other compounds. In this regard, the sum of the concentrations of m-xylene, o-xylene, and p-xylene is considered as a single xylene substance. From the findings, the maximum concentration of toluene, ethylbenzene, and xylene over WWTP are 5.6, 6.1, and 1.4 mg/m3, respectively. In addition, their values are less than 0.01 mg/m3 for distances longer than 1,000 m far from emission points.
Due to the high emission of benzene than other BTEX compounds, greater levels of its concentration can be found in farther distances, and correspondingly its risk is greater. However, benzene concentration is less than 0.001 mg/m3 in distances greater than 4,000 m in x-direction, which is less than the recommended annual thresholds for ambient air benzene in several Asian countries [50]. Thus, pollutants cannot travel long distances from sources, and considerable ground-level concentrations are found up to 4,000 m in the x-direction.
Now, the results of the reactive and non-reactive simulations are compared to evaluate the effect of the reaction on BTEX concentration. In this regard, the conversion of BTEX compounds over different simulation times is calculated (Fig. S3). As shown, similar to the 1D simulation results, an increase in the simulation time leads to an increase in the reaction time, which, in turn, results in more consumption and conversion of the BTEX compounds. However, conversion in the 3D simulation state is greater than that of 1D. Thus, BTEX compounds have a greater chance of reacting with the fresh OH radicals because they occupy more space in the 3D state. It is worth noting that OH radical is considered a limiting reactant around the petrochemical company due to the high concentration of pollutants. As a result, the reaction effect is negligible for the smaller domains, while it is considered a significant phenomenon for the greater space domain. Furthermore, similar to the preliminary conclusion in the 1D state, xylene is affected more among the BTEX compounds, while benzene and toluene are not affected significantly by the reaction. The conversion of benzene, toluene, ethylbenzene, m-xylene, o-xylene, and p-xylene is 0.3, 1, 1.4, 5, 3.1, and 3.2%, respectively.
Finally, the model result is validated for benzene using the measured data in four different spatial points obtained from the GC method. Table 2 shows the comparison of simulation results and measured data for the same conditions in terms of temperature, wind velocity, spatial point, etc. Based on the results, there is a good agreement between simulation results and provided data, indicating the accuracy of both the Toxchem and 3D dispersion models.

### 4. Conclusions

The result of Toxchem indicates that a huge amount of around 868 kg of BTEX pollutants is released from WWTP into the environment every day, in which the Equalization basin and API unit with more than 80% emissions are considered as the major contributors to the release of BTEX compounds. As another result, based on the 3D dispersion model, the ground-level concentration of pollutants exceeds world standards for a distance of less than 2,000 m, which threatens the health of employees and surrounding residents. However, they are insignificant for distances more than 3,000 m in x-direction. Besides, the study of chemical reactions revealed that a limited number of reactions are a good representative of all reactions. In fact, there is no need to consider all of the possible reactions and species for BTEX concentration calculations. Furthermore, the reaction is not effective in closer places due to the limitation of OH concentration rather than BTEX and NO2 compounds. However, a greater space domain should be highlighted. Among the BTEX compounds, benzene accounts for about 75% of pollutants, this is not affected significantly by the reaction. Thus, special attention should be paid to avoid its adverse effects. In conclusion, the Toxchem and presented 3D dispersion model give realistic results and can be used for similar cases.

### Notes

Author Contributions

N.E. (Ph.D. student) conducted the experimental and modeling parts of the manuscript, and wrote the manuscript. E.F. (Professor) wrote and revised the manuscript. M.A.K. (Professor) revised the manuscript.

### References

1. Słomińska M, Konieczka P, Namieśnik J. The fate of BTEX compounds in ambient air. Crit Rev Environ Sci Technol. 2014;44:455–472.

2. Yu CWF, Kim JT. Building pathology, investigation of sick buildings—VOC emissions. Indoor Built Environ. 2010;19:30–39.

3. El-Hashemy MA, Ali HM. Characterization of BTEX group of VOCs and inhalation risks in indoor microenvironments at small enterprises. Sci Total Environ. 2018;645:974–983.

4. Pyta H. BTX Air Pollution in Zabrze, Poland. Pol J Environ Stud. 2006;15(5)785–791.

5. Jenkin M, Watson L, Utembe S, Shallcross D. A Common Representative Intermediates (CRI) mechanism for VOC degradation. Part 1: Gas phase mechanism development. Atmos Environ. 2008;42:7185–7195.

6. Behnami A, Benis KZ, Shakerkhatibi M, Derafshi S, Chavoshbashi MM. A systematic approach for selecting an optimal strategy for controlling VOCs emissions in a petrochemical wastewater treatment plant. Stoch Environ Res Risk Assess. 2019;33:13–29.

7. Shakerkhatibi M, Mosaferi M, Zorufchi Benis K, Akbari Z. Performance evaluation of a full-scale ABS resin manufacturing wastewater treatment plant: a case study in Tabriz Petrochemical Complex. Environ Health Eng Manage J. 2016;3:151–158.

8. Zwain HM, Vakili M, Faris AM, Dahlan I. Modeling the Fate of Phenol in Moving Bed Biofilm Reactor Sewage Treatment Plant. In : Proceedings of AICCE’19 2019; Springer International Publishing; 2020;

9. Oskouie AK, Lordi DT, Granato TC, Kollias L. Plant-specific correlations to predict the total VOC emissions from wastewater treatment plants. Atmos Environ. 2008;42:4530–4539.

10. Benis KZ, Shakerkhatibi M, Yousefi R, Kahforoushan D, Derafshi S. Emission patterns of acrylonitrile and styrene around an industrial wastewater treatment plant in Iran. Int J Environ Sci Technol. 2016;13:2353–2362.

11. Zhang C, Geng X, Wang H, Zhou L, Wang B. Emission factor for atmospheric ammonia from a typical municipal wastewater treatment plant in South China. Environ Pollut. 2017;220:963–970.

12. Calvo MJ, Prata AA, Hoinaski L, Santos JM, Stuetz RM. Sensitivity analysis of the WATER9 model: emissions of odorous compounds from passive liquid surfaces present in wastewater treatment plants. Water Sci Technol. 2018;2017:903–912.

13. Santos J, Sá L, Reis N, Gonçalves R, Siqueira R. Modelling hydrogen sulphide emission in a WWTP with UASB reactor followed by aerobic biofilters. Water Sci Technol. 2006;54:173–180.

14. Baawain M, Al-Mamun A, Omidvarborna H, Al-Jabri A. Assessment of hydrogen sulfide emission from a sewage treatment plant using AERMOD. Environ Monit Assess. 2017;189:263

15. Dinçer F, Dinçer FK, Sarı D, Ceylan Ö, Ercan Ö. Dispersion modeling and air quality measurements to evaluate the odor impact of a wastewater treatment plant in İzmir. Atmos Pollut Res. 2020;11(12)2119–2125.

16. Nemalipuri P, Das HC, Pradhan MK. Simulation of Emission from Coal-Fired Power Plant. Advances in Mechanical Engineering. Springer; 2020. p. 975–986.

17. Sagan V, Pasken R, Zarauz J, Krotkov N. SO2 trajectories in a complex terrain environment using CALPUFF dispersion model, OMI and MODIS data. Int J Environ Sci Technol. 2018;69:99–109.

18. Amoatey P, Omidvarborna H, Affum HA, Baawain M. Performance of AERMOD and CALPUFF models on SO2 and NO2 emissions for future health risk assessment in Tema Metropolis. Hum Ecol Risk Assess. 2019;25:772–786.

19. Abdul-Wahab SA. SO2 dispersion and monthly evaluation of the industrial source complex short-term (ISCST32) model at Mina Al-Fahal refinery, Sultanate of Oman. Environ Manage. 2003;31:0276–0291.

20. Fatehifar E, Elkamel A, Osalu AA, Charchi A. Developing a new model for simulation of pollution dispersion from a network of stacks. Appl Math Comput. 2008;206:662–668.

21. Xiang J, Li R, Wang G, et al. Modeling Urban PM2.5 Concentration by Combining Regression Models and Spectral Unmixing Analysis in a Region of East China. Air Qual Atmos Health. 2017;228:250

22. Ragland KW. Multiple Box Model for Dispersion of Air Pollutants from Area Sources. Atmos Environ. 1973;7:1017–1032.

23. Nourbakhsh H, Mowla D, Esmaeilzadeh F. Predicting the three dimensional distribution of gas pollutants for industrial-type geometries in the south pars gas complex using computational fluid dynamics. Ind Eng Chem Res. 2013;52:6559–6570.

24. Issakhov A, Baitureyeva A. Modeling of a passive scalar transport from thermal power plants to atmospheric boundary layer. Int J Environ Sci Technol. 2019;16:4375–4392.

25. Baker J, Walker HL, Cai X. A study of the dispersion and transport of reactive pollutants in and above street canyons—a large eddy simulation. Atmos Environ. 2004;38:6883–6892.

26. Moradpour M, Afshin H, Farhanieh B. A numerical study of reactive pollutant dispersion in street canyons with green roofs. Build Simul. 2018;11:125–138.

27. Kitabayashi K, Konishi S, Katatani A. A NOx plume Dispersion Model with Chemical Reaction in Polluted Environment. JSME Int J. 2006;49:42–47.

28. Tetzlaff G, Dlugi R, Friedrich K, et al. On Modeling Dry Deposition of Long-Lived and Chemically Reactive Species over Heterogeneous Terrain. J Atmos Chem. 2002;42:123–155.

29. Baik J-J, Kang Y-S, Kim J-J. Modeling reactive pollutant dispersion in an urban street canyon. Atmos Environ. 2007;41:934–949.

30. Joelsson LMT, Pichler C, Nilsson EJK. Tailored reduced kinetic mechanisms for atmospheric chemistry modeling. Atmos Environ. 2019;213:675–685.

31. Bright VB, Bloss WJ, Cai X. Urban street canyons: Coupling dynamics, chemistry and within-canyon chemical processing of emissions. Atmos Environ. 2013;68:127–142.

32. Fernández E, Vidal L, Canals A. Zeolite/iron oxide composite as sorbent for magnetic solid-phase extraction of benzene, toluene, ethylbenzene and xylenes from water samples prior to gas chromatography mass spectrometry. J Chromatogr A. 2016;1458:18–24.

33. Mcclenny WA, Colón M. Measurement of volatile organic compounds by the US Environmental Protection Agency Compendium Method TO-17: Evaluation of performance criteria. J Chromatogr A. 1998;813:101–111.

34. Santos JM, Kreim V, Guillot J-M, Reis NC, De Sá LM, Horan NJ. An experimental determination of the H2S overall mass transfer coefficient from quiescent surfaces at wastewater treatment plants. Atmos Environ. 2012;60:18–24.

35. Zwain HM, Nile BK, Faris AM, Vakili M, Dahlan I. Modelling of hydrogen sulfide fate and emissions in extended aeration sewage treatment plant using TOXCHEM simulations. Sci Rep. 2020;10:22209

36. Melcer H, Bell JP, Thompson DJ, Yendt CM. Modeling volatile organic contaminants’ fate in wastewater treatment plants. J Environ Eng. 1994;120:588–609.

37. Melcer H, Bell J, Thompson D, Yendt C, Kemp J, Steel P. Modeling volatile organic contaminants’ fate in wastewater treatment plants. J Environ Eng. 1994;120:588–609.

38. Monteith H, Parker W, Bell J, Melcer H. Modeling the fate of pesticides in municipal wastewater treatment. Water Environ Res. 1995;67:964–970.

39. Parker WJ, Shi J, Fendinger NJ, Monteith HD, Chandra G. Pilot plant study to assess the fate of two volatile methyl siloxane compounds during municipal wastewater treatment. Environ Toxicol Chem. 1999;18:172–181.

40. Fatehifar E, Elkamel A, Taheri M. A MATLAB-based modeling and simulation program for dispersion of multipollutants from an industrial stack for educational use in a course on air pollution control. Comput Appl Eng Educ. 2006;14:300–312.

41. Alexandrov VN, Owczarz W, Thomson P, Zlatev Z. Parallel runs of a large air pollution odel on a grid of Sun computers. Math Comput Simul. 2004;65:557–577.

42. Baker J, Walker HL, Cai X. A study of the dispersion and transport of reactive pollutants in and above street canyons—a large eddy simulation. Atmos Environ. 2004;38:6883–6892.

43. Alexandrov V, Sameh A, Siddique Y, Zlatev Z. Numerical integration of chemical ODE problems arising in air pollution models. Environ Model Assess. 1997;2:365–377.

44. Zlatev Z, Dimov I. Computational and Numerical Challenges in Environmental Modelling. 13:1st edoxford: 2006.

45. Jenkin ME, Saunders SM, Wagner V, Pilling MJ. Protocol for the development of the Master Chemical Mechanism, MCM v3 (Part B): tropospheric degradation of aromatic volatile organic compounds. Atmos Chem Phys. 2003;3:181–193.

46. Atkinson R. Our present understanding of the gas-phase atmospheric degradation of VOCs. Barnes I, Kharytonov MM, editorsSimulation and Assessment of Chemical Processes in a Multiphase Environment NATO Science for Peace and Security Series C: Environmental Security. Dordrecht: Springer; 2008. p. 1–19.

47. Seinfeld JH, Pandis SN. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd edNew York: John Wiley & Sons; 2006.

48. Jenkin ME, Saunders SM, Derwent RG, Pilling MJ. Development of a reduced speciated VOC degradation mechanism for use in ozone models. Atmos Environ. 2002;36:4725–4734.

49. Jenkin M, Khan M, Shallcross D, et al. The CRI v2. 2 reduced degradation scheme for isoprene. Atmos Environ. 2019;212:172–182.

50. Rad HD, Babaei AA, Goudarzi G, Angali KA, Ramezani Z, Mohammadi MM. Levels and sources of BTEX in ambient air of Ahvaz metropolitan city. Air Qual Atmos Health. 2014;7:515–524.

##### Fig. 1
A description of the mechanism reduction algorithm.
##### Fig. 2
The fate of entering BTEX to WWTP: (a) Benzene, (b) Toluene, (c) Ethylbenzene, and (d) Xylene.
##### Fig. 3
A comparison of the emission rate of BTEX from different treating units.
##### Fig. 4
A comparison of conversion of BTEX compounds for 1D simulation.
##### Fig. 5
Dispersion patern and concentration (mg/m3) of benzene in ground level for (a) t = 0.4 s, (b) t = 10 min, (c) t = 6 h, and (d) t = 12 h after starting emission.
##### Fig. 6
Dispersion pattern and concentration (mg/m3) of volatile compounds in ground level for (a) toluene, (b) ethylbenzene, and (c) xylene for t = 12 h.
##### Table 1
A Comparison of the Measured Values with Toxchem Results for the Concentration of BTEX Compounds
Plant outlet (ppm) Equalization outlet (ppm) API outlet (ppm) Plant inlet (ppm)

Error (%) Measured (ppm) Estimated (ppm) Error (%) Measured (ppm) Estimated (ppm) Error (%) Measured (ppm) Estimated (ppm)
0.023 trace 5.3 29.26 30.81 4.1 69.83 72.69 187.47 benzene
0.003 trace 7.9 4.26 4.59 6.4 10.32 10.98 23.19 toluene
trace trace 8.1 4.16 4.5 5.7 12.56 13.28 30.52 ethylbenzene
trace trace 8.4 1.18 1.27 8.1 2.22 2.40 4.34 xylene
##### Table 2
A Comparison between Simulation Results of Benzene Concentration and Measured Data (mg/m3)
benzene
measured simulated Error
Point (1) (X = 50, Y = 1,080) 0.051 0.047 8%
Point (2) (X = 50, Y = 890) 0.169 0.158 7%
Point (3) (X = 200, Y = 890) 14.15 13.07 8%
Point (4) (X = 200, Y = 1080) 1.67 1.48 11%
TOOLS
Full text via DOI
Supplement
E-Mail
Print
Share:
METRICS
 0 Crossref
 0 Scopus
 697 View