Optimal Production Planning and Pollution Control in Petroleum Reﬁneries Using Mathematical Programming and Dispersion Models

: Oil reﬁneries, producing a large variety of products, are considered as one of the main sources of air contaminants such as sulfur oxides (SO x ), hydrocarbons, nitrogen oxides (NO x ), and carbon dioxide (CO 2 ), which are primarily caused by fuel combustion. Gases emanated from the combustion of fuel in an oil reﬁnery need to be reduced, as it poses an environmental hazard. Several strategies can be applied in order to mitigate emissions and meet environmental regulations. This study proposes a mathematical programming model to derive the optimal pollution control strategies for an oil reﬁnery, considering various reduction options for multiple pollutants. The objective of this study is to help decision makers select the most economic pollution control strategy, while satisfying given emission reduction targets. The proposed model is tested on an industrial scale oil reﬁnery sited in North Toronto, Ontario, Canada considering emissions of NO x , SO x , and CO 2 . In this analysis, the dispersion of these air pollutants is captured using a screening model (SCREEN3) and a non-steady state CALPUFF model based on topographical and meteorological conditions. This way, the impacts of geographic location on the concentration of pollutant emissions were examined in a realistic way. The numerical experiments showed that the optimal production and pollution control plans derived from the proposed optimization model can reduce NO x , SO x , and CO 2 emission by up to 60% in exchange of up to 10.7% increase in cost. The results from the dispersion models veriﬁed that these optimal production and pollution control plans may achieve a signiﬁcant reduction in pollutant emission in a large geographic area around the reﬁnery site. whose production and emission reduction strategies are designed based on the proposed optimization model. The reported estimations of pollutant concentrations, dispersion, and transport from CALPUFF and SCREEN3 veriﬁed that the optimal production and pollution control plans derived by the proposed MINLP model may signiﬁcantly reduce the CO 2 , SO x , and NO x emission around the area of study signiﬁcantly. These results show that mathematical programming and dispersion models can be used simultaneously to derive e ﬃ cient and e ﬀ ective production and pollution control plans for oil reﬁneries.


Introduction
Global warming and the associated risks have been debated in recent decades, and climate change has been raised as a global ecological concern [1]. The United Nations Framework Convention on Climate Change (UNFCCC) [2] has set mandatory targets and timelines for greenhouse gas emissions of 38 industrial countries, including Canada, which has a 6% reduction target. Oil refineries are one of the significant sources of air contaminants including, sulfur oxides (SO x ), hydrocarbons, nitrogen oxides (NO x ), particulate matter, volatile organic compounds, and carbon dioxide (CO 2 ) [1]. Reducing gases such as CO 2 , NO x , and SO x , released from burning fuel to supply heat to different units in oil refineries, is a priority for the welfare of the society.

Figure 1.
A schematic diagram of the refinery using state-equipment network (SEN) representation (adapted from Elkamel and Al-Qahtani [28], see Appendix A for details).
The plant processes in Figure 1 are inspired by the work of Al-Qahtani and Elkamel [28,29]. They explain these plant processes in detail in their work. Therefore, for brevity, readers are referred to Al-Qahtani and Elkamel [28,29] about details of overall process flow, operation units, and the reactions and the technologies used in each unit. In Figure 1, the crude oil is processed in an atmospheric distillation unit to produce different fractions, such as liquefied petroleum gas (LPG), naphtha, kerosene, gas oil, and residues. The heavy residue is processed in a vacuum unit to produce vacuum gas oil and vacuum residues. Heavy naphtha is processed in a catalytic reformer unit to produce high octane reformates for gasoline blending. The middle distillates are processed in hydro-treating and blending units to produce jet fuels and gas oils.
The refinery model adapted in this study is for a single refinery site and is a subset of the multirefinery modeling framework presented in Al-Qahtani and Elkamel [28]. The material streams used in this model pertain to a single refinery only and include raw materials, intermediates, final products, and fuel. General mass balances are applied to all streams, along with volumetric flow rates carried out for some streams where the quality attributes blend by volume. The refinery model is formulated as a mixed integer nonlinear program (MINLP) using binary variables for selecting the ideal pollutant control strategies, which add nonlinearity to the model. In order to solve the model efficiently, non-linear terms were linearized by defining component flows instead of individual flows [30]. Further details about the model components presented below can be found in Alnahdi [31]. The plant processes in Figure 1 are inspired by the work of Al-Qahtani and Elkamel [28,29]. They explain these plant processes in detail in their work. Therefore, for brevity, readers are referred to Al-Qahtani and Elkamel [28,29] about details of overall process flow, operation units, and the reactions and the technologies used in each unit. In Figure 1, the crude oil is processed in an atmospheric distillation unit to produce different fractions, such as liquefied petroleum gas (LPG), naphtha, kerosene, gas oil, and residues. The heavy residue is processed in a vacuum unit to produce vacuum gas oil and vacuum residues. Heavy naphtha is processed in a catalytic reformer unit to produce high octane reformates for gasoline blending. The middle distillates are processed in hydro-treating and blending units to produce jet fuels and gas oils.

Optimal Production Planning Model for an Oil Refinery
The refinery model adapted in this study is for a single refinery site and is a subset of the multi-refinery modeling framework presented in Al-Qahtani and Elkamel [28]. The material streams used in this model pertain to a single refinery only and include raw materials, intermediates, final products, and fuel. General mass balances are applied to all streams, along with volumetric flow rates carried out for some streams where the quality attributes blend by volume. The refinery model is formulated as a mixed integer nonlinear program (MINLP) using binary variables for selecting the ideal pollutant control strategies, which add nonlinearity to the model. In order to solve the model efficiently, non-linear terms were linearized by defining component flows instead of individual flows [30]. Further details about the model components presented below can be found in Alnahdi [31].

Optimal Production Planning Model for an Oil Refinery
Before introducing the main mathematical model selecting both (i) the best pollution control strategies given emission reduction targets, and (ii) the best production/process plans for processing units given the level of end-product demand to be satisfied for an oil refinery, an initial model for the latter task only is presented first to improve the tractability of the proposed model for the readers. The production/process optimization model represents the refinery operations and associated material flows between process units; thus, the equations in this section are inspired by those in Al-Qahtani and Elkamel [28]. The modeling contribution of the paper lies in the way that pollution control dynamics are incorporated into the analysis as described in Section 2.2. All decision variables, variable indices, and parameter descriptions are given in the nomenclature at the end of the paper.
In order to make the operation of a refinery profitable, optimization of production plans for different intermediate and final products is required. The refinery modeling problem addressed here can be stated as 'Given product demands and quality specifications for an oil refinery, determine the production plan with minimal annual cost'.
This model aims to minimize the annualized cost, including crude oil costs (Rcost cr ) and operating costs of processing units (Ocost p ), minus the profit from the export of final products (Xpr c f r ), as expressed in Equation (1). The operating cost is assumed to be proportional to the throughput of the process and is expressed on per annum basis. Note that the revenue from the satisfied internal/local (i.e., non-export) demand is not included in the objective function because (1) a constraint is added to satisfy the internal demand level and (2) export brings significantly more revenue than satisfying the internal demand. Therefore, all feasible good solutions would have the same revenue contribution from the satisfied internal demand, which could be skipped on the objective function for the brevity of representation.
Equation (2) expresses the intermediate material balances for all oil refinery processes, where the material-balance coefficient α cr,ip,p is positive for inputs to a unit and negative for outputs. Equation (2) ensures that the intermediate streams are either consumed in the refinery fuel system (denoted with w cr,ip,r f ) or final product pool (denoted with w cr,ip, c f r ). p Z cr,p ≤ S cr , ∀cr The previous constraint is related to refinery raw material balance as the throughput is bounded by the maximum supply S cr . Basically, it limits the amount of crude oil cr to be allocated to various processing units p (i.e., Z cr,p ) based on the available crude oil supply.
The material balance for final products (i.e., x c f r ) is defined by the difference between the flow rates of intermediate streams contributing to the final product (i.e., w cr,ip, c f r ) and to the fuel system (i.e., w cr, c f r,r f ) as given in Equation (4). In Equation (10), some portions of these final products are allocated for export (e c f r ), which contributes to the objective function in Equation (1). Since a few quality attributes blend by volume, in Equation (5), the mass flow rate x c f r is converted to volumetric flow rate xv c f r by dividing it by individual specific gravity for the associated intermediate stream (i.e., sg cr,ip ). The volumetric flow rates are used in Equations (7) and (8), to specify the quality constraints for the final products.
Fuel system material balance is represented by Equation (6) using calorific value equivalent for each stream (cv ip,r f ). The fuel system may be composed of a single or a combination of different streams. In Equation (6), β cr,r f ,p denotes the consumption coefficient of fuel rf from crude stream cr in process p.
The lower and upper bounds on quality constraints for refinery products are given in Equations (7) and (8) based on products which blend by volume (q ∈ Q v ) or by mass (q ∈ Q m ). Naturally, only one of the lower quality bound multipliers q L c f r, q∈Q v and q L c f r, q∈Q m and one of the upper quality bound multipliers q U c f r, q∈Q v and q U c f r, q∈Q m are non-zero for each q.
Maximum (C max u ) and minimum (C min u ) crude oil flow rates for each processing unit are expressed in Equation (9), where zero-one coefficients γ u,p represents the assignment of unit u to an operating mode p. For instance, a reformer unit can be operated at low or high severity modes.
x c f r − e c f r ≥ D c f r , ∀c f r The final product supply from the refinery for the local market as stipulated in Equation (10) is x c f r minus the exported amount e c f r for each product. This equation ensures that the final product supply covers the local demand.
An upper and lower bound is set by the imports or resources as given in Equation (11) as per available feed-stock to the refineries. Its lower bound is useful for a situation where there is a protocol to exchange or supply oil (crude) between different countries. The upper bound may indicate limits on the feed-stock availability and refinery capacity.

Emission Control Dynamics
In this section, the selection of emission strategies is incorporated into the model for planning processing unit operations. The resulting model optimizes the overall costs of an oil refinery incurred for controlling pollutants that are generated during processing unit operations. In this study, three different mitigation techniques (n) are covered for air emission reduction: (i) fuel switching to reduce emissions Sustainability 2019, 11, 3771 7 of 31 from one type of fuel to another type of fuel (typically shifting from fuel oil to natural gas), (ii) process load shifting to adjust the production across the refinery units for reducing emissions, and (iii) implementing various air emission capture technologies. The following equation specifies emission flows based on the selection of mitigation technologies: In particular, Equation (12) formulates the emission flow rate of a production unit u ∈ U of a certain pollutant t ∈ T over multiple mitigation methods n ∈ N, where G u,t,n is a binary variable representing the selection of mitigation scheme n for process u to mitigate the emission of pollutant t.
The emission of each unit is computed as the product of the emission factor of each fuel and fuel consumption in that unit, which is related to the inlet flow rate. Hence, this formulation gives an MINLP due to the multiplication of binary and continuous variables. Thus, the above Equation (12) can be replaced by inequality constraints for the three mitigation options as shown in Constraints (13)- (16): where EF r f , CF u , and E u are the emission factor of each fuel, fuel consumption in Unit u, and an upper bound on emission, respectively. The relevant constraints can be expressed as in Equations (15) and (16): Applying a process for capture of a pollutant for a given production unit can be written as: where ε cap represents the efficiency of the capturing process. It can be noted that for each production unit, only one of the mitigation methods can be applied for each pollutant as shown in Equations (19) and (20). Equations (19) and (20) ensure that up to one fuel switching option and one emission capturing technology is used, respectively. Alternately, it can be specified that selection of one of the fuel switching options is mutually exclusive with the selection of one of the air emission capturing technologies, as given in Equation (21), which guarantees that only one technology, in total, may be selected. n∈switch G u,t,n ≤ 1, ∀u, t n∈cap G u,t,n ≤ 1, ∀u, t Once the cost of emission control is included, nonlinear terms may appear in the objective function due to the multiplication of process emissions E u,t with binary variable G u,t,n . This is prevented by introducing additional variables and a new set of constraints.
Equations (22) and (23) express the cost terms for emission capturing by defining a set of bounds for a continuous variable φ u,t defined to represent the annual cost of emission capturing in the objective function, where HC − u,t and HC + u,t are properly selected for the cost, cost n∈cap E u,t . Below a similar method is used to express the cost associated with fuel switching.
Equations (24) and (25) express the cost of fuel switching by defining a set of bounds for a continuous variable ϕ u,t , defined to represent the annual cost of fuel switching in the objective function, where lower and upper bounds HS − u,t and HS + u,t are properly selected for the cost, cost n∈switch E u,t . After adding these set of costs, the objective function may be reformulated as shown in Equation (26). The main model for selecting both the best emission control strategies and operation plans for processing units can be formulated as follows: The main model was coded and solved using GAMS [30] software and the computations were performed on a Pentium 4, 3.0 GHz processor.

Air Pollution Dispersion Models
This section provides the details of the two software models used to describe how air pollutants, emitted by a source, disperse in the ambient atmosphere. CALPUFF and SCREEN3 are used as modeling tools to estimate the overall concentration of SO x , NO x , and CO 2 within the area of study.

CALPUFF Dispersion Model
CALPUFF is a computer-based tool for air dispersion modeling which has been developed by the United States Environmental Protection Agency [10]. It consists of a meteorological, non-steady-state puff dispersion and post-processing modules, and it simulates the effects of temporally and spatially changing meteorological conditions on air pollutant movement. CALPUFF can be used for modeling areas that are 5 km to 300 km away from the source. CALPUFF makes provision for point, area, line, and volume sources and assesses the mesoscale transport of pollutants as well as their dispersion in the surrounding complex terrain. For instance, the puffs emitted from a stack point are modeled individually based on conditions of wind direction and speed on an hourly basis, to estimate the pollutant concentrations.
As shown in Figure 2, the CALPUFF modeling system is divided into three modules namely: CALMET, CALPUFF, and CALPOST. CALMET is a meteorological model requiring inputs of terrain elevation, land use, atmospheric temperature, wind velocity, cloud cover, relative humidity, and atmospheric pressure (i.e., surface, upper air, precipitation, and over water) [10]. The purpose of CALMET is to generate 3D wind fields that are used in CALPUFF and CALPOST modules. CALPUFF is a Gaussian puff model with different effects such as chemical removal, wet and dry deposition, and complex terrain algorithms. CALPOST is a post-processing package to process the outputs generated by CALMET and CALPUFF for plotting on modeling domain maps [10]. CALMET is to generate 3D wind fields that are used in CALPUFF and CALPOST modules. CALPUFF is a Gaussian puff model with different effects such as chemical removal, wet and dry deposition, and complex terrain algorithms. CALPOST is a post-processing package to process the outputs generated by CALMET and CALPUFF for plotting on modeling domain maps [10].

SCREEN3 Dispersion Model
Screen View version 3.0 (SCREEN3) is among the programs recommended by the United States Environmental Protection Agency and is freely available. SCREEN3 is designed to provide a simple way to estimate the concentration of pollutants based on simple tracking information available to a large number of users. SCREEN3 uses a Gaussian plume model, taking into account meteorological factors to calculate the concentration and dispersion of contaminants from stationary sources. The model examines a number of classes of stability conditions and assumes stable wind speed to identify the turbulence of the atmosphere, which greatly influences the spread of pollutants. Stability conditions are divided into six classes: class A-extremely unstable, class B-unstable, class Cslightly unstable, class D-neutral, class E-slightly stable, and class F-stable. The stability condition occurs when there is an absence of solar radiation, no clouds, and the presence of mild winds. This condition is less favorable for the dispersion of pollutants. Mixing time is a measure of the atmosphere being at higher process turbulence, thus favoring the dispersion of pollutants.
One of the limitations of the Screen model is that it is unable to determine the impacts from multiple sources and merge them into a single representation. To make such a representation, it is necessary to make separate simulations, and then manually superimpose the results obtained. Another limitation is that the program presents the results only linearly and at a maximum distance of 50 km. These limitations can be overcome by using other commercial modeling tools such as Industrial Source Complex Short Term Version 3 (ISCST3).

Case Study Settings and Results for Production Planning and Emission Reduction in an Oil Refinery
The oil refinery flow diagram considered in the case studies is shown in Figure 1, and the major unit capacity limits are shown in Table 1. The planning period considered is one year, and the

SCREEN3 Dispersion Model
Screen View version 3.0 (SCREEN3) is among the programs recommended by the United States Environmental Protection Agency and is freely available. SCREEN3 is designed to provide a simple way to estimate the concentration of pollutants based on simple tracking information available to a large number of users. SCREEN3 uses a Gaussian plume model, taking into account meteorological factors to calculate the concentration and dispersion of contaminants from stationary sources. The model examines a number of classes of stability conditions and assumes stable wind speed to identify the turbulence of the atmosphere, which greatly influences the spread of pollutants. Stability conditions are divided into six classes: class A-extremely unstable, class B-unstable, class C-slightly unstable, class D-neutral, class E-slightly stable, and class F-stable. The stability condition occurs when there is an absence of solar radiation, no clouds, and the presence of mild winds. This condition is less favorable for the dispersion of pollutants. Mixing time is a measure of the atmosphere being at higher process turbulence, thus favoring the dispersion of pollutants.
One of the limitations of the Screen model is that it is unable to determine the impacts from multiple sources and merge them into a single representation. To make such a representation, it is necessary to make separate simulations, and then manually superimpose the results obtained. Another limitation is that the program presents the results only linearly and at a maximum distance of 50 km. These limitations can be overcome by using other commercial modeling tools such as Industrial Source Complex Short Term Version 3 (ISCST3).

Case Study Settings and Results for Production Planning and Emission Reduction in an Oil Refinery
The oil refinery flow diagram considered in the case studies is shown in Figure 1, and the major unit capacity limits are shown in Table 1. The planning period considered is one year, and the feedstock to the refinery is of a single type (e.g., Arabian light crude) with a total flow rate of 12,000 kton/y in order to produce different final products. These parameters are inspired by the work of Al-Qahtani [32], and further details about parametrization and the following results are available in Alnahdi [31]. The pollutants considered in this study are SO x and NO x , whereas CO 2 emission is incorporated into the analysis as a greenhouse effect. Several emission mitigation alternatives are considered. Among these, fuel switching represents switching from the current fuel (fuel oil) to natural gas. The SO x capture process considered in this study is wet scrubber (WS) with fuel-gas desulphurization (FGD) technology, and the considered NO x capture process is based on retrofitting current burners with low NO x burner technology (will be denoted as LNB hereafter). For the CO 2 capture process, MEA absorption technology is considered. CO 2 emissions parameters are estimated based on the work by Ritter et al. [33], whereas emissions data of the pollutants are based on Kassinis [34]. Cost parameters associated with fuel switching and MEA capture process are estimated based on the work by Elkamel et al. [35]; while, those associated with SO x and NO x capture processes are taken from a report by the World Bank [36].
Three scenarios are considered in this study in order to illustrate the validity of the model discussed in the previous section. For the first scenario, the process planning model in Section 2.1 is solved without any emissions mitigation, which is considered as the base case scenario. The second scenario is related to reducing one particular emission at a time (e.g., CO 2 ) without considering other emissions (e.g., SO x and NO x ) using a special case of the model in Section 2.2. That is, each pollutant is reduced in an independent manner. Finally, the third scenario is related to the examination of reducing different emissions together in a general/dependent manner using the general model in Section 2.2.

Results of Base Case Scenario
The objective of this base case scenario is to minimize the overall annualized cost minus the revenue from export while meeting the internal demand for each product with quality specifications. The market demands and specifications for different products that the refinery has to meet are shown in Table 2, and this is applied for all case studies. The emission results of the base case scenario, where no emission reduction plans are considered, are shown in Figure 3. These results imply a total annual production cost minus export revenue of $3,295,058 for the considered refinery. The total emissions from all the units for SO x , NO x , and CO 2 were 8170.6, 2826.9 and 1342.3 kton/y, respectively.

Results for Independent Emission Reduction
In this scenario, the focus is on reducing one emission at a time without considering the other pollutants, in order to analyze the impact of reducing each emission independently and to figure out the overall refinery cost correspondingly. Comparisons of cost increment when reducing SOx, NOx, and CO2 emissions are given in Tables 3-5, respectively.

Results for Independent Emission Reduction
In this scenario, the focus is on reducing one emission at a time without considering the other pollutants, in order to analyze the impact of reducing each emission independently and to figure out the overall refinery cost correspondingly. Comparisons of cost increment when reducing SO x , NO x , and CO 2 emissions are given in Tables 3-5, respectively. For the SO x reduction scenario, it has been observed that with only a 6.6% increase in cost, almost 60% of SO x releases can be alleviated by installing a wet scrubber capture process (WS) and/or through switching of fuels. On the other hand, NO x discharge may be reduced by retrofitting the existing burners with LNB. There may be a maximum emissions reduction of 60% at a cost increment of 6.9% according to Table 4. The NO x release is higher when no fuel switching was selected over the present fuel. Table 5 shows that a reduction of 60% of CO 2 emissions can be achieved in exchange for a 6.7% increase in cost. Furthermore, it has been observed that a decrease in the CO 2 discharge from 25% to 10% has been achieved without shifting units to the natural gas or by setting up extra MEA capturing practices. This decrease has been attained through load shifting, which is the recommended option for 10% emissions reduction target with a minor increase in cost. These scenarios demonstrate the flexibility of the model for proposing diverse emission reduction actions with different targets.

Results for Simultaneous Emissions Reduction
The goal here is to mitigate SO 2 , NO 2 , and CO 2 emissions simultaneously to various levels, and to compare the overall costs of the resulting mitigation plan. Table 6 provides a summary of selected scenarios when varying reduction targets for all emissions including SO 2 , NO 2 , and CO 2 and the corresponding annual cost for each plan. For instance, reducing SO 2 , NO 2 , and CO 2 by 10%, 30%, and 60% each, results in increasing the annual cost by 3.2%, 9.2%, and 10.7%, respectively. Figures 4a and 5b show the annual cost when reducing SO 2 from 0% to 60% with varying CO 2 and NO 2 reduction levels. These figures show that overall cost increases almost linearly with the emissions reduction level. 5b show the annual cost when reducing SO2 from 0% to 60% with varying CO2 and NO2 reduction levels. These figures show that overall cost increases almost linearly with the emissions reduction level.
(a) (b)    5b show the annual cost when reducing SO2 from 0% to 60% with varying CO2 and NO2 reduction levels. These figures show that overall cost increases almost linearly with the emissions reduction level.

Case Study Details and Results for Air Pollution Dispersion
In this section, CALPUFF and SCREEN3 modeling tools are used to estimate the overall concentration of SO x , NO x , and CO 2 dispersed from an oil refinery site located in North Toronto, Ontario-Canada. This refinery has a similar process design to the one illustrated in Figure 1. The location of this refinery is 50 km away from the urban Toronto; thus, there is a need for a screening analysis gauging to what extent the emission from the refinery under the proposed mitigation strategies would affect the surrounding area.
The first step in air pollution analysis with CALPUFF model is identifying the meteorological domain information for the case study region, which is given in Table 7. The surface data for the Toronto area were acquired from the weather records in the Canadian government website [37]. The surface stations were chosen based on their proximity to the source point and upper air stations. For each station, the hourly data includes the date, time, temperature, wind speed and direction, ceiling height, cloud cover, and station pressure. The hourly data for two modeling periods from (i) 1 January 2014 at 00:00 h to 31 January 2014 at 23:00 h; (ii) 1 May 2014 at 00:00 h to 31 May 2014 at 23:00 h were extracted and organized in a certain layout that is suitable for use in 'SMERGE' to create a formatted file 'SURF.DAT', which is compatible for usage with CALMET. The relevant surface station data is shown in Table 8.  The meteorological data of the upper air for the location was obtained from the radiosonde station records in the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (NOAA/ESRL) radiosonde database [38]. These data records contain the station ID number, date and time, and information of the sounding level followed by pressure, temperature, elevation, and wind direction and speed for each sounding level. The hourly data for the Toronto area was taken from one radiosonde station that is close to the study site for the two modeling periods mentioned above and was then prepared in a format suitable to use in "READ62" to generate the "UP.DAT" file that was used in the CALMET program. The relevant information about the radiosonde station is shown in Table 9. The geophysical data such as land use and terrain were obtained from the website of Geographic Information Systems [39], and used as input files in CTG-PROC and TERREL to produce "LU.DAT" and "TERREL.DAT", respectively. All this data is compressed by a MAKEGEO program to generate the output file 'GEO.DAT', which was used in the CALMET program. The proposed mathematical model in Section 2.2 was used to find the pollutant emission rates for the case study. The source parameters for the case study are shown in Table 10. For different reduction plans, the emission rates of the three pollutants are given in Table 11, which were used in the CALPUFF model and specified for the two modeling periods. Initially, the CALPUFF model was used to estimate the pollutant concentrations emanated from the oil refinery over the surrounding region for a two-month modeling period of January and May of 2014. The CALPOST post-processor was used to determine their spatial distribution. Figure 6 presents the characteristics of the study area along with the ability of the CALPUFF model to simulate the geographical condition of the area of interest. Table 12 illustrates the maximum and minimum monthly average concentrations of CO 2 for January 2014.
Initially, the CALPUFF model was used to estimate the pollutant concentrations emanated from the oil refinery over the surrounding region for a two-month modeling period of January and May of 2014. The CALPOST post-processor was used to determine their spatial distribution. Figure 6 presents the characteristics of the study area along with the ability of the CALPUFF model to simulate the geographical condition of the area of interest. Table 12 illustrates the maximum and minimum monthly average concentrations of CO2 for January 2014.    Figure 7 shows the CO 2 plume distributed significantly in the Toronto area in January 2014. The highest monthly CO 2 concentration calculated by CALPUFF view model is 216 µg/m 3 and the lowest is 2 µg/m 3 when considering the base case scenario as shown in Figure 7a. The plume affecting Toronto residential area concentration is 20 µg/m 3 . For the 10% reduction target, the maximum monthly average concentration of CO 2 was 194 µg/m 3 , and the minimum average concentration was 2 µg/m 3 , as seen in Figure 7b. The pollutant is dispersed in all the directions, especially in the northwest direction as seen in Figure 7c. For the 50% reduction target of CO 2 , the maximum monthly concentration was 127 µg/m 3 at a distance of 2 km from the source, as shown in Figure 7e. The dispersion of CO 2 was heading significantly in the south and northeast direction as presented in Figure 7d, thus affecting people in that area. As displayed in Figure 7f, the plume dispersed significantly towards the northwest and it shows that the maximum monthly average concentration of CO 2 obtained for 60% reduction target reached 115 µg/m 3 within 2 km northeast of the source. It is clear how the plumes are covering the Toronto residential area with the plumes concentration ranging from 216 to 1 µg/m 3 represented by the color code. Figure 7 shows the CO2 plume distributed significantly in the Toronto area in January 2014. The highest monthly CO2 concentration calculated by CALPUFF view model is 216 µ g/m 3 and the lowest is 2 µ g/m 3 when considering the base case scenario as shown in Figure 7a. The plume affecting Toronto residential area concentration is 20 µ g/m 3 . For the 10% reduction target, the maximum monthly average concentration of CO2 was 194 µ g/m 3 , and the minimum average concentration was 2 µ g/m 3 , as seen in Figure 7b. The pollutant is dispersed in all the directions, especially in the northwest direction as seen in Figure 7c. For the 50% reduction target of CO2, the maximum monthly concentration was 127 µ g/m 3 at a distance of 2 km from the source, as shown in Figure 7e. The dispersion of CO2 was heading significantly in the south and northeast direction as presented in Figure 7d, thus affecting people in that area. As displayed in Figure  7f, the plume dispersed significantly towards the northwest and it shows that the maximum monthly average concentration of CO2 obtained for 60% reduction target reached 115 µ g/m 3 within 2 km northeast of the source. It is clear how the plumes are covering the Toronto residential area with the plumes concentration ranging from 216 to 1 µ g/m 3 represented by the color code. For the second period of May 2014, Table 13 shows the maximum and minimum monthly average concentrations of CO2 for the six reduction plans. Figure 8 shows the dispersion of CO2 emission to the surrounding area of Toronto for different mitigation plans. The maximum monthly CO2 concentration calculated by the CALPUFF model for the base case scenario is 175 µ g/m 3 and the lowest is 4 µ g/m 3 . The plume affecting the Toronto residential area has a concentration of 20 µ g/m 3 , as shown in Figure 8a. The maximum CO2 monthly average concentration became 158 µ g/m 3 once the oil refinery site reduces its emission by 10%, as shown in Figure 8b. For the case of 30% CO2 reduction, Figure 8c shows that the maximum average concentration reduced to 128 µ g/m 3 and the minimum concentration became 3 µ g/m 3 . Furthermore, Figure 8f shows that reducing the CO2 emission by 60% results in reducing the monthly average concentration to 93.1 µ g/m 3 . The maximum For the second period of May 2014, Table 13 shows the maximum and minimum monthly average concentrations of CO 2 for the six reduction plans. Figure 8 shows the dispersion of CO 2 emission to the surrounding area of Toronto for different mitigation plans. The maximum monthly CO 2 concentration calculated by the CALPUFF model for the base case scenario is 175 µg/m 3 and the lowest is 4 µg/m 3 . The plume affecting the Toronto residential area has a concentration of 20 µg/m 3 , as shown in Figure 8a.
The maximum CO 2 monthly average concentration became 158 µg/m 3 once the oil refinery site reduces its emission by 10%, as shown in Figure 8b. For the case of 30% CO 2 reduction, Figure 8c shows that the maximum average concentration reduced to 128 µg/m 3 and the minimum concentration became 3 µg/m 3 . Furthermore, Figure 8f shows that reducing the CO 2 emission by 60% results in reducing the monthly average concentration to 93.1 µg/m 3 . The maximum monthly average concentration of CO 2 for all reduction plans that affect the Toronto area is between 2 and 30 µg/m 3 . No health symptoms are associated with this range of concentration values according to air quality standard and guideline [40]. Table 13. The range of average CO 2 concentrations in May 2014 under the optimal production and emission control strategies for various targeted emission reduction levels. Similarly, the following set of results are related to the dispersion of SO x and NO x for the periods of January and May 2014, the tables and figures for which are excluded here for brevity and can be found in Appendix B. For January 2014, the maximum monthly average concentration of SO 2 when no reduction plan is applied to the refinery was 1312 µg/m 3 and the plumes dispersed 5 km northeast of the source. The maximum and minimum average monthly concentrations of SO 2 were 1181 and 12 µg/m 3 for 10% reduction target. It was observed that when the pollutants from the source were reduced by certain percentages, the concentration of the pollutants in the receptor area reduced accordingly. Therefore, reducing the SO 2 by 30%, 50%, and 60% results in reducing the maximum monthly average concentration to 956.5, 775, and 697 µg/m 3 , respectively. For May 2014, the maximum and minimum monthly average concentrations of SO 2 are 1066 and 22 µg/m 3 for the base case scenario, whereas the plume affecting the Toronto residential area has a concentration of 100 µg/m 3 . Reducing 10% of SO 2 in the oil refinery decreases the maximum and minimum monthly average concentration of SO 2 to 960 and 20 µg/m 3 , respectively. The highest monthly SO 2 concentration computed by CALPUFF when reducing the SO 2 emission by 30% is 777 µg/m 3 and the lowest is 16 µg/m 3 . The monthly concentration of plume covering the Toronto residential area is 80 µg/m 3 . reduction, Figure 8c shows that the maximum average concentration reduced to 128 µ g/m 3 and the minimum concentration became 3 µ g/m 3 . Furthermore, Figure 8f shows that reducing the CO2 emission by 60% results in reducing the monthly average concentration to 93.1 µ g/m 3 . The maximum monthly average concentration of CO2 for all reduction plans that affect the Toronto area is between 2 and 30 µ g/m 3 . No health symptoms are associated with this range of concentration values according to air quality standard and guideline [40].  Similarly, the following set of results are related to the dispersion of SOx and NOx for the periods of January and May 2014, the tables and figures for which are excluded here for brevity and can be found in Appendix B. For January 2014, the maximum monthly average concentration of SO2 when no reduction plan is applied to the refinery was 1312 µ g/m 3 and the plumes dispersed 5 km northeast of the source. The maximum and minimum average monthly concentrations of SO2 were 1181 and 12 µ g/m 3 for 10% reduction target. It was observed that when the pollutants from the source were reduced by certain percentages, the concentration of the pollutants in the receptor area reduced accordingly. Therefore, reducing the SO2 by 30%, 50%, and 60% results in reducing the maximum For January 2014, the maximum NO 2 concentration for the oil refinery was about 454 µg/m 3 when no reduction plan was applied. For the 10% reduction plan, it was observed that NO 2 was dispersed from the northwest to the northeast of the Toronto area with maximum and minimum monthly average concentrations of 409 and 4 µg/m 3 , respectively. The maximum monthly average concentration of NO 2 for the plant was about 331 µg/m 3 when reducing the emission by 30%. The monthly average of SO 2 concentration for 40%, 50%, and 60% reduction plans was 298, 268, and 241 µg/m 3 , respectively. For May 2014, the maximum monthly average concentration of NO 2 was 332 µg/m 3 for the 30% reduction plan.

Reduction Plan
For SCREEN3, the input data used is for the atmospheric emission rates of SO x , NO x, and CO 2 from the oil refinery. The results of SCREEN3 are expressed in concentration units (µg/m 3 ) and distance (m). The simulations were performed considering the option of 'full meteorology'; that is, defining the type of atmospheric stability class where the program assumes the 'C' class was omitted. The maximum concentration was calculated for SO x for the base case and various mitigation plans. In all scenarios, the maximum concentration was found at 1200 m away from the source and the emission concentrations were 1.6 × 10 5 , 1.3 × 10 5 , 9.5 × 10 4 , and 6.4 × 10 4 µg/m 3 for the base case, 20%, 40%, and 60% reduction plans, respectively.

Conclusions
This study addresses the problem of selecting the best pollution control strategies for an oil refinery given specific values of emission reduction targets. The problem has been formulated MINLP based on three mitigation options and chooses the optimal set to meet a certain emission reduction goal with the minimum annual cost minus export revenues, while ensuring the satisfaction of the demand levels and quality specifications. The model was illustrated on an industrial scale refinery case study, considering three pollutants (SO x , NO x , and CO 2 ) with different mitigation alternatives of fuel switching and capturing pollutants. The results showed that with only a 6.6% increment in price, almost 60% of SO x releases can be alleviated by installing a wet scrubber capture process WS and through switching of fuels. Furthermore, reducing NO x by 60% increases the cost by 6.9%. Our results showed that a reduction of 60% of CO 2 emissions can be gained in exchange for a 6.7% increase in the total cost. Nevertheless, when three pollutants considered together, reducing SO x , CO x , and NO x emission by 60% simultaneously requires a 10.7% increase in cost.
This comparison illustrates that emission rates can be reduced with reasonable cost using the aforementioned strategies in an industrial-scale facility. However, emission reduction strategies should be analyzed collectively for all considered pollutants. This is because the strategies for individual pollutants may not be the same and the strategy for one pollutant may not achieve the desired effect on the others. As indicated our results, decreasing the emission of all three pollutants would be close to 40% (e.g., (10.7−6.6%)/10.7%) more expensive. Therefore, relying on an independent analysis of each pollutant could be quite misleading.
Furthermore, two air dispersion models are used to investigate the dispersion of the pollutants released from a potential oil refinery located in North Toronto, Ontario, Canada, whose production and emission reduction strategies are designed based on the proposed optimization model. The reported estimations of pollutant concentrations, dispersion, and transport from CALPUFF and SCREEN3 verified that the optimal production and pollution control plans derived by the proposed MINLP model may significantly reduce the CO 2 , SO x , and NO x emission around the area of study significantly. These results show that mathematical programming and dispersion models can be used simultaneously to derive efficient and effective production and pollution control plans for oil refineries.
applied to the refinery was 1312 µ g/m 3 as seen in Figure A2a and the plumes dispersed 5 km northeast of the source. Figure A2b illustrates that the maximum and minimum average monthly concentrations of SO2 were 1181 and 12 µ g/m 3 for 10% reduction target. One can clearly notice that when reducing the pollutants from the source by certain percentages, the concentration of the pollutants in the receptor area will reduce accordingly. Therefore, reducing the SO2 by 30%, 50%, and 60% results in reducing the maximum monthly average concentration to 956.5, 775, and 697 µ g/m 3 , respectively, as shown in Figure A2c,e,f.   Figure A3 shows the SO2 overall plume dispersion for May 2014. For the base case scenario, the maximum and minimum average monthly concentrations of SO2 as shown in Figure A3a is 1066 and 22 µ g/m 3 , whereas the plume affecting the Toronto residential area has a concentration of 100 µ g/m 3 . Reducing 10% of SO2 in the oil refinery decreases the maximum and minimum monthly average concentration of SO2 to 960 and 20 µ g/m 3 , respectively, as seen in Figure A2b. From Figure A3c, the highest monthly SO2 concentration computed by CALPUFF when reducing the SO2 emission by 30% is 777 µ g/m 3 and the lowest is 16 µ g/m 3 and the monthly concentration of plume covering the Toronto residential area is 80 µ g/m 3 . Figure A3c,e,f illustrate the maximum and minimum average monthly concentrations of SO2 for other reduction plans of 40%, 50%, and 60%, respectively.   Figure A3 shows the SO 2 overall plume dispersion for May 2014. For the base case scenario, the maximum and minimum average monthly concentrations of SO 2 as shown in Figure A3a is 1066 and 22 µg/m 3 , whereas the plume affecting the Toronto residential area has a concentration of 100 µg/m 3 . Reducing 10% of SO 2 in the oil refinery decreases the maximum and minimum monthly average concentration of SO 2 to 960 and 20 µg/m 3 , respectively, as seen in Figure A2b. From Figure A3c, the highest monthly SO 2 concentration computed by CALPUFF when reducing the SO 2 emission by 30% is 777 µg/m 3 and the lowest is 16 µg/m 3 and the monthly concentration of plume covering the Toronto residential area is 80 µg/m 3 . Figure A3c,e,f illustrate the maximum and minimum average monthly concentrations of SO 2 for other reduction plans of 40%, 50%, and 60%, respectively. different reduction plans for May 2014. Figure A3 shows the SO2 overall plume dispersion for May 2014. For the base case scenario, the maximum and minimum average monthly concentrations of SO2 as shown in Figure A3a is 1066 and 22 µ g/m 3 , whereas the plume affecting the Toronto residential area has a concentration of 100 µ g/m 3 . Reducing 10% of SO2 in the oil refinery decreases the maximum and minimum monthly average concentration of SO2 to 960 and 20 µ g/m 3 , respectively, as seen in Figure A2b. From Figure A3c, the highest monthly SO2 concentration computed by CALPUFF when reducing the SO2 emission by 30% is 777 µ g/m 3 and the lowest is 16 µ g/m 3 and the monthly concentration of plume covering the Toronto residential area is 80 µ g/m 3 . Figure A3c,e,f illustrate the maximum and minimum average monthly concentrations of SO2 for other reduction plans of 40%, 50%, and 60%, respectively.   Table A3 presents the maximum and minimum monthly average concentrations of NO2 for different reduction plans in January 2014. Figure A4 shows the typical NO2 dispersion of a one-month average from the stack of the oil refinery for the base case scenario and various reduction plans. The maximum NO2 concentration for the month of January for the oil refinery was about 454 µ g/m 3 when no reduction plan is applied, as shown in Figure 11a. From Figure A4b, NO2 was dispersed from the northwest to the northeast of the Toronto area with maximum and minimum average monthly concentrations of 409 and 4 µ g/m 3 , respectively, when applying the 10% reduction plan. The maximum monthly average concentration of SO2 for the plant was about 331 µ g/m 3 when reducing the emission by 30%, as shown in Figure 11c. The monthly average of SO2 concentration for 40%, 50%, and 60% reduction plans was 298, 268, and 241 µ g/m 3 , respectively, as illustrated in Figure A4d-f, respectively.   Table A3 presents the maximum and minimum monthly average concentrations of NO 2 for different reduction plans in January 2014. Figure A4 shows the typical NO 2 dispersion of a one-month average from the stack of the oil refinery for the base case scenario and various reduction plans. The maximum NO 2 concentration for the month of January for the oil refinery was about 454 µg/m 3 when no reduction plan is applied, as shown in Figure A3a. From Figure A4b, NO 2 was dispersed from the northwest to the northeast of the Toronto area with maximum and minimum average monthly concentrations of 409 and 4 µg/m 3 , respectively, when applying the 10% reduction plan. The maximum monthly average concentration of SO 2 for the plant was about 331 µg/m 3 when reducing the emission by 30%, as shown in Figure A3c. The monthly average of SO 2 concentration for 40%, 50%, and 60% reduction plans was 298, 268, and 241 µg/m 3 , respectively, as illustrated in Figure A4d-f, respectively.  For the month of May in 2014, Table A4 shows the maximum and minimum monthly average concentrations of NO2 for the base case scenario and all reduction plans. Figure A5 shows the typical NO2 dispersion of one month (May) from the oil refinery for the base case scenario and all reduction plans. The typical NO2 dispersion of 1-month average in the Toronto area was shown in Figure A5b when applying the 10% reduction plan. The maximum monthly average concentration of NO2 was Figure A4. NO x dispersion around North Toronto for January 2014: (a) base case; (b) 10% reduction; (c) 30% reduction; (d) 40% reduction; (e) 50% reduction; (f) 60% reduction. Table A3. The range of average NO 2 concentrations in January 2014 under the optimal production and emission control strategies for various targeted emission reduction levels. For the month of May in 2014, Table A4 shows the maximum and minimum monthly average concentrations of NO 2 for the base case scenario and all reduction plans. Figure A5 shows the typical NO 2 dispersion of one month (May) from the oil refinery for the base case scenario and all reduction plans. The typical NO 2 dispersion of 1-month average in the Toronto area was shown in Figure A5b when applying the 10% reduction plan. The maximum monthly average concentration of NO 2 was 332 µg/m 3 for the 30% reduction plan (see Figure A5c). Figure A5d-f, illustrate the maximum and minimum average monthly concentration of SO 2 for other reduction plans of 40%, 50%, and 60%. 332 µ g/m 3 for the 30% reduction plan (see Figure A5c). Figure A5d-f, illustrate the maximum and minimum average monthly concentration of SO2 for other reduction plans of 40%, 50%, and 60%.  For SCREEN3, the necessary input data includes the atmospheric emission rates of SOx, NOx, and CO2 from the oil refinery. The results of SCREEN3 are expressed in concentration units (µ g/m 3 ) and distance (m). The first group of simulations was performed considering the option of Full  For SCREEN3, the necessary input data includes the atmospheric emission rates of SO x , NO x , and CO 2 from the oil refinery. The results of SCREEN3 are expressed in concentration units (µg/m 3 ) and distance (m). The first group of simulations was performed considering the option of Full Meteorology. That is, defining the type of atmospheric stability class, where the program assumes the 'C' class, was omitted. Figure A6 shows the maximum concentration calculated for SO x for the base case and various mitigation plans. The figure depicts that the maximum concentration was found 1200 m away from the source for all scenarios and the emission concentrations were 1.6 × 10 5 , 1.3 × 10 5 , 9.5 × 10 4 and 6.4 × 10 4 µg/m 3 for the base case, 20%, 40%, and 60% reduction plans, respectively. The second set of simulations was performed for each type of atmospheric stability. Figure A7 shows the concentration of emissions of SO 2 for Type A, D, and F stability classes, for base case scenario and 40% reduction plans.

Reduction Plan
Meteorology. That is, defining the type of atmospheric stability class, where the program assumes the 'C' class, was omitted. Figure A6 shows the maximum concentration calculated for SOx for the base case and various mitigation plans. The figure depicts that the maximum concentration was found 1200 m away from the source for all scenarios and the emission concentrations were 1.6 × 10 5 , 1.3 × 10 5 , 9.5 × 10 4 and 6.4 × 10 4 µ g/m 3 for the base case, 20%, 40%, and 60% reduction plans, respectively. The second set of simulations was performed for each type of atmospheric stability. Figure A7 shows the concentration of emissions of SO2 for Type A, D, and F stability classes, for base case scenario and 40% reduction plans.