Interactive comment on “ Soot and SO 2 contribution to the supersites in the MILAGRO campaign from elevated flares in the Tula Refinery ”

This paper describes an investigation of the contribution of flare emissions to local and regional air quality in Mexico. The relevance of this study extends beyond this specific application since flare emissions are poorly understood and contribute significant uncertainty to air quality modeling in areas with intensive petrochemical operations. The introduction does a good job describing the relevance of this work, reviewing previous studies on flame emissions, and giving background on the Tula area and its emission activities. While Sections 1 and 2 read well and are clearly presented the later results sections are rather disjointed, rambling, and difficult to follow. Also, grammatical


Introduction
Flaring is an important source of greenhouse gases, particulate matter and Volatile Organic Compounds (VOCs) in both upstream and downstream operations in the oil and gas industry.Worldwide, about 134 billion cubic meters (bcm) of gas were flared in 2010 according to the World Bank Global Gas Flaring Reduction Initiative (GGFR) estimates of flared volume from satellite data (GGFR, 2012).This represents up to 400 million metric tons of CO 2 eq global greenhouse gas emissions.In 2010 Mexico was the eleventh emitting country with 2.5 bcm (World Bank, 2012).The US Energy Information Administration (EIA) estimates about 17 million metric tons of CO 2 eq emitted from flaring of natural gas for Mexico in 2009 (US EIA, 2012).
The composition of the gas to be flared depends on process operations.It may contain heavy hydrocarbons, water vapor, hydrogen sulfide, nitrogen, and carbon dioxide principally (GE Energy, 2011).As a result, harmful by-products such as NO x , reactive organic gases (ROG), polycyclic aromatic hydrocarbons (PAHs), sulfur dioxide and soot are emitted to the atmosphere.
Black carbon (a component of soot) emissions from flaring facilities are of particular interest because of its harmful effect to air quality and its contribution to both global and regional climate change (Chung and Seinfeld, 2005).It warms the atmosphere, reduces surface albedo on ice surfaces accelerating snow melt, affects clouds, slows convection and impacts the hydrological cycle (Skeie et al., 2011;McMeeking et al., 2011;Kahnert et al., 2011).In terms of global warming potential, black carbon is second to CO 2 (Moffet and Prather, V. H. Almanza et al.: Soot and SO 2 contribution to the supersites in the MILAGRO campaign 2009); since its atmospheric lifetime is short (days to weeks) with respect to CO 2 , it is considered a short-lived climate forcer (SLCF) (UNEP, 2011).This feature promotes a rapid climate response when changing emissions.Thus, control of black carbon particles by emission reduction measures would have immediate benefits on air quality and global warming (Bond, et. al, 2011;UNEP, 2011).Therefore, having robust estimates of emission rates of soot can be helpful in improving emission inventories for air quality modeling.
In Mexico, flaring operations are basically associated with exploration, exploitation, refining and petrochemical operations of PEMEX, the national oil company.Geographically, these are located in the eastern coast along the Gulf of Mexico (Tamaulipas, Veracruz and Tabasco); in the southeastern part of the country (the Campeche sound), and in the central region of the country (Guanajuato, Hidalgo and Oaxaca states).In this work we focus on flaring emissions from the Tula refinery (Hidalgo State) during MILAGRO campaign (Molina et al., 2010).The motivation to focus on flaring operations is twofold: (i) to infer the contribution of the emissions by elevated flares on pollution levels; (ii) to assess the viability to further extend the present investigation to different flaring activities of PEMEX.

Tula region
The city of Tula is located in the Mezquital Valley, in southwest Hidalgo with a total population of nearly 94 000 inhabitants and more than 140 industries.The region is semi-arid with average temperatures of 17 • C and precipitation ranging from 432 to 647 mm, increasing from north to south.In this region, the Tula Industrial Complex (TIC) is settled in an area of 400 km 2 .The major industries of the city are located within this region, including the Miguel Hidalgo Refinery (MHR), the Francisco Perez Rios power plant (FPRPP), several cement plants and limestone quarries.Other minor industries include metal manufacturing, food processing, chemical synthesis, and incineration of industrial waste.The emission of pollutants from combustion processes of these industries impacts the regional air quality.In addition, the inflow of untreated sewage water from Mexico City promotes severe pollution problems to soil and water resources (Cifuentes et al., 1994;Vazquez-Alarcón et al., 2001).According to current environmental regulations, this region is classified as a critical area due to its high emissions of SO 2 and particulate matter (SEMARNAT-INE, 2006).
The Mexican Ministry of the Environment and Natural Resources (SEMARNAT) brought the attention to the Tula area for two main reasons: (1) its relative closeness to Mexico City (60 km NW); and (2) environmental authorities of the State of Mexico sued the two major industries in the region at the beginning of 2006.They claimed that the high levels of pollutants measured by the Automatic Atmospheric Monitoring Network (RAMA) stations resulted from their emissions.In this respect, the RAMA has reported high concentrations of SO 2 in the northern area of Mexico City that can exceed the national air quality standard (0.13 ppm, 24-h average).According to the latest emission inventory in this region, the emissions rates are about 323 000 tons yr −1 of SO 2 ; 109 321 tons yr −1 of CO; 44 265 tons yr −1 of NO x ; 22 671 tons yr −1 of Particulate Matter (PM) and 44 632 tons yr −1 of VOCs (SEMARNAT-INE, 2006).
This situation has motivated several research studies evaluating the influence of Tula emissions on regional air quality, particularly in Mexico City (Querol et al., 2008;Rivera et al., 2009;de Foy et al., 2006).For instance, de Foy et al. (2009a) analyzed the influence of the SO 2 plumes of both the TIC and Popocatepetl volcano on the air quality of the Mexico City Metropolitan Area (MCMA).They used SO 2 values reported in Rivera et al. (2009) together with satellite retrievals of vertical SO 2 columns from the Ozone Monitoring Instrument (OMI).They found that during MILAGRO the volcano contributed about one tenth of the SO 2 emissions in the MCMA whereas the rest is roughly half local and half from TIC.In addition, Rivera et al. (2009) also simulated the plume from TIC with forward particle trajectories considering both the MHR and the FPRPP.They used average emissions of SO 2 obtained with a mobile mini-DOAS system.They represented each one of them with a single stack and focused on 26 March and 4 April of 2006.The model correctly captured the plume transport.The agreement of model study with column measurements supported the estimates of SO 2 and NO 2 for the TIC plume.However, all these studies do not consider explicitly the plume of elevated flares, thus the apportionment of one of the main emitters in the TIC to regional air pollution remains unknown.
The emissions from elevated flares are difficult to quantify due to the intrinsic features of their operation; such as high stacks, strong heat radiation and the shifting of flame position by wind speed changes, among other parameters (Gogolek et al., 2010).Strosher (1996) characterized the flare emissions from two medium-height oil field battery flare systems in Alberta (Canada) for both hydrocarbons and sulfur compounds.A sampling system was designed and mounted for this purpose.He reports, among others, that benzene and fluorene were abundant compounds in the field flare; and that tiophenes were present for sour solution gas flares.Hydrogen sulfide content in the sour gas ranged from 10 % to 30 %.Although he reports concentrations of combustion by-products, there were no estimates of soot emission rates in his work.
For regulatory purposes, most estimates are based on a set of emission factors published by the United States Environmental Protection Agency (EPA).However, these rely on enclosed flares burning landfill gases or on plume opacity estimates by trained personnel (Johnson et al., 2011).Technological differences like flare designs, applications and operating conditions can introduce more uncertainty to the estimates by these factors (McEwen and Johnson, 2011).In addition, these methodologies do not handle intermediate species formed in the combustion process that are either difficult to measure or to estimate, and which contribute to air pollution, such as volatile organic compounds, sulfur compounds and nitrogen compounds principally.
Other approaches to estimate emission rates include the application of stoichiometry calculations for flares in Nigeria (Sonibare and Akeredolu, 2004;Abdulkareem et al., 2009); the development of quantification methods based on novel measurement techniques for flares in Houston and Canada (Mellqvist et al., 2010;Johnson et al., 2011), and numerical studies with Computational Fluid Dynamics (CFD) techniques to lab-scale flares (Castiñeira and Edgar, 2008;Lawal et al., 2010), and field scale flares (Galant et al., 1984).Each of these methods provides estimates of emissions by considering in their methodology the physical and chemical aspects of a flare, thus yielding more representative rates.
The present work is related to the estimation of emission rates of combustion by-products released by a sour gas flare, in order to obtain information regarding the regional impact of flaring emissions from Tula Refinery, and to further extend this methodology for nationwide flaring activities.The main objectives are: (a) to estimate mass flow rates of primary pollutants emitted by a sour gas flare with a CFD combustion model which handles the contribution of soot radiation; (b) to model the plume of flaring operations; (c) to discuss possible contribution of the plume to MILAGRO supersites; and (d) to provide complementary information for related studies in the campaign regarding Tula emissions of SO 2 and soot.
The paper is organized as follows: Sect. 2 presents the modeling procedure; Sect. 3 results and discussion, and finally Sect. 4 summarizes the conclusions.

Combustion model
The CFD toolbox OpenFOAM® (Open Field Operation And Manipulation) version 1.5 is used to simulate the combustion part of this work.It is an open source finite volume library designed for continuum mechanics applications (Kassem et al., 2011;Marzouk and Huckaby, 2010;Weller et al., 1998 (Nordin, 2001;D'Errico et al., 2007;Marzouk and Huckaby, 2010).The following transport equation, which handles the turbulence-chemistry interaction, together with the equations of mass, momentum and energy is solved for each chemical species: In Eq. (1) the bar denotes a time-averaged value whereas the tilde represents a Favre-averaged quantity (Bilger, 1975).ỸJ is the mass fraction of the jth species, ρ is the density of the mixture and µ eff is the effective dynamic viscosity.The PaSR model splits the computational cell into a reacting and a no reacting zones.With this approach, mass exchange with the reaction zone drives the change in composition.The chemical source term, RR j , is scaled by κ, the reactive volume fraction which takes a value from 0 to 1.It is calculated as: In this expression, τ res is the residence time; τ c is the chemical reaction time and τ mix is the mixing time which depends on µ eff and the turbulent dissipation rate.

reactingFoam solver set-up
The flame of the flare is simulated as an unconfined nonpremixed combustion under crosswind conditions.In this work, the term crosswind is used in the context of open flames subjected to the influence of a cross flow of atmospheric air.In this case, the fuel stream exits at right angle with respect to the crossflow.The computational domain is a 2-D structured hexahedral mesh with dimensions of 1600 m wide and 500 m high.The flare is located at the center of the domain.The mesh is finer near the flare in order to refine the flame region.Imposed profiles of pressure, temperature and velocity are used as initial conditions.The velocity varied from 0 m s −1 at ground to 5 m s The data composition of the flared gas stream is scarce; however, in the model it was set to an available composition based on information provided by PEMEX.Mass fractions of 0.7, 0.2, and 0.1 were assigned to methane, hydrogen sulfide, and nitrogen respectively.Methane is considered to represent natural gas.The total mass flow rate for the model was set to a value reported by the emission inventory of IMP (IMPei).

V. H. Almanza et al.: Soot and SO 2 contribution to the supersites in the MILAGRO campaign
The Glassman chemical mechanism is used to solve the combustion chemistry, since it considers reactions of postcombustion gases (Glassmann, 2008).It considers 83 species and 512 elementary reactions organized in 11 mechanisms, which includes the formation of both sulfur and nitrogen oxides.
Soot was calculated with the phenomenological model of Moss (Moss et al., 1995) as implemented by D' Errico et al. (2007).It includes simplified terms describing the nucleation, coagulation, surface growth and oxidation (R ox ) as they apply to the balance between transport and production of soot volume fraction, f v , and particle number density, n, in non-premixed flames.
where ρ s denotes soot particulate density and N o is the Avogadro's number.The reader is referred to the references for further details regarding the expressions of the reaction rates, parameters and derivation.In this work, acetylene is taken as the nucleation species, since the nucleation rate of the model is assumed to be in proportion to the local concentration of acetylene (Brookes and Moss, 1999).Flame radiation is calculated with the P1 model (Sazhin et al., 1996;Morvan et al., 1998).The In-Situ Adaptive Tabulation (ISAT) (Pope, 1997), as implemented by D' Errico et al. (2007), is used to speed up the computing time.The mass flow of the pollutants was estimated with a slice method (see Sect. 3).The criterion for placing the slice was to consider low values of the OH radical in order to represent far-field conditions.
It is important to mention that the combustion model is being studied with more detail in a separate work (Almanza and Sosa, 2012).In that paper, the influence of wind field, gas composition and fuel mass flow on the emission rate of combustion by-products is being investigated.The Gas Research Institute chemical mechanism of natural gas, GRI3.0, that includes NO x formation is being considered for the chemistry of hydrocarbons.In addition, the chemical mechanisms of Lutz (Lutz et. al, 1988) and Leung (Leung et al, 1991) were tested for C1-C3 hydrocarbon combustion; however, better performance and numerical stability were obtained with the Glassman mechanism.For this reason, we use the Glassman mechanism in the present work (Glassman, 2008).

Air quality model
WRF-Chem version 3.2.1 is used for the air quality simulation.It is an online chemistry model fully coupled to the Weather Research and Forecasting (WRF) model (Ska-marock, 2005) developed at the National Oceanic and Atmospheric Administration (NOAA) (Grell et al., 2005).Aside from the gas-phase chemistry mechanisms, it includes several aerosol modules and photolysis schemes.

WRF-Chem set-up
A 5-day simulation period, from 00:00 UTC 22 March to 00:00 UTC 27 March of 2006 was conducted using three domains in one-way nesting configuration.It started at 00:00 UTC on 20 March with two days of spin-up.The grid cell sizes for the domains are 27, 9 and 3 km each with 100 x 100 nodes with 35 vertical levels (Fig. 1a).The parameterizations used in this work include the Purdue Lin microphysics (Lin et al., 1983;Chen and Sun, 2002), the NOAH Land Surface Model (Chen and Dudhia, 2001), the Yonsei University (YSU) scheme for the Planetary Boundary Layer (PBL) (Hong et al., 2006), the Monin-Obukov model for the surface layer (Skamarock et al., 2005), the RRTM and Dudhia schemes for the longwave and shortwave radiation respectively (Mlawer et al., 1997;Dudhia, 1989).The gravity wave drag option is used for the first domain only.Six hourly Global Final Analysis data with 1 • resolution were used for initial and lateral boundary conditions.
Four-dimensional data assimilation (FDDA) (Stauffer and Seaman, 1990) was used to nudge meteorology during the first 24 h of the simulation period.Only analysis nudging was applied for the first two domains.After sensitivity tests (not shown) we found it is better to turn off nudging in the PBL and in the lower ten levels as well, for wind, temperature and water vapor mixing ratio.
The model was run in concurrent mode for the first two domains with the chemistry turned off in order to save computing time.In this step, the Grell-Devenyi scheme was used for the convective parameterization (Grell and Devenyi, 2002).The output of the second domain was used for boundary and initial conditions of the third domain in an hourly basis.The chemistry of the plume was solved in the third domain.
WRF-Chem was run with the Regional Acid Deposition Model version 2 (RADM2) chemical mechanism (Stockwell et al., 1990;Chang et al., 1989) for the gas phase and the Modal Aerosol Dynamics for Europe coupled with the Secondary Organic Aerosol Model (MADE/SORGAM) mechanism for the aerosol phase (Ackermann et al., 1998;Schell et al., 2001), together with the FAST-J (Wild et al., 2000) photolysis scheme.Cumulus parameterization was turned off following the work of Zhang et al. (2009b), Doran et al. (2008), andLi et al. (2010).The chemical boundary conditions are set with the ideal profile provided by WRF-Chem.It consists of idealized, northern hemispheric, mid-latitude, clean environmental profiles of trace gases based upon the results from the NOAA Aeronomy Lab Regional Oxidant Model (NALROM) (Grell et al., 2005;Tuccella et al., 2012).Zhang and Dubey (2009a), report small sensitivity of forecast concentrations using the default profiles.Plume rise of the flaring plume was calculated as suggested in Beychok (1995), which considers the flame length and assumes a tilted flame by 45 • .
Since this work focuses only on flaring emissions, the national emission inventory was not included; therefore no other emission was set for WRF-Chem.

Sour gas flare
Miguel Hidalgo refinery has three elevated sour gas flares.In this work, a single flare equivalent to the three real flares is used in the combustion simulation, with the aim to represent, as best as possible, the total feed stream that is flared at Tula Refinery.This allows to simplify the assumption of gas composition and to save computing time.The total mass flow used in the combustion simulation is based on information provided by IMPei which reports 4.65 kg s −1 .However, it was decided to slightly increase the total mass flow in order to represent flow variation at the refinery, without losing generality with respect to IMPei information.For this reason, the total mass flow rate for the combustion simulation was rounded to 5 kg s −1 .Finally, the emission rates obtained with the combustion model are scaled to each of the three flares in order to compare with IMPei.
The 2-D sour gas flare simulation results are presented in Fig. 2 and Table 1.The table shows mass flow estimates of different species at two slices.These slices are located at different heights with respect to the tip of the simulated flame.S1 corresponds to the region near the flame and S2 to the region far from the flame.The slices are included mainly for two reasons (1) to track how the chemical species and physical properties of the combustion process evolve along the flame and the plume, and (2) to graphically depict the method for the mass flow calculation.Buoyancy and air entrainment vortices are relevant for open flames, particularly in the near field.Thus, the flow field is different outside and inside the plume.As a result, wind vectors are not aligned with the mean flow field of the plume, as shown in the velocity profile of Fig. SM3 of the Supplement.For this reason, there is inflow when the direction of the velocity vector with respect to the slice is greater than 180 degrees and less than 360 degrees, and outflow when the direction of the velocity vector with respect to the slice is greater than 0 and less than 180 degrees.The angle is calculated for each point along the slice taking a Cartesian coordinate system as reference.Thus, we refer to a quadrant in this context, where outflow corresponds to quadrant I and II, and inflow to quadrant III and IV.The mass flow was computed with 500 points of the mixture density, velocity and species mass fraction average fields.These points were sampled for each slice at intervals of 0.1 s, for the last 6 seconds of simulation time.Data within 2 standard deviations of a variable's profile are used in order to account for plume spread.Total flow through the slice is obtained as the difference of outflow minus inflow.Finally, the product of mixture density, velocity and species mass fraction is obtained.The result is further numerically integrated.
As expected, results indicate that the variation of the mass flow rate depends on the location of the slice.The interest in the region far from the flame is to represent as far as possible the conditions of measurements by Rivera et al. (2009), which considered the plume downwind of the point sources.These conditions can be achieved when soot oxidation by combustion is not significant; i.e OH radical concentration is similar to background levels.Slices at 30 and 40 meters above the flame had OH levels slightly higher than soot levels (not shown).At 60 m above the flame, it was considered that the influence of the atmospheric chemistry could start to be important, since the OH levels have decreased about two orders of magnitude.In addition, at this distance NO 2 levels are higher than those of NO, since NO 2 concentration increases at lower temperatures (see Fig. SM6 of the Supplement).However, the levels of both SO 2 and CO 2 are more variable suggesting an increase in the domain height and hence a longer simulation time.In addition, the combustion chemical mechanism is more accurate at high temperatures, so that inherent inaccuracies are present in the far flame region estimates, mainly related with the cooling of species as they are transported.

Emission rates of primary species
The estimates of flaring emissions mass flow rates obtained with the combustion model are presented in Table 2.It shows the rates at the two slices for the following species: SO 2 , NO x , VOCs, soot, CO 2 and CO.The table lists information of IMPei in the TIC (IMP 2006a;IMP 2006b), as well as the measurements by Rivera et al. (2009).The IMPei emission rates for SO 2 , CO 2 and CH 4 were calculated with a mass balance approach based on Asociación Regional de Empresas Petroleras de Latinoamérica y el Caribe (ARPEL) methodology (IMP, 2006b).A combustion efficiency of 98 % was considered for CH 4 according to the recommendations of British Petroleum.The H 2 S to SO 2 conversion rate was set to 100 % and to 99.5 % for CO 2 according to recommendations of the Intergovernmental Panel on Climate Change (IPCC).IMPei particle emission estimations are based on ARPEL methodology.
Total SO 2 mass flow rate obtained with the combustion model is higher than the IMPei estimate.It ranges from 2.64 kg s −1 to 3.86 kg s −1 according to the slice location.This represents a difference of 0.7 kg s −1 and 1.92 kg s −1 with respect to IMPei.On the other hand, the estimates are comparable with the RdF reported value.The rates of this work are in agreement with the confidence interval of the measurements; however, the rate at S2 is relatively high.It is important to note that aside from the elevated flares, Tula Refinery has other important sources, particularly heaters, ground flares, boilers and oxidizers.Based on previous studies conducted by IMP, elevated flares represent roughly 60 % of total SO 2 emissions from the Refinery, and about a quarter of the total emissions of the TIC.This implies that in the reported impacts by de Foy et al. (2009), flaring emissions could be comparable to the Popocatepetl volcano emissions.With respect to Table 2 reported in Rivera et al. (2009), an emission rate of about 2.82 kg s −1 can be assigned to elevated flares when taking the 1999 emission inventory of SEMARNAT-INE, and of 2.57 kg s −1 when taking the emission inventory of PEMEX.However, Fast et al. (2009) report in Table 3 of their paper a slightly higher emission rate for the TIC which would correspond to about 3 kg s −1 .The results of the combustion model of this work are comparable to these rates.For this reason, the combustion model mass flow rates can be considered representative for Tula Refinery.Rivera et al. (2009) also report high variability on the fluxes when acquiring measurements, and associate them to emission sources variation in the Tula region.In the case of flares, the contribution of wind to the variability of emissions is important, since it can diminish the combustion efficiency of flares, resulting in an increase of pollutant emissions.Perhaps this contributed to the main peak reported by these authors on 26 March, which was about 12.47 kg s −1 for SO 2 at TIC.In particular wind direction variations could have promoted flame shifts that led to the increment of emissions.This highlights the relevance of including flame-wind interaction when estimating emissions from flaring sources.
The results also suggest to increase simulation time and to refine the mesh downwind of the combustion plume.Coarse regions on a grid can promote artificial dispersion of the plume, so that small-scale structures are better resolved as the mesh is refined.For instance, turbulent eddies smaller than the plume radius enhance turbulent diffusion (Zhang and Ghoniem, 1993).In this work, artificial spreading can occur downwind, since the mesh refined region is upstream of the flame.This can promote that eddies generated after the ignition of the flame, which present relatively high concentration of SO 2 , can be slowly dissipated by the crosswind flow.As a result, species profiles along the slices are wider.Besides, the composition and velocity of the feed stream sent to the flare are important, so that our assumptions also contribute to the uncertainty of the estimates.This implies that hydrogen sulfide concentration should be lower.In addition, the 2-D domain lacks the dynamics and flame width that can be obtained with a 3-D setup.Work is in progress to account for this in both transient and steady state.
Regarding NO x , the estimate of this work is higher than IMPei but in agreement with RdF with respect to NO 2 .The reported NO 2 emission rate for the main peak of 26 March is about 0.614 kg s −1 and the upper limit of RdF is 0.36 kg s −1 .
Since IMPei reports only NO x estimates for just one flare, the comparison is made with RdF estimate.The total emission rate ranges from 0.018 kg s −1 to 0.07 kg s −1 according to the slice, both lower than the reported peak value and in agreement with the limits of RdF estimate.The influence of the combustion model chemical mechanism is important.Currently, it only includes a sub-mechanism for thermal NO formation.As methane is considered to represent natural gas, fuel NO formation is excluded to contribute to the emissions.However, it is possible that prompt NO formation could contribute in fuel-rich conditions.In addition, the presence of sulfur in the fuel stream affects flame dynamics because it can influence fuel oxidation and thermal NO formation (Alzueta et al., 2001).
As for VOCs, we obtained a lower value than the IMPei estimate.Since the chemical mechanism only considers hydrocarbons up to C3, it is not possible at this stage to take into account higher hydrocarbons present in the real stream, like C4, C5 or C6, because the computing time could be prohibitive.
Soot mass flow is lower than IMPei estimate.It ranges from 2.24 × 10 −4 kg s −1 to 2.46 × 10 −5 kg s −1 according to the slice location.IMP information suggests that elevated flares represent about 6.5 % of total emissions from the Refinery and roughly 1.45 % of total emissions coming from TIC.These fractions are based on the assumption of high soot content in the PM 2.5 estimate.Taking again Table 3 of Fast et al. (2009), results in an emission rate of roughly 5.65 g s −1 , 1.6 times greater than the total rate of IMPei, 3.37 g s −1 (Table 2).
Considering only methane for the fuel stream can underestimate the soot generated in the flame, since this hydrocarbon has the tendency to produce low levels of soot (Richter and Howard, 2000;Woolley et al., 2009).The Nagle-Strickland-Constable model for soot oxidation was used for R ox in Eq. (3).It assumes oxidation only by O radical; however, OH radical can be relevant in the oxidation step (Puri et al., 1994).Moreover, the influence of sulfur in the flame chemistry can be important in the oxidation step, since SO 3 can   mann, 2008).This could influence the location of the slice in far field conditions.Thus, it is possible that the combined presence of nitrogen and sulfur affect soot formation.In addition, the original values of soot model parameters are used in this work.Moss et al. (1995) report that model parameters are sensitive to the kind of fuel in diffusion flames.In this work, perhaps the parameters are sensitive to the presence of sulfur.This is being studied in the corresponding paper (Almanza and Sosa, 2012).Johnson et al. (2011) developed the Sky-LOSA technique, an in-situ method to quantify mass emission rates of soot from flares.They applied the method on a large-scale flare at a gas plant in Uzbekistan and determined an average rate of 2 g s −1 with an uncertainty of 33 %.According to Fig. 1 of their paper, the flare is visibly sooty.Moreover, their estimate is representative for the region near the flame.In this work, the corresponding soot emission rate is 0.22 g s −1 , roughly 1 g s −1 less than the lower limit of Sky-LOSA estimate.However, the aforementioned flare at Uzbekistan has a diameter of 1.05 m, about two times the diameter of MHR flares.
CO emission rate is about 2 times higher than IMPei at S1 slice, and about 4 times higher at S2 slice.However, IMPei estimate is only available for F3.The high variability in the stream composition of each flare as reported by IMPei can be attributed to the Refinery configuration, since the main difference between the flares relies on the process to which they are related within the facilities.This implies different fuel streams, different compositions and variable flow rates.Based on the IMPei emission rates of SO 2 , soot and CO 2 for each flare, it can be inferred that there exists a different amount in mass of carbon and hydrogen sulfide in the feed stream sent to each flare, with the possibility of including hydrogen.Unfortunately, we have no information with this respect to confirm the amount of species other than carbon and sulfur.However, the estimates of this work are comparable to measurements reported by Rivera et al. (2009), and IMP emissions inventory.
It is important to mention that the applicability of the combustion model could be extended to estimate other important species for atmospheric chemistry, in particular nitrous acid (HONO) and formaldehyde, even though they are not related with the purpose of this work.Previous research reports the importance of HONO in the morning photochemistry of the MCMA (Li et al., 2010), and current research have reported measurements of primary emissions of formaldehyde at Mont Belvieu, in the Houston Galveston Bay facilities (Parrish et al. 2011).

Air quality model performance
The performance of the meteorological fields obtained with WRF-chem for the third domain is assessed by means of the Mean Absolute Error (MAE), the Root Mean Squared Error (RMSE), the mean bias (BIAS), and the Index of Agreement (IOA) (Willmott et al., 1985;Willmot and Matsuura, 2005).The model surface variables considered for this purpose are the temperature at 2 m above ground level (T), wind speed (WS) at 10 m above ground level and wind direction (WD).With respect to the wind field, the RMSE of the vector wind difference (RMSEvec) is calculated.This statistic considers both speed and direction errors (Fast, 1995).
Surface temperature presents MAE and RMSE slightly greater than 2 • C for the considered stations, except in MER, PLA and XAL with almost 3 • C. The IOA is at least 0.9 in all stations, except in XAL where it is 0.88.A cold bias is present in most of the stations with XAL being the coldest.Wind speed presents MAE and RMSE close to 2 m s −1 in all stations, except in TPN where it is 2.20 m s −1 and 2.72 m s −1 respectively.The IOA for TPN, VIF and XAL is 0.53, 0.69 and 0.70 respectively, which suggests that the model is not capturing the dynamics with enough accuracy, particularly at TPN where the wind speed bias was the highest.RMSEvec ranges from 2 to 3 m s −1 in all stations except CHA, TPN and TAH.TPN presented the highest value with 3.93 m s −1 and a relatively high wind direction bias.Although the IOA for wind direction is rather low in some of the stations, especially at TAH with 0.54, it is greater than 0.7 and less than 0.86 in most of them.VIF presented the highest wind direction bias.
The values of RMSE for temperature and wind speed are comparable with those reported by Zhang and Dubey (2009a).In addition, the highest RMSE these authors report corresponds for wind speed at TPN with a value of 2.68 m s −1 .According to these authors the correlation coefficient for wind direction at TPN was higher than at TAH.In this work a similar behavior of the IOA for the same stations is obtained.Nevertheless, they modeled the whole MILA-GRO period and this work focuses only after 22 March as mentioned in Sect.2.2.1.Fast et al. (2009) report reasonable predictions of the wind fields at CHA station for the period from 6 to 30 March.They used FDDA with observation nudging.They reported inaccuracies at XAL and VIF because the model tended to over-predict the extent of the gap flow.Moreover, the downslope flows in their simulation did not propagated at night in stations near the western side of the basin rim, except in CUA.Similar behavior is obtained in this study, where the IOA for wind speed and direction is relatively high for CHA, and slightly lower at CUA station.In addition, wind direction bias is low for these stations.
On the other hand, the model captured relatively well the dynamics of temperature at the three supersites, with slightly better performance at T0.The MAE is less than 2 • C and the IOA is greater than 0.9 in all the sites.A warm bias is present at T0 and T1, whereas T2 had a cold bias.The model better reproduced wind speed behavior at T0, with both MAE and RMSE about 1 m s −1 , and with an IOA of 0.82.For the other sites the performance decreased; however, at T2 an IOA of 0.76 and a MAE less than 2 m s −1 were obtained.With respect to the IOA, the model better resolved wind direction at T2.This value is similar to T1 with 0.81.At T0 the model apparently had the lowest performance; however, the wind direction bias is low.In addition, the RMSEvec is 2.61 m s −1 , lower than T1 and T2.
Aside from the inherent variability in the parameterizations employed by WRF, possible reasons for the inaccuracy of model performance can be attributed to the vertical resolution of the grid, since the first full level in this work is located about 50 m above surface.In this respect, Zhang et al. (2009b) seem to place the first level below 50 m since their higher resolution within the boundary layer is around 10 to 100 m.This can give higher wind speeds with respect to surface observations since the height of the monitoring equipment is even less than the first mass level of the model.Nevertheless, Li et al. (2010) also place the first model level at around 50 m without visible impact in the performance of their simulations.Moreover, the local topographic and thermal effects are not well captured (Doran et al., 2009) and can influence the modeled wind direction at the surface.In addition, de Foy et al. (2009b) employ high resolution satellite remote sensing data for the land surface model and without data assimilation to improve the model performance.It is worth to mention that the simulation period in this study covers the last two cold surge events (NORTE2 and NORTE 3) described in (Fast et al., 2007).These add more variability to the dynamics of the simulation period.Since no convective parameterization was included for the innermost domain in this work, this could have contributed to the accuracy of the model as well.de Foy et al. (2009b) turned it on, whereas Zhang and Dubey (2009a), Zhang et al. (2009b), Li et al. (2010) and Doran et al. (2008) turned it off.
Although the statistics suggest improvement in the model set-up, the meteorological fields can be considered reliable for this study given the moderately high values of IOA for the surface wind field together with the performance of surface temperature.

Flaring plume transport
This section presents the WRF-Chem simulations of the plume by flaring activities at Tula Refinery (Fig. 1b).Only the estimates of SO 2 , soot, NO x , CO, C 2 H 2 and C 2 H 4 are considered as inputs to WRF-Chem.The latter two hydrocarbons were considered since acetylene is important to soot formation and ethylene is an important by-product in methane rich combustion (Law, 2006).Since the exact chemical composition is not known at all for each of the three flares in the refinery, the result of the combustion model for each combustion by-product was assigned proportionally to each flare according to the mass fluxes reported by the IMPei.Although IMPei reports soot emission rates for the three flares, the result of the combustion model was set equal for the three flares since it is lower than the IMPei value.The same was done with acetylene and ethylene.All of the emissions are held constant in all the simulation period.Table 4 presents the mass rates considered as inputs to WRF-Chem.

SO 2
In Fig. 3 surface measurements of SO 2 at T0 and T1 are compared against model results.At T0, the model suggests a contribution after midnight on 22 March until 05:00 LT.
There was a northerly wind on the previous day that transported the plume at this site, with a strong lateral transport from the east after midnight.On 23 March northerly winds prevailed most of the day, and the plume reached T0 again in the early morning.The model reproduced the gradual increase of the observations after 02:00 and until 07:00, with a small decrease after this hour.There were missing data in this period to confirm this, but the results suggest that it is related with topographic effects of Sierra de Guadalupe together with stability conditions.Although the timing of the main peak was slightly underpredicted, nudging was impor-tant to model relatively well the phase of the observations peaks.As the plume continued to be further transported by the northwesterly component, the model also reproduced the evolution of measured SO 2 levels for the rest of the day.According to the model, the greatest contribution of the plume in SO 2 levels at T0 occurred on this day.The plume practically spent most of the day in the basin.
On 24 March, a northeasterly component prevented the plume to be transported to the basin and practically there was no contribution at T0.The model suggests that the small peak at 19:00 is originated by downslope flow that promoted the recirculation of a residual mass of the plume located at the western side of the basin.However, the measurements show two peaks along the day with higher concentration than on 23 March, from 08:00 to 21:00 LT.A peak was measured at T1 with a similar timing to the first peak at T0.It is possible that this was a contribution from Tula, since on 23 March peaks with a similar evolution were registered at both T0 and T1 (Fig. 3).Nevertheless the model did not reproduce this since the plume moved farther to the west without reaching the basin.de Foy et al. (2009a) attributed this to subtle changes in the strength of the down-valley flow from the northeast, and the up-valley flow to Tula also from the northeast, which totally change the resulting plume transport.
Recirculation was present in the early morning of 25 March until a northerly flow gradually transported the plume back to the basin.The model suggests that the plume reached T0 after midday.26 March presents a similar behavior, with the contribution mostly by recirculation before midday and with a late northerly flow in the afternoon.Recirculation was also important in early hours of 27 March.Nearby Tula, there was flow from both the southeast and southwest that transported the plume further to the north, supporting the increment of concentration by recirculation.The gradual decrease of SO 2 levels as shown by observations is relatively well reproduced by the model.
At T1 the behavior is similar in all the simulation period.The main difference is related to the influence of local sources on SO 2 levels, since the diurnal cycle is not as visible as at T0.Nevertheless, measured values are about two times greater than at T0. On 23 and 25 March the model reproduced relatively well, with respect to observations, the timing of the plume first reaching T1 and then T0.
In order to infer about the influence of the emission rate on the modeled concentration values at the supersites, a similar simulation (LE) was conducted with a total emission rate of 1.97 kg s −1 .For 23 March results suggest maximum differences with respect to the aforementioned simulation, of about 4 ppb at T0 and 27 ppb at T1.The peak for this day was also reproduced as in the previous simulation, and in general the same behavior was obtained for the time series.For the other days the maximum differences are 7 ppb at T0 and 8 ppb at T1.These values suggest that even though the estimate of SO 2 obtained with the combustion model is high, the differences in concentration levels are relatively low in the whole simulation period at both supersites.
With these two scenarios, the average contribution of flaring emissions from MHR to the total SO 2 levels was calculated for the simulation period.Results are presented in Table 5.When using the combustion model estimate the contribution is about 37 % at T0 and 39 % at T1.With the LE mass flow rate, the average contribution is 23 % at T0 and 29 % at T1.This yields a global difference between the two scenarios of about 14 % at T0 and 10 % at T1.The corresponding average concentration is of 2.22 ppb and 1.42 ppb at T0.At T1 a concentration of 2.02 ppb and 1.48 ppb is obtained for each emission rate.
The percentage contribution to the supersites was estimated following de Foy et al. (2009a), by taking the ratio of mean model to observation values for the entire simulation period of this work.Although not shown, the contribution was also estimated by taking the area under the curve for both observations and model results for the entire simulation period and comparing them.In this case, there were variations of up to 5 % with respect to taking the mean values.For simplicity we retained the contribution obtained with the mean values.
de Foy et al. (2009a) present impact fractions in the MCMA by emissions from TIC for the whole MILAGRO campaign.They are about 40 % to 57 %, according to different configurations in their model set-up.This suggests the possibility of periods when TIC emissions can impact more than local sources.According to results of this study, it is feasible that on 23 March, emissions from MHR contributed more than local sources to the total SO 2 levels at T0.However, it is important to note that we did not include the urban emissions in this study.

Soot transport
Elemental carbon (EC) of WRF-Chem model is considered to represent soot.Model results are compared with surface measurements of elemental carbon obtained during the MI-LAGRO campaign (Molina et al., 2010).Even though meteorological fields are similar for both species, the time series between these species are different.Model concentrations of EC at T1 are compared with observations in Fig. 4. Important to note is the difference of about three orders of magni- tude with respect to observations.Possible reasons for such a difference can be attributed to the influence of local sources together with the inherent uncertainty of the estimate.
Another simulation was conducted using estimates of PM 2.5 obtained by the IMP (IMP, 2006b).The emission rate considered all the possible combustion sources for MHR and FPRPP.The purpose to include other sources of soot rather than just elevated flares is to determine if dilution accounts for the low modeled values of soot.The mass flow rate was set to 0.23 kg s −1 .Similar to de Foy et al. (2009a), one stack for FPRPP was considered and used the three stacks for the flares to represent all the MHR emissions.Results are presented in Fig. 5.
The diurnal variation at T0 and T1 is not reproduced, mainly because the national emission inventory is not included.Instead, model time series reveals the days with most possible contribution from TIC in EC levels at the three supersites.For instance, on 22 March there was small contribution at the three supersites, principally after midnight.In contrast, after 24 March most of the contribution possibly occurred late in the morning and early in the afternoon.Fast et al. (2009) shows biomass burning sources close to T0.However, after 23 March the frequency and intensity of fires diminished considerably with the NORTE3 (Fast et al., 2007).
At T1 a similar contribution of local sources is also present.In this case, the emissions from the highway connecting Mexico City and Pachuca can be the most relevant (Fast et al., 2007).Because T2 is a remote site, daily variations are more frequent so that diurnal patterns are not as pronounced as in the other sites, and dilute plumes from nearby sources are more important (Fast et al., 2009).This implies that a contribution of the plume from TIC could present a similar timing with the peak of EC observations at this site.Consequently, the influence of local sources at T0 and T1 can induce a difference in the timings between the plume of TIC and observations.For instance, on 22 March the main peak of observations at T0 occurs at 06:00 LT, while the model peak is around midnight.
This influence of local sources is clearer on 23 March.On this day the model suggests a transport from T2 to T1 to T0.At T2 the observations show that after 04:00, a peak started to develop and ended later at 13:00.The model reproduced this peak relatively well, including its gradual diminishing, and the timing of the maximum value at 05:00 LT.However, the concentration was overpredicted and the peak ended earlier, at about 11:00.As the plume continued to be transported, it reached T1 at 09:00.In contrast, the observed peak was at 06:00.Model surface fields suggest southerly wind at T1 in the early morning, with a gradual flow developing from the north later in the morning.Therefore, it is possible that the peak of the observations can be related to local sources, and that TIC emissions were more important before midday at T1.The plume reached T0 at 10:00, with the peak at 11:00.It roughly coincided with the minimum of the observations, suggesting that around midday of 23 March part of the EC levels at T0 came from TIC.The gradual decrease in EC concentration measurements can be attributed to northerly flow.In a similar way, results show that the plume impacted the three supersites on 25 and 26 March.On 26 March, most of the model EC levels at T2 are due to recirculation.In contrast, the plume directly impacts T1 first and later T0.The observations minimum value at T0 is close to the model maximum, like on 23 March.The absence of sharp peaks in the morning can be related to calendar day, since it corresponds to Sunday.On 27 March, there was a slight contribution at T1 in the morning; at T2 after midday, and no visible contribution at T0 before 17:00.Fast et al. (2009), report considerable underestimation of EC at T0 and T1 in the period from 05:00 to 10:00.In this work, the model suggests contribution by TIC at T1 in the period after 05:00 until 10:00 inclusive, on 23, 25, 26, and 27 March.Nevertheless, it was clearer on 23 March.However, since at this stage the model is not considering aqueous reactions, possible scavenging by precipitation on rainy days as suggested by Doran et al. (2008), can be relevant.
These results showed that even though the plume can dilute as it is transported, EC levels are comparable in terms of order of magnitude when taking into account all the combustion sources at TIC.For this reason, if we take the IMP estimate for flares and compare it with the total of TIC, it is about 68 times lower.This suggests that the low values of EC obtained with WRF-Chem when taking the combustion model estimate are feasible, and thus the contribution of local sources is rather more significant (Fig. 6).
Contribution of soot was estimated at the three supersites and results are presented in Table 5.When considering the total emissions of the power plant and the refinery, more than half of the total levels at T2, and less than 10 % at T0, can be attributed to TIC in the simulation period.At T1, the contribution is roughly double of T0 in all scenarios.However, when taking the estimate of the combustion model it is less than 0.1 %.
With this information a plot was constructed in order to depict the surface impact of the TIC flaring plume.It was obtained by tracking the points at which the plume exceeded a threshold value within a region encompassing Tula and MCMA.Basically it shows how much time the plume spent in this region during all the simulation period.First, a threshold value is set in this region based on the detection limits of measuring instruments.A value of 0.01ug m −3 is used for soot whilst 1 ppb is set for SO 2 .Then the number of hours in which this threshold was exceeded in all the simulation period was counted.The plot is shown in Fig. 7. Please refer to Fig. SM1 in the Supplement for the procedure.
The figure shows the time, as percent of hours for the simulation period, that the flaring plume spent in representative locations within this region.It considers total TIC soot emissions and MHR SO 2 emissions.The spatial dis-tribution is similar to that reported by Zambrano García et al. ( 2009), with a tendency of transport towards the northeast.For these primary pollutants the distribution is similar, but it can change for secondary pollutants.The plot suggests that the plume spent more time at T0 than at the other supersites.This implies that the north region of the basin was the most exposed to emissions from flaring operations.
Since the flares operate continuously, there exists the potential of a constant exposure, even though the concentration is small.These contribution plots can give further information when considering the emission inventory, and may provide supporting information for exposure in specific regions which can include vegetation, crops, soils and population.

Conclusions
This work presents simulations of the plume emitted by the three elevated flares of Miguel Hidalgo Refinery at Tula, Mexico, in order to study the contribution of flaring emissions at the MILAGRO supersites.This was accomplished in two steps.First, mass flow rates of combustion by-products were estimated with OpenFOAM® based on an equivalent sour gas flare in a 2-D configuration.This model considered the crosswind interaction with the flame and the content of hydrogen sulfide in the flared stream.The considered byproducts were C 2 H 2 , C 2 H 4 , NO x , SO 2 , CO and soot.The emission rates were calculated with a slice method.Second, the atmospheric evolution of these combustion by-products www.atmos-chem-phys.net/12/10583/2012/was simulated with WRF-Chem, for the period from 22 to 27 March.The regional simulations focused on SO 2 and soot.
The soot emission rate estimates ranged from 0.024 g s −1 to 0.2 g s −1 , compared to 2.63 g s −1 of the emission factor by the IMP.As for SO 2 , the mass flow estimates at slices S1 and S2 ranged from 2.64 kg s −1 to 3.86 kg s −1 , which are in the range of measurements obtained by Rivera et al. (2009), 4.90 ± 3.80 kg s −1 .In addition, the calculated rate for NO 2 is 7.07 × 10 −2 kg s −1 at S2 and measurements by Rivera et al. (2009) suggest 2.77 × 10 −1 ± 8.10 × 10 −2 kg s −1 .Although the combustion model requires further improvement, particularly in selecting the location of the slice and the gas stream composition, the estimates of this work can be helpful for delimiting the combustion emissions by the flares of MHR.The applicability of the combustion model can be extended to estimate other important species relevant for atmospheric chemistry, particularly nitrous acid and formaldehyde, a highly reactive organic compound.
The impact of SO 2 and soot on the air quality at the three MILAGRO supersites was further evaluated.Given the relatively good values of MAE, RMSE, BIAS, RMSEvec and IOA, the performance of WRF-Chem was reliable enough to compare results with surface measurements at MILAGRO supersites.Analysis nudging was important to capture relatively well the timing when the plume reached the three supersites, given the relative agreement of model peaks with observation peaks, especially on 23 March.However, further improvement is recommended, particularly the refinement of the vertical levels within the PBL, and to consider convective parameterization in the inner domain.
The results from the present study suggest a more feasible contribution of TIC in SO 2 levels on 23 March in most of the basin and at the three supersites, with similar behavior on 25 March of 2006, and with a potential contribution on 24 March according to measurements.The estimated contribution of elevated flares to total SO 2 levels at MILAGRO supersites, is about 37 % at T0 and 39 % at T1.The high contribution values can be attributed to the persistence of the plume in the basin on 23 March.Results showed a transport of emissions from T2 to T1 to T0 on that day, when meteorology favored northeasterly winds due to the presence of the third cold surge.At T2 the model peaks compared relatively well with observations.
As for soot, the estimated contribution of flares to total soot levels was less than 0.1 % when taking the combustion model estimate.Nevertheless, when considering all the possible emission sources at TIC, including the power plant, the contribution is about 7 %, 17 % and 59 % at T0, T1 and T2 respectively.Concentration values of EC less than 1 µg m −3 are feasible at the supersites as well as in the basin.The model suggested that with respect to local sources, the flaring plume has less influence at the supersites, so that background concentration was higher than maximum model concentrations in almost all the simulation period, particularly on 23 March.Since the national emission inventory was not included, this could have contributed to the difference between observed and modeled concentration of soot.Even though results suggest that EC emissions from flaring at MHR are not significant for urban locations within the MCMA, they can exert a greater contribution to both urban and rural areas near the refinery.
For the simulation period of this study, the flaring plume spent more time at T0 than at the other supersites.This implies that the north region of the basin could have had higher exposure to TIC pollutants, which is in agreement with previous studies.It also showed that when considering emissions of the power plant, the plume can reach part of the south region of the basin.
These results complement previous findings of studies related to TIC by other research groups, and at the same time give the possibility to extend this work in studying the contribution of flaring activities to levels of secondary pollutants both in the MCMA and near the refinery.In addition, it is feasible to apply these tools to a country-wide analysis of the impact of flaring activities in Mexico.For instance, it can extend the modeling of air quality emissions of offshore flares in Campeche Sound previously conducted by IMP (Villaseñor et al., 2003).In addition, this information can support related studies regarding the possible recovery of the gas in the refinery.

Table 4 .Fig. 3 .
Fig. 3. Comparison of modeled (blue) and observed (red) time series of SO 2 at T0 and T1, for the simulation period.Black dashed lines are plotted every 6 h.

Fig. 4 .
Fig. 4. Time series of EC obtained with the soot mass rate obtained with the combustion model (dashed blue) and measurements (red) at T1.

Fig. 5 .
Fig. 5. Model (blue) and measurements (red) time series of EC concentration for the simulation period.Results correspond to total TIC emissions (Power plant plus Refinery).T0 (upper panel); T1 (middle panel) and T2 (bottom panel).Black dashed lines are plotted every 6 h.

Fig. 7 .
Fig. 7. Suggested spatial contribution of SO 2 and EC in the period from 22 to 27 March of 2006.

Table 1 .
Mass flow rates obtained with the combustion model.The capital S and the number at the right represent the number of the slice.The number in parenthesis shows the height of the slice with respect to the tip of the simulated flame.All the estimates are in kg s −1 .
* SO x emissions.

Table 3 .
Performance statistics of the surface variables in the third simulation domain.Location with respect to the center of the city is in front the station name.