**Contents**


Reprinted from: *Int. J. Environ. Res. Public Health* **2020**, *17*, 8349, doi:10.3390/ijerph17228349 . . . **151**

## **About the Editors**

#### **Isidro A. P´erez**

His scientific career began with photochemical oxidant measurements and modelling. The following research line was the management of RASS sodar data, which led to different papers. This device analyses the low boundary layer beyond the heights investigated by any meteorological tower. After that, CO<sup>2</sup> concentrations were studied at a rural site and several papers that combined RASS sodar observations with concentrations measured present the results of this research. Moreover, parametric and nonparametric procedures were used to investigate the daily and annual cycles of greenhouse gases. The analysis of air parcel trajectories was later considered, and the contrast between Atlantic and continental trajectories was observed on the greenhouse gas concentration recorded. Finally, his research is oriented towards the trend determination by parametric and nonparametric procedures and the origin of greenhouse gas outliers and their influence on this trend.

#### **M. Angeles Garc´ıa ´**

Her research has included photochemical oxidant measurements and modelling. Special attention was devoted to ozone analysing temporal variations, air quality and the influence of meteorological variables. Another research topic in which she has participated emphasizes the greenhouse gases included in some projects with the aim of studying the climate variability and the need to quantify the sources and sinks of CO<sup>2</sup> of natural ecosystems. Different procedures have been applied to investigate temporal cycles and trends of greenhouse gases, the influence of the air masses and the atmosphere evolution on the concentrations recorded. The research has also included the study of CO<sup>2</sup> fluxes in an agricultural soil.

## **Preface to "Air Pollution Meteorology"**

A knowledge of human influence on pollutant concentrations in the atmosphere is both a current and a major issue, with the growth of these concentrations proving key to quantifying possible risks for both the population and materials. The air pollution system consists of three parts: the first is the sources that release these substances, the second is the atmosphere, which is a substratum where they evolve, while the receptors that experience the consequences of air pollution are the third part of this system. The atmosphere is beyond human control and is where chemical reactions as well as pollutant transport and removal processes, such as deposition, occur. These processes may be significantly conditioned by meteorological variables.

Air pollution covers a wide spatial range, and certain episodes present a local scope. For instance, ozone formation is favoured under high-pressure systems, with clear skies, a hot and dry atmosphere, strong solar radiation, and weak convection. Moreover, meteorological variables, such as relative humidity, may determine the type of air pollution, such as wet haze, which is linked to high relative humidity, and dry haze, when relative humidity is low.

However, pollution transport may extend the spatial range of this problem. Although transport may help to decrease air pollution at a local level due to the inverse relationship between concentrations and wind speed, certain substances, such as particulate matter, are sometimes transported from sources like deserts to sites where they could determine pollution episodes. They may even be accompanied by toxic substances, since some sites swept by air parcel motion are war zones. In some situations, microorganisms have even been transported, and avoiding outdoor activities for at least two days after these episodes is recommended. The air parcel path conditions this transport. For instance, photochemical ozone production is amplified by polluted and hot air coming from certain directions, whereas low ozone concentrations are linked to clean, cool, and dry air from a different site. Moreover, although transport is determined by orographic features, the wind may sometimes overcome these barriers.

Various procedures are used to quantify the influence of meteorology on air pollution, such as the direct relationship between meteorological variables and the concentrations of certain substances. However, a more elaborate treatment of observations is sometimes required, such as the analysis of synoptic meteorological patterns. Chemical-transport models are complex and help to supplement procedures aimed at investigating air pollution events under varied conditions. In particular, sensitivity analyses can be performed to investigate model response to changes in input variables. Moreover, meteorological conditions that favour high concentrations may be investigated.

This link between meteorology and air pollution is usually considered to be of secondary importance, yet merits additional study in order to enhance current knowledge. The result of this research will be useful not only from a theoretical point of view, i.e., in terms of gaining insights into atmospheric processes, but also with regard to implementing actions in the field of air pollution control. Consequently, exploring what impact meteorology has on air pollution processes is a promising research field, with varied study options and whose results have a direct social application. This book is aimed at improving current knowledge of this relationship. Its content includes a range of meteorological variables and synoptic patterns to describe atmospheric evolution. Air pollution is characterised by particulate matter, harmful substances and indoor air quality. Finally, the relationship between them is established by statistical procedures such as cluster analysis or multiple linear regression, while calculations such as air parcel trajectories are also considered. The result is a book that presents current research on the subject and which is geared towards increasing and sharing experience in this field with scientists and managers responsible for taking decisions to prevent or deal with air pollution situations.

> **Isidro A. P´erez, M. Angeles Garc´ıa ´** *Editors*

International Journal of *Environmental Research and Public Health*

## *Article* **Pollution Sources and Carcinogenic Risk of PAHs in PM<sup>1</sup> Particle Fraction in an Urban Area**

**Ivana Jakovljevi´c 1,\* , Zdravka Sever Štrukil <sup>1</sup> , Ranka Godec <sup>1</sup> , Ivan Bešli´c <sup>1</sup> , Silvije Davila <sup>1</sup> , Mario Lovri´c <sup>2</sup> and Gordana Pehnec <sup>1</sup>**


Received: 30 September 2020; Accepted: 17 December 2020; Published: 21 December 2020 -

**Abstract:** Airborne particles are composed of inorganic species and organic compounds. PM<sup>1</sup> particles, with an aerodynamic diameter smaller than 1 µm, are considered to be important in the context of adverse health effects. Many compounds bound to particulate matter, such as polycyclic aromatic hydrocarbons (PAH), are suspected to be genotoxic, mutagenic, and carcinogenic. In this study, PAHs in the PM<sup>1</sup> particle fraction were measured for one year (1/1/2018–31/12/2018). The measuring station was located in the northern residential part of Zagreb, the Croatian capital, close to a street with modest traffic. Significant differences were found between PAH concentrations during cold (January–March, October–December) and warm (April–September) periods of the year. In general, the mass concentrations of PAHs characteristic for car exhausts (benzo(ghi)perylene (BghiP), indeno(1,2,3-cd)pyrene (IP), and benzo(b)fluoranthene (BbF)) were higher during the whole year than concentrations of fluoranthene (Flu) and pyrene (Pyr), which originated mostly from domestic heating and biomass burning. Combustion of diesel and gasoline from vehicles was found to be one of the main PAH sources. The incremental lifetime cancer risk (ILCR) was estimated for three age groups of populations and the results were much lower than the acceptable risk level (1 × 10−<sup>6</sup> ). However, more than ten times higher PAH concentrations in the cold part of the year, as well as associated health risk, emphasize the need for monitoring of PAHs in PM1. These data represent a valuable tool in future plans and actions to control PAH sources and to improve the quality of life of urban populations.

**Keywords:** BaP; HPLC; carcinogenic; diagnostic ratio

#### **1. Introduction**

Particulate matter (PM) is assumed to be among the most hazardous of all ambient pollutants. Particle pollution contains "inhalable coarse particles" with diameters larger than 2.5 µm and smaller than 10 µm and "fine particles" with diameters of 2.5 µm or smaller. One of the most significant organic groups bound to PM in terms of health risk are polycyclic aromatic hydrocarbons (PAH). These compounds can exist in the atmosphere in the vapor phase (PAHs with low molecular weight), whereas heavier ones (PAHs with high molecular weight) are mostly adsorbed on the particle phase [1]. PAHs generally occur as complex mixtures, products of incomplete combustion processes, and originate from natural and anthropogenic sources [2]. High PAH levels in the ambient air of large metropolitan cities are usually linked with traffic, as well as diesel and gasoline automobiles [1,3]. Catalytic converters have shown a significant effect on reducing levels of the PAH concentrations in exhaust gases, but PAH emission levels continue to increase due to the contribution of other sources, such

as traffic congestion [4]. PAHs are always emitted as a mixture, and the molecular concentration ratios are considered to be typical for a given emission source [5]. The toxicity, carcinogenicity, and mutagenicity of aromatic hydrocarbons have led to increased concerns for human populations. Long-term exposure to PAHs can cause toxic effects such as breathing problems, lung function abnormalities, decreased immune function, kidney and liver damage, skin irritation, and inflammation. The most significant health effect to be expected from inhalation of PAHs is increased risk of lung cancer primarily in occupations with heavy exposure to traffic-related air pollution, such as policeman [6] and newsagents [7]. Piccardo et al. [8] noticed that taxi drivers in Genoa were exposed to significantly higher daily BaP concentrations in comparison to workers from another occupational category. Fuel and biomass combustions, traffic emissions, and use of lubricant oils were identified as the main sources of PAH exposures at five Portuguese fire stations. During a normal shift in non-fire situations, levels of light PAHs were predominant and may promote some adverse health outcomes [9]. Bladder cancer is also linked to exposure to PAHs. Boada et al. [10] conducted a case study in the Canary Islands measuring PAH serum levels of 140 patients diagnosed for bladder cancer. Their results showed a difference in PAH contamination profile between patients and the control group, leading to the conclusion that specific PAH mixtures play a role in bladder cancer. People are persistently exposed to PAHs. Benzo(a)pyrene (BaP) is the most investigated PAH, and most information on the toxicity and manifestation of PAHs is related to this compound, which is why it was used as an indicator of carcinogenic hazard in polluted environments [11–13].

Although particle composition has been determined by many authors, there is still not much research related to the study of PAHs in airborne fine particles. Fine particles, with an aerodynamic diameter of less than 2.5 µm (PM2.5), and especially particles with aerodynamic diameter of less than 1 µm (PM1), may play an important role in affecting human health. Squizzato et al. [14] stated the following reasons: they penetrate more effectively into the deep lung; they can penetrate more readily into indoor environments; they can remain suspended for longer periods of time in the atmosphere than coarse particles; they may be transported over long distances; they tend to carry higher concentrations of more toxic compounds, including acids, heavy metals, and organic compounds; and they have a larger surface area per unit mass compared to larger particles and can thus absorb larger amounts of semi-volatile compounds. Measurements of PM<sup>1</sup> and its content have so far not been included in routine measurement programs, although they would provide more information about potential PM sources and could be used to improve PM control strategies and health protection. A recent study by Yang et al. [15] found that both PM<sup>1</sup> and PM2.5 levels were associated with poorer lung function in children, with stronger associations for PM<sup>1</sup> compared to PM2.5, pointing to the importance of regulating finer PM fractions. Furthermore, a lot of previous health risk estimations were carried out taking into account larger fractions of particulate matter (PM<sup>10</sup> or even total suspended particles), which do not enter deeply into the respiratory system [16–18].

Studies regarding PAH in PM<sup>1</sup> were carried out only at a few urban locations in Europe [19–25]. It was found that the levels varied significantly between heating and non-heating season [22,23,26]. The highest concentrations were observed during winter in areas where coal and wood burning were used for heating [22,24]. Traffic-loaded sites showed a large contribution of PAHs with larger molecular weight [21,22,27,28].

In this study, the mass concentrations of individual PAHs in the PM<sup>1</sup> particle fraction were determined in an urban area together with some gaseous pollutants. The relationship between measured species and meteorological parameters (temperature, relative humidity, atmospheric pressure, wind direction, and velocity) was determined. Potential pollution sources were assessed using PAH diagnostic ratios, Spearman's regression, and principal component analysis (PCA). In a previous study at the same location [19], factor analysis (FA) was used to identify pollution sources. The aim of both PCA and FA is to determine the main relationships between the observed variables and to reduce the large number of observed variables to a smaller number of factors. Mathematically, in the PCA, the entire variance in the observed variables is analyzed, while in the FA only the mutual variance is

analyzed. Attempts are made to identify and eliminate variance due to error and variance specific to each individual variable. The assumption of FA is that factors "create" variable values, while in PCA it is assumed that variables "create" components. In PCA, there is no basic theory that would define how variables are grouped into factors, the variables are simply empirically related within components. As atmospheric PAH levels are affected by different weather conditions and numerous atmospheric reactions, and it is not possible to determine their interdependencies for each day, we estimated that the PCA method is more applicable for this analysis. In addition, the adverse health effect of PAHs on a human organism was assessed using toxic equivalent concentrations. The incremental lifetime cancer risk (ILCR) was estimated for three age group of populations, which presents the first such study carried out for PAHs in PM<sup>1</sup> in this part of Europe.

#### **2. Materials and Methods**

#### *2.1. Location and Sampling*

Concentrations of 11 PAHs in PM<sup>1</sup> particle fraction were measured continuously from January to December 2018, together with gaseous pollutants (CO, NO2, SO2, and O3). The measuring station (45◦50′6.83" N, 15◦58′42.12" E, 168 m a.s.l.) was situated in the northern residential area of Zagreb, the Croatian capital (~800.000 inhabitants). It is a low-rise building area with small inhabitant density and mild traffic density. Residential (domestic) heating relies mostly on gas, but some households still use oil or wood for heating and cooking.

Samples of PM<sup>1</sup> particle fraction were collected on quartz filters with low-volume sequential automatic sampler (Sven Leckel) from about 55 m<sup>3</sup> of air per day. Filters were collected 24-h a day, and the total number of samples was 363. The sampler inlet was located approximately 1.5 m above ground and 15 m away from the road. The PM<sup>1</sup> samples were kept frozen in aluminum foil at −18 ◦C until PAH analysis to avoid PAH losses and sample degradation. The filters were extracted no later than two months after sampling.

#### *2.2. Measurements of Gaseous Pollutants and Meteorological Data*

Gaseous pollutants were measured using automatic devices. SO<sup>2</sup> was determined using a HORIBA APSA 370 device according to the norm (EN 1421:2014). NO<sup>2</sup> was measured by a HORIBA APNA 370 device (according the norm EN 14211:2012), O<sup>3</sup> was measured using a HORIBA APOA 370 device (according the norm EN 14625:2012), and CO was measured using a HORIBA APMA 370 device (according the norm (EN 14626:2012). Measurements of SO2, CO, NO2, and O<sup>3</sup> are part of the local air quality monitoring network funded by the City of Zagreb, City Office for Economy, Energetics and Environment Protection.

Meteorological data (temperature, relative humidity, pressure, amount of rainfall, wind direction, and velocity) were obtained from the nearest meteorological station (Maksimir) of the Croatan Hydrological and Meteorological Service (www.meteo.hr).

#### *2.3. Analysis of PAHs*

Extraction of PAHs from filters was performed with a solvent mixture (toluene:cyclohexane 7:3) in an ultrasonic bath. Further preparation included centrifugation (10 min, 3000 rpm) and evaporation to dryness. After that, samples were redissolved in acetonitrile. Concentrations of PAHs were determined by Agilent Infinity 1260 high-performance liquid chromatography (HPLC) with a fluorescence detector. Zorbax Eclipse PAH column (100 × 4.6 mm) was used for separations of PAHs. The mobile phase was acetonitrile and water (60:40), with a flow rate of 1 mL min−<sup>1</sup> [29,30]. Calibration curves were prepared with a commercial PAH standard (Supelco EPA 610 PAHs Mix). The calibration range was from 0.005 ng µL −1 to 0.08 ng µL −1 for pyrene, benzo(a)anthracene, chrysene, benzo(k)fluoranthene, benzo(a)pyrene, and indeno(1,2,3,cd)pyrene, but for fluoranthene, benzo(j)fluoranthene, benzo(b)fluoranthene, dibenzo(ah)anthracene, and benzo(ghi)perilene the

calibration range was from 0.01 ng µL −1 to 0.16 ng µL −1 . Samples were analyzed for the following PAHs: fluoranthene (Flu), pyrene (Pyr), benzo(a)anthracene (BaA), chrysene (Chry), benzo(j)fluoranthene (BjF), benzo(b)fluoranthene (BbF), benzo(k)fluoranthene (BkF), benzo(a)pyrene (BaP), dibenzo(ah)anthracene (DahA), benzo(ghi)perylene (BghiP), and indeno(1,2,3-cd)pyrene (IP). The limit of detection (LoD) was from 0.001 ng m−<sup>3</sup> for BaA to 0.03 ng m−<sup>3</sup> for BjF. The quantification limit (QL) varied from 0.002 ng m−<sup>3</sup> for BaA to 0.1 ng m<sup>3</sup> for BjF. The accuracy of the method was determined by analyzing the standard reference material (SRM 1649a, urban dust) provided by the National Institute of Standards and Technology (NIST) and ranged from 88% for Flu to 109% for BkF. Samples of SRM 1649a were processed the same way as real samples. The details of detection and quantification limits, as well as accuracy are shown in Table S1 of the Supplementary Data section.

#### *2.4. Statistical Analysis*

The statistical results were processed by Microsoft Excel and the Statistica 13 (Tibco Software Inc.) program. Statistical significance was set at 5% (*p* < 0.05). The Shapiro-Wilk test was used to test the normality of variables. Seasonal differences between warm and cold period for each PAH concentrations, heavy and light PAHs, and PM<sup>1</sup> particles were tested by Mann–Whitney U test.

The concentration ratios between selected PAHs were investigated to find out more about the nature of pollution sources. Diagnostic ratios are common tools for the identification of pollution sources [30–33]. However, some papers have presented their restrictions. Katsoyiannis et al. [34] demonstrated that diagnostic ratios cannot be effective as markers of sources, because they are influenced by weather condition and atmospheric reaction of PAHs such as photodegradation, and reaction with ozone and other atmospheric pollutants. To minimize confounding factors such as dissimilarities in volatility, water solubility, and adsorption, diagnostic ratio calculations usually are restricted to PAHs of similar molecular mass [35]. Due to the restrictions of the diagnostic ratio, the multivariate principal component analysis (PCA) was performed on data and calculated anomalies for the studied samples. Principal components (PCs) with eigenvalues greater than 1 for both periods were calculated, and their contributions in the total variance were determined. We graphically displayed the PCs loading as well as the orientation of the variables and samples with respect to these principal components.

#### *2.5. Carcinogenic Activity and Population Exposure*

Benzo(a)pyrene was used as an indicator of health impact of the PAH mixture to human health. Many studies have shown that BaP was present at more than 50% in total carcinogenic activity and that made it a good indicator of carcinogenic activity of the PAH mixture. The carcinogenic activity of individual PAHs was estimated on the basis of toxic equivalency factors (TEFs) from the literature. Different TEF schemes are developed by different authors, based on experiments in animals [36–39]. Nisbet and LaGoy [39] completed a new list of TEFs, which appears to better reflect the current state of knowledge on the relative potency of individual PAHs. In this study we used Nisbet and LaGoy's [39] toxic equivalent factors.

For calculating the risk of the PAH mixture in ambient air, the carcinogenic potencies of individual PAHs are expressed relative to the potency of BaP. BaP equivalents (BaPeq) were calculated by multiplying the mass concentration of an individual PAH with its respective toxic equivalency factor. Total carcinogenic potency (TCP) was calculated by summing up the BaPeq of each measured PAH. The equation used to calculate the TCP is presented below (1):

$$\text{TCP} = \Sigma \text{BaP}\_{\text{eq}} = \Sigma \gamma \mathbf{i} \cdot \text{TEFi} \tag{1}$$

TCP—total carcinogenic potency/ng m−<sup>3</sup> BaPeq—BaP equivalents concentration

TEFi—toxic equivalency factor of a particular PAH

γi—mass concentration of individual PAH/ng m−3.

To determine the daily population exposures, total carcinogenic potency was used to calculate the daily dose according to Equation (3). In this study we tried to estimate the most probable scenario for three age groups: infant (0–1 year), children (5–19 year), and adult (20–70 year). If we assumed that people spent an average of ten hours at their job/school and eight to ten hours at home (including sleeping), we can assume that they were elsewhere for the rest of their day, due to the fact that people spend approximately 25% of their time outdoors, that is 6 h per day.

The incremental lifetime cancer risk (ILCR) posed by exposure to PM1-bounded PAHs was computed following Equation (2) [20].

$$\text{ILCR} = \left( (\text{SF}\_{\text{inh}} \times \text{IEL} \times \text{EF} \times \text{ED}) (\text{BW} \times \text{AT} \times \text{cf}) \right) \times \text{y} \tag{2}$$

SFinh—inhalation cancer slope factor of BaP/kg day mg−<sup>1</sup> IEL—BaPeq daily dose/ng day−<sup>1</sup> EF—exposure frequency/day year−<sup>1</sup> ED—daily exposure level/µg g−<sup>1</sup> BW—body weight/kg AT—average time/day Cf—conversion factor (10−<sup>6</sup> ) y—age

In this study, SFinh, EF, and ED was used as derived by Chen and Liao [17]. Parameters were different for infants, children, and for adults, with average body weights of 6.79 kg, 36.24 kg, and 59.78 kg, respectively. Parameters that are used in calculation of ILCR are shown in Table S3 of the Supplementary Data section. Results of ILCR shows incremental cancer risk per year and for lifetime cancer risk this result should be multiplied by the age of the person.

In this study, a SFinh of 3.14 kg day mg−<sup>1</sup> was used, and the proper description of how this value was derived was explained by Chen and Liao [17]. The BaPeq daily dose was calculated by multiplying the TCP concentrations (ng m−<sup>3</sup> ), inhalation rate (IR, m<sup>3</sup> day−<sup>1</sup> ), and daily exposure time span (t, h). The equation used to calculate the IEL is presented below (3):

$$\text{IEL} = \text{TCP} \times \text{t} \times \text{IR} \tag{3}$$

#### **3. Results and Discussion**

#### *3.1. PAH Concentrations*

The monthly average mass concentrations of the PM<sup>1</sup> particle fraction during the whole calendar year are shown in Figure 1. This result shows seasonal differences (*p* < 0.001) of PM<sup>1</sup> concentrations with high values during the cold period of the year (January–March; October–December) and lower values during the warm period (April–September). The average 24-h concentrations of PM<sup>1</sup> varied from 0.7 to 55.3 µg m−<sup>3</sup> with an annual average of 13.6 µg m−<sup>3</sup> . Monthly average concentrations were also calculated for the PAH sum and for individual PAHs, and the results are presented in Figures 2 and 3. Sums of PAHs concentrations (ΣPAHs) ranged from 0.437 ng m−<sup>3</sup> during June to 21.497 ng m−<sup>3</sup> during December, and the annual mean mass concentration was 6.354 ng m−<sup>3</sup> . The highest mass concentration during the cold period was measured for BbF, while in the warm period, the highest mass concentrations were determined for BghiP. High BghiP mass concentrations in the warm period in comparison to other hydrocarbons were probably due to BghiP stability at high temperatures. The lowest mass concentrations for both periods (cold and warm) were determined for DahA. The range of the monthly average concentrations for DahA was from 0.010 ng m−<sup>3</sup> (May) to 0.386 ng m−<sup>3</sup> (December). The average monthly mass concentrations of BaP ranged from 0.038 ng m−<sup>3</sup> (June) to 2.826 ng m−<sup>3</sup> (December), while the annual mass concentrations were 0.765 ng m−<sup>3</sup> . The target value in the European Union set by Directive 2004/107/EC for BaP content in the PM<sup>10</sup> fraction is

only 1 ng m−<sup>3</sup> averaged over a calendar year. The reason for higher concentration exclusively during the cold part of the year can be explained because of the long temperature inversion periods that occur in Zagreb. Such weather conditions are characterized by high pressure and impaired air mixing, which favors the pollutants accumulation. These stable atmospheric conditions combine with increased emissions from heating are the basic origins of elevated concentrations of PAHs during the cold season. In the study area, gas, oil or, wood are mostly used for domestic heating, while coal has not been in use for more than thirty years.

The relevant literature comprises a very limited number of papers related to PAHs in the PM<sup>1</sup> fractions [19–24,27,28,40,41]. Table S2 of the Supplementary Data section shows the summary results of PAH mass concentrations obtained in this study in comparison with similar studies. Rogula-Kozłowska et al. [40], Kozielska et al. [22], and Majewski et al. [20] reported much higher mass concentrations of PAHs in the PM<sup>1</sup> particle fraction during both the heating and non-heating seasons at different locations in Poland compared with Zagreb. Much higher winter concentrations were also found in Ostrava-Radvanice (industrial site), in Kladno-Švermov and Brno (urban sites) in Czech Republic, while levels of PM<sup>1</sup> and ΣPAHs were similar in Košetice and Celakovice. In a coastal area in Poland ˇ (Baltic Sea), BaP levels were much higher than in this study, while in a coastal area in Greece PM<sup>1</sup> the levels were similar as in Zagreb but particle-bound PAHs were significantly lower. Similar PAH concentrations were also found in Guadalajara Metropolitan Area in Mexico while PAH concentrations determined in the Metropolitan Area of Porto Alegre, Brazil, were much lower during winter but higher during summer [28].

**Figure 1.** Monthly mass concentrations of the PM<sup>1</sup> particle fraction during one calendar year.

**Figure 2.** Monthly mass concentrations of ΣPAHs (polycyclic aromatic hydrocarbons) measured during one calendar year.

**Figure 3.** Monthly mass concentrations of PAHs measured during one calendar year. \* The target value (TV) in the European Union set by Directive 2004/107/EC for BaP content in PM<sup>10</sup> fraction is 1 ng m−<sup>3</sup> averaged over a calendar year.

As significant differences were found between cold and warm months, we divided the results into two groups, separately for the cold (January–March; October–December) and the warm (April–September) period of the year. Summary statistical parameters for all of the measured PAHs for these two periods are shown in Table 1. Concentrations of PAH characteristic for domestic heating or biomass burning (Flu, Pyr) were lower in the cold period than those of PAHs characteristic for car exhausts (BghiP, BbF, and IP). This indicated that the PAHs in the PM<sup>1</sup> particle fraction could originate predominantly from car exhausts [42–44].


**Table 1.** Average, minimum, and maximum values and standard deviations of the PAH 24-h mass concentrations during the cold and warm period (ng m−<sup>3</sup> ).

Cmin—minimum value; Cmax—maximum value; C—arithmetic mean; SD—standard deviation; and N—number of samples. \* Seasonal differences were tested by Man–Whitney U test (*p* < 0.001).

According to the number of rings and molecular weight, PAHs can also be classified into two groups: heavy and light. Heavy PAH concentrations were calculated as the sum of PAHs with five or more aromatic rings, and light PAHs represent the sum of PAHs with four aromatic rings. Heavy PAHs are usually characteristic for car exhausts, while light PAHs originate mostly from domestic heating or biomass burning [43,44] (Figure 4). Figure 4 shows the monthly mass concentration for these two groups. During both measurement periods, the contributions of heavy PAHs were much higher (>69%) than those of light PAHs (Table 2). Concentrations of PM1-bounded PAHs during warm periods can easily evaporate from particles to gas phase, so their concentrations in the gas phase increased, but in the particle phase their concentrations decreased. The reason behind this could be that for these PAHs the dominant sources were biomass burning and during the warm period the effect of these sources were minimal. These results also indicated that traffic (diesel or gasoline) could be the main pollution source of PAHs in PM<sup>1</sup> in this area. To confirm the assumption, we performed diagnostic ratio and principal component analysis to identify possible pollution sources.

**Figure 4.** Mass concentrations of light and heavy PAHs measured during one calendar year.

**Table 2.** Contribution of light and heavy PAHs in the sum of measured PAHs during cold period and warm period.


\* Seasonal differences between warm and cold period (test by Mann–Witney U test; *p* < 0.001).

#### *3.2. Diagnostic Ratios*

Specific ratios of individual PAHs are characteristic for some combustion processes and many authors used diagnostic ratios of PAHs to identify potential pollution sources [10–13]. In this study, the following ratios were selected: IP/(IP + BghiP), BaP/BghiP, Flu/(Flu + Pyr), and BaP/(BaP + Chry). The values determined in this study were compared with the same diagnostic ratios computed from characteristic emission sources based on the relevant literature. IP/(IP + BghiP) ratio values between 0.35 and 0.7 are characteristic for diesel, while some authors reported a value of 0.56 that indicates coal combustion [28,45]. BaP/BghiP values between 0.3 and 0.4 are characteristic for traffic, 0.46–0.81 for diesel combustion, and 0.9–6.6 for coal combustion [28,31,46]. A Flu/(Flu + Pyr) ratio between 0.2 and 0.5 indicates diesel, a ratio between 0.4 and 0.5 liquid fossil fuel, while a ratio >0.5 suggests wood combustion [32,33]. Finally, a BaP/(BaP + Chry) ratio of <0.5 indicates diesel but >0.5 gasoline [4,28]. Average PAH diagnostic ratios for the cold and warm period are presented in Table 3.


**Table 3.** Comparison of PAH average diagnostic ratios in PM<sup>1</sup> during cold and warm period and the main emission sources.

In this paper, the average IP/(IP + BghiP) ratios were 0.5 during both the cold and warm periods, suggesting that the produced PAHs stemmed from the emission of diesel vehicles. Those average values (0.5) for cold and warm period also corresponded to the 97th and 99th percentile value, respectively, i.e., during the cold period 97% of days had a ratio value lower than 0.5, while during the warm period 99% days were with a ratio of <0.5. The average BaP/(BaP + Chry) and Flu/(Flu + Pyr) ratios were 0.5 during both periods. The average ratio BaP/(BaP + Chry) in the same time represents the 60- and 67-percentile for cold and warm periods, respectively. The other ~40% of results were higher than average value, and the individual value did not exceed 0.7, indicating both diesel and gasoline combustion as potential sources. The average ratio value Flu/(Flu + Pyr) was marginal between biomass burning, emissions from gasoline or diesel vehicles, and combustion of other liquid fossil fuels, but during the warm period there were more days (15%) with a ratio value characteristic for wood burning (during the cold period 96% of results was <0.5). Because of that, it could be concluded that during the warm period mixed sources are probably present (wood combustion and emissions from gasoline or diesel vehicles). The BaP/BghiP ratio was 0.9 during the cold and 0.6 during the warm period. In the warm period, the value was characteristic for emission from diesel (Table 3), but in the cold period it was characteristic for mixed sources (coal, wood combustion, and emission from diesel). The BaP/BghiP ratio included 70% and 60% results during the cold and warm period, respectively. The other 30% during the cold period were higher than 0.9, which indicated that coal or wood combustion were present as pollution sources. During the warm period, 40% of results were higher than 0.6 but did not exceed 0.8, which is the highest value for diesel combustion. Similar results were reported by Agudelo-Castañeda and Teixeira [28] and Hanedar et al. [46]. They also found that PAHs in PM<sup>1</sup> originated predominately from diesel or gasoline emission. The contributions of the individual PAHs to the sum of total PAH mass concentration were calculated as well and are presented in Figure 5.

**Figure 5.** The contributions of the individual PAHs to the sum of total PAH mass concentration during: (**a**) cold and (**b**) warm period.

From Figure 5 it is evident that the PAH mixture was characterized by a high contribution of 6-ring (BghiP, IP) and 5-ring PAHs (BaP, BjF, BbF, BkF) characteristic for vehicle exhaust emission. Contributions of Flu and Pyr (markers for biomass burning) in the cold period were similar to warm period and it was 5% in cold and 7% in warm period. Previous investigations at the same measuring site but in the PM<sup>10</sup> particle fractions showed similar results during heating season [26]. These results suggest that four-ring PAHs concentrations have similar dominant sources during the warm and cold period (probably emission from gas and vehicle exhaust and domestic heating). A higher result than those measured in this study was found in Sarajevo [47] and in Poland [40]. Results in Sarajevo and Poland showed a higher contribution of four–ring PAHs (Flu, Pyr, and BaA), which were emitted from coal and wood combustion. The principal compound was BbF (16%) followed by BghiP (13–16%), IP (12–14%), and Chry (9–12%) in both periods. All of them were of pyrolytic origin, which suggests that PAHs reacted at similar extents and the dominant sources were similar.

However, Katsoyiannis et al. [34] presented that the diagnostic ratio cannot be effective as markers of pollution sources, because they can be influenced by atmospheric reaction of PAHs and weather conditions. Because of the fact that the diagnostic ratio can lead to mistakes in estimating the possible pollution sources, in this paper, PCA analysis was used in order to determine the most probable main pollution sources.

#### *3.3. Principal Component Analysis*

Furthermore, for investigating the similarities and differences between samples, PCA was applied, both for the cold and warm period. In the cold and warm period, the eigenvalues of the first two PCs were larger than 1, which indicated their significance. In the cold period, the first two principal components represented 98.08 cum.%, with the first and the second PCs contributing to the total variance of the data set with 89.09 cum.% and 8.99 cum.%, respectively. Results of the PCA during the cold period are presented on the PC1-PC2 loading (Figure 6a) and score plots (Figure 6b), illustrating the orientation of the variables and samples with respect to these principal components.

All variables had a negative effect on PC1, while all PAHs had positive values on PC2 except for Flu, Pyr, and Chry. Chry showed the lowest negative values. The highest negative effect on PC2 was for Pyr and Flu. In the warm period the first two principal components represent 97.74 cum.% whereas the first and the second PCs contributed with 95.20 cum.% and 2.54 cum.%, respectively, to the total variance of the data set. Results of the PCA during the warm period are presented on the PC1-PC2 loading (Figure 7a) and score plots (Figure 7b), illustrating the orientation of the variables and samples with respect to these principal components.

**Figure 6.** PC1-PC2 loading (**a**) and score plot (**b**) during cold period. I-1 to I-365 are the numbers of days.

**Figure 7.** PC1-PC2 loading (**a**) and score plot (**b**) during warm period. I-1 to I-365 are the numbers of days.

All of the variables also had a negative effect on PC1 as in the cold period, while Flu, Pyr, and Chry exhibited positive values on PC2, but the rest of PAHs showed negative values. Based on these parameters, two different factors were extracted during both periods. In the cold period, Flu and Pyr were grouped separately because their common origin was different from the rest of the PAHs. These are PAHs characteristic for biomass burning, and their sources probably could have been domestic heating [4,28,48,49]. In the warm period, these two PAHs (Flu and Pyr) were extracted separately probably because they evaporated easily from particles, which makes them more sensitive to meteorological conditions than the other PAHs [50]. However, the results obtained by PCA analysis confirmed those obtained by the diagnostic ratio, but with limitations in the warm period when processes such as meteorological condition and evaporation come into play. PCA analysis secluded some samples during the warm and cold period. From the PC1-PC2 loading, it is evident that six samples during the warm period (Figure 7b) and seven samples during the cold period (Figure 6b) were secluded from other samples. Results of contribution of the individual PAHs to the sum of total PAH for only those selected days showed some differences from the contribution made for all measured warm and cold periods. During warm period samples, I-271 and I-272 had the highest contribution of Flu and Pyr, but samples I-93, I-269, I-270, and I-273 had the highest contribution of BaP except for the I-273 sample, which also had the highest contribution of Chry. During the cold period, samples I-55-59 had the highest contribution of Flu and Pyr, and samples I-354 and I-362 had the highest contributions of BaP, Chry, and BaA. At the end of February and beginning of March (samples I-58 to I-59), the air temperature was extremely lower (from -2 to -10 ◦C) and probably emission from domestic heating was higher, which caused a higher contribution of Flu and Pyr to the total PAH concentration. For days at the end of December (I-354 and I-362), the air temperature was not as extremely low as in February but nevertheless below zero, and these days are also celebratory days, which can cause high concentrations of PAHs. For the warm period, at the end of September (samples I-269 to I-273) the air temperature was approximately 10 ◦C, which is much lower than the average value for the warm period (20 ◦C). Results of PCA analysis were in very good agreement with diagnostic ratios, and for those days, the diagnostic ratios showed mixed sources (diesel/gasoline and wood combustion).

PCA analysis extracted two factors during the cold and warm periods: the first factor separated Flu and Pyr, and the second factor the rest of the PAHs. Because of that, wind roses were shown for two groups, the sum of Flu and Pyr concentrations (ΣFlu+Pyr) and the sum of the remaining nine PAHs (ΣRest PAHs). The dependence of PAH concentrations on wind direction is shown on Figure S1 of the Supplementary Data section together with wind frequencies and wind velocities.

Winds coming from ENE were the most frequent during the cold and warm periods, followed by winds from NNE and SE-ESE sector, while the winds from other directions were relatively rare. On the other hand, winds of high velocities came from north-western and south-western directions, as well as from NE-ENE and SE sections. In the cold months, the highest concentrations of ΣFlu+Pyr primarily came from the west (residential part with domestic heating and street with modest density) and then from the south (center of the town, with dense traffic, while in the warm months the highest concentrations came from the east (industrial part of the town). The difference in concentrations in the cold and warm part of the year is greater than 80%. For the ΣRest PAH wind roses look similar as for ΣFlu+Py, but in the cold months of the year the dominant direction was SSE and then W, while in the warm months the dominant source of pollution was the east (Figure S2). In the cold periods of the year, pollution concentrations from the east were completely absent, which indicated that in the cold months the dominant source of pollution came mostly from the west for ΣFlu+Pyr while for ΣRest PAH it can be concluded that the dominant source of pollution was located both in the west and in the SSE direction. Taking into account the general absence of winds from west direction it can be concluded that the dominant sources of PAHs are probably of local origin. In the warm months of the year, there was no direction from which there was no pollution.

#### *3.4. Carcinogenic Activity of the PAH Mixture and Health Impact*

The BaP equivalent concentration and total carcinogenic potency were calculated according to Equation (1) and are shown in Table 4. In the cold period of the year, the TCP was more than ten times higher than in the warm period. The average TCP in the warm period was 0.141 ng m−<sup>3</sup> and in the cold 2.139 ng m−<sup>3</sup> . These results were similar to the TCPs obtained previously in Zagreb [26]. In a study by Pehnec and Jakovljevi´c [26], TCP values were shown for 30 samples collected per season (spring, summer, autumn, and winter) and TCP values varied from 0.063 ng m−<sup>3</sup> during summer to 4.503 ng m−<sup>3</sup> during winter. TCP values similar or lower than the ones in this study were reported by other authors who used the same TEF scheme [20,28], while higher TCP values were reported for some locations in Poland [22,51] and the Czech Republic [24].


**Table 4.** Equivalent BaPeq (ng m−<sup>3</sup> ) concentrations for individual PAHs measured during the warm and cold period.

(a) according to Nisbet and LaGoy (1992) (b) total carcinogenic potency.

Contributions of BaP to the carcinogenic potency in this study exceeded 60%. These values were similar for both measuring periods and were on average 62% and 67% for the warm and cold season, respectively. This indicates that benzo(a)pyrene had an important role in the carcinogenic activity of the PAH mixture. The same contribution was reported by Jakovljevic et al. [29] for the same location. Many authors reported similar BaP contributions but in other (PM<sup>10</sup> or PM2.5) particle fractions [46,52].

To determine the daily population exposures, total carcinogenic potency was used to calculate the daily dose according to Equation (3). In this study, we tried to estimate the most probable scenario for three age groups: infant (0–1 year), children (5–19 year), and adult (20–70 year). If we assumed that people spent an average ten hours on their job/school and eight to ten hours at home (including sleeping), we can assume that they were elsewhere for the rest of their day, due to the fact that people spend approximately 25% of their time outdoors.

The incremental lifetime cancer risk (ILCR) posed by exposure to PM-bounded PAHs was computed following Equation (2) [20].

We used Equation (3) for calculations of IEL for the warm and cold period, for three age group and the results are shown in Table 5.


**Table 5.** The BaPeq daily dose (IEL) and the incremental lifetime cancer risk (ILCR) for three age groups.

As the variables in the ILCR calculation belong to certain distributions taking the average value into account can be misleading for yielding the risk estimation. Therefore, we employed a probabilistic Monte Carlo (MC) simulation to feed the distributions to the ILCR model for infants, children, and adults in cold and warm seasons. All of the sampled variables were sampled randomly and independently. The information on variable distributions was taken from Chen and Lia [17] and Liu et al. [18]. We sampled the TCP from our own data after log-transforming them, separately for warm and cold seasons. The MC was iterated 10,000 times per sample variable and finally for the IEL

and ILCR models. The parameter values used for MC are show in Table S3 of Supplementary Data. The average ILCR value and the MC simulation results for ILCR were than compared.

IEL was lower during the warm than cold period; 5.04 ng day−<sup>1</sup> , 21.0 ng day−<sup>1</sup> , and 27.7 ng day−<sup>1</sup> in the warm period for infants, children, and adults, respectively. In the cold period, the BaPeq daily dose was higher; 76.5 ng day−<sup>1</sup> , 319.2 ng day−<sup>1</sup> , and 420.2 ng day−<sup>1</sup> for infants, children, and adults, respectively. Exposure time is a parameter that influences the inhalation risk most strongly. The incremental lifetime cancer risk (ILCR) was calculated according to Equation (2) and it was 3.35 × 10−<sup>10</sup> for infants, 4.56 × 10−<sup>9</sup> for children, and 1.08 × 10−<sup>8</sup> for adults in the cold period, but in the warm period the ILCR was 15 times lower than in the cold. According to US EPA [53], a one in a million chance of developing an additional human cancer over a lifetime (ILCR = 10−<sup>6</sup> ) is considered to be an acceptable level risk. A lifetime risk of one in a thousand (ILCR = 10−<sup>3</sup> ) was considered to be a serious health threat. The results determined in this study were much lower than the acceptable risk level of 10−<sup>6</sup> . Higher ILCR values were found by Majewski et al. [20], who calculated the ILCR for students and lecturers at Gliwice and Warsaw University. They found that the ILCR of students exposed to PAHs bounded to the PM<sup>1</sup> particle ranged from 5.49 × 10−<sup>8</sup> (Warsaw) to 1.43 × 10−<sup>7</sup> (Gliwice). However, these results were also below the acceptable risk level of 1 × 10−<sup>6</sup> .

#### *3.5. Influence of Meteorological Conditions and Gaseous Pollutants on PAH Concentrations*

Spearman's correlation matrix (Tables 6 and 7) showed that most of the PAHs measured in the cold and warm periods correlated well with each other (*p* < 0.05). The correlation between PAHs was very strong and coefficients ranged from 0.52 to 0.99 in the cold and from 0.81 to 0.99 in the warm period. These very high correlation coefficients suggested their shared sources and activities. Negative correlations were established between all of the PAHs and temperature and between PAHs and rainfall in both measuring periods (cold and warm). Humidity and pressure were negatively correlated with concentrations of PAHs in the warm period. During the cold part of the year, long temperature inversion periods may occur. Such weather conditions are characterized by high pressure and impaired air mixing, which favors pollutant accumulation. These stable atmospheric conditions combined with increased emissions are the basic origins of elevated concentrations of PAHs during the cold season.

The negative correlation with temperature can be explained with the ability of PAHs to easily vaporize from particle to gas phase [50,54]. A significant correlation was found in the cold period between CO, SO2, and all PAHs, as was a positive correlation between NO<sup>2</sup> and PAHs with high molecular weight, while a negative weak correlation but still statistically significant was established between heavy PAHs and O3. Similar correlations were found during the warm period except between PAH and SO2, where the correlation was not significant. As there is no industry at the location, SO<sup>2</sup> was mostly emitted from vehicles [55]. Previous studies have shown that SO<sup>2</sup> concentrations have been continuously decreasing over the last ten years, the reason for which is the fact that coal has not been in use for more than thirty years at this location [56]. Positive correlations between PAHs and NO<sup>2</sup> were in good agreement with other studies [55,57]. Previous studies have shown that NO<sup>2</sup> is one of the main traffic-related air pollutants, and in reaction with PAHs, it can produce nitro PAHs in the atmosphere [55]. In both measuring periods, a negative weak correlation but still significant was found between ozone and PAHs. Ozone is a reactive atmospheric pollutant, and it is able to, in the presence of heat and sunlight, react with PAHs so this was probably one of the reasons for the low PAHs concentrations during summer [55,58]. Positive correlations between CO and PAH suggested their shared sources and activities in both the cold and warm period. Similar results were found in other studies [59].


**Table 6.**Spearman's correlation coefficients between PAHs, gaseous pollutants, and some meteorological parameters during the cold period.

t—temperature, RH—relative humidity, and P—atmospheric pressure. Statistically significant correlation coefficients (*p*<0.05) are underlined.


**Table 7.**Spearman's correlation coefficients between PAHs, gaseous pollutants, and some meteorological parameters during the warm period.

t—temperature, RH—relative humidity, and p—atmospheric pressure. Statistically significant correlation coefficients (*p*<0.05) are underlined.

#### **4. Conclusions**

Measurements of PAHs in the PM<sup>1</sup> particle fraction at an urban location in Zagreb, Croatia, showed seasonal differences with higher mass concentrations of PAH in the cold than in the warm period. The annual mass concentration for BaP was 0.765 ng m−<sup>3</sup> , which indicated that the value set by Directive 2004/107/EC for BaP of 1 ng m−<sup>3</sup> has not been exceeded. During both measurement periods, the contributions of heavy PAHs, characteristic for vehicle emissions, were much higher than those of light PAHs. Results of diagnostic ratios and PCA showed that the main emission source of PAHs associated with PM<sup>1</sup> in this study area was engine combustion of diesel or gasoline during the warm period but that in the cold period emission from domestic heating is, in addition to diesel, the dominant source. The total carcinogenic potency was estimated using toxic equivalency factors. The average TCP was 0.14 ng m−<sup>3</sup> and 2.14 ng m−<sup>3</sup> for the warm and cold period, respectively. The incremental lifetime cancer risk (ILCR) was determined for three age groups of the population: infants (0–1 year), children (5–19 year), and adults (20–70 year) and was below the maximum acceptable level (1 × 10−<sup>6</sup> ), revealing that carcinogenic risk posed to the three age groups via inhalation is acceptable. However, more than ten times higher PAH values in the cold part of the year, as well as associated health risk, emphasize the need for regular monitoring of PAH levels in smaller particle fractions, such as PM1. These data are a valuable tool in future plans and actions to control PAH sources and to improve the quality of life of urban populations.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1660-4601/17/24/9587/s1, Table S1: Detection (DL) and quantification limits (QL) and recoveries (R) of PAHs analyzed by high-performance liquid chromatography (HPLC) with a fluorescence detector, Table S2: Comparison of PAH mass concentrations in PM1 particle fraction obtained in this study with other studies, Table S3: Risk parameters for different age groups, Figure S1: Wind roses. Dependence of PAH concentrations (ng m−<sup>3</sup> ) on wind direction, average wind velocities (m s−<sup>1</sup> ) and wind direction frequencies (%) for cold and warm measuring periods, Figure S2: Distribution of incremental lifetime cancer risk for adults, children and infants, derived using Monte Carlo simulation in cold and warm period.

**Author Contributions:** Conceptualization, I.J. and G.P.; data curation, I.J. and G.P.; validation, S.D., I.B. and M.L.; formal analysis, I.J. and Z.S.Š.; writing—original draft, I.J.; writing—review and editing, I.J., G.P., R.G. and M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** These measurements were conducted within the internal scientific project "Organic content of PM<sup>1</sup> particle fraction" funded by Institute for Medical Research and Occupational Health (PI: R. Godec).

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Statistical Forecast of Pollution Episodes in Macao during National Holiday and COVID-19**

**Man Tat Lei 1,2,\* , Joana Monjardino <sup>3</sup> , Luisa Mendes <sup>1</sup> , David Gonçalves <sup>2</sup> and Francisco Ferreira <sup>3</sup>**


Received: 27 May 2020; Accepted: 14 July 2020; Published: 15 July 2020

**Abstract:** Statistical methods such as multiple linear regression (MLR) and classification and regression tree (CART) analysis were used to build prediction models for the levels of pollutant concentrations in Macao using meteorological and air quality historical data to three periods: (i) from 2013 to 2016, (ii) from 2015 to 2018, and (iii) from 2013 to 2018. The variables retained by the models were identical for nitrogen dioxide (NO2), particulate matter (PM10), PM2.5, but not for ozone (O3) Air pollution data from 2019 was used for validation purposes. The model for the 2013 to 2018 period was the one that performed best in prediction of the next-day concentrations levels in 2019, with high coefficient of determination (R<sup>2</sup> ), between predicted and observed daily average concentrations (between 0.78 and 0.89 for all pollutants), and low root mean square error (RMSE), mean absolute error (MAE), and biases (BIAS). To understand if the prediction model was robust to extreme variations in pollutants concentration, a test was performed under the circumstances of a high pollution episode for PM2.5 and O<sup>3</sup> during 2019, and the low pollution episode during the period of implementation of the preventive measures for COVID-19 pandemic. Regarding the high pollution episode, the period of the Chinese National Holiday of 2019 was selected, in which high concentration levels were identified for PM2.5 and O3, with peaks of daily concentration exceeding 55 µg/m<sup>3</sup> and 400 µg/m<sup>3</sup> , respectively. The 2013 to 2018 model successfully predicted this high pollution episode with high coefficients of determination (of 0.92 for PM2.5 and 0.82 for O3). The low pollution episode for PM2.5 and O<sup>3</sup> was identified during the 2020 COVID-19 pandemic period, with a low record of daily concentration for PM2.5 levels at 2 µg/m<sup>3</sup> and O<sup>3</sup> levels at 50 µg/m<sup>3</sup> , respectively. The 2013 to 2018 model successfully predicted the low pollution episode for PM2.5 and O<sup>3</sup> with a high coefficient of determination (0.86 and 0.84, respectively). Overall, the results demonstrate that the statistical forecast model is robust and able to correctly reproduce extreme air pollution events of both high and low concentration levels.

**Keywords:** air pollution; air quality forecast; modelling; pollution episodes; national holiday; COVID-19

#### **1. Introduction**

The development of air quality forecast models is essential for cities with high population density, including Macao, one of the most densely populated cities in the world. It is extremely important to predict pollution episodes so the authority can provide a warning to the local community in advance to avoid the adverse air quality, which may lead to severe health consequences. In order to predict

next-day concentrations of nitrogen dioxide (NO2), particulate matter (PM<sup>10</sup> and PM2.5), and maximum hourly concentration of ozone (O3 MAX) for roadside, ambient, and residential stations in Macao, a forecast model was developed based on statistical methods using multiple linear regression (MLR) and classification and regression tree (CART) analysis.

There are three forms of total suspended particles (TSPs), which include coarse, fine, and ultrafine particles. Coarse particles, also known as PM10, are derived from suspension of dust, soil, sea salts, pollen, mold, and other crustal materials. Fine particles, also known as PM2.5, are derived from emissions from combustion process, including vehicles powered by petrol and diesel, wood burning, coal burning, and other industrial processes. Ultrafine particles are derived from combustion related sources such as vehicle exhausts and atmospheric photochemical reactions [1].

O<sup>3</sup> is the most important index substance for photochemical smog, one of the major air pollutants [2]. The formation of ground-level O<sup>3</sup> heavily depends on the concentration levels of volatile organic compounds (VOCs) and nitrogen oxides (NOx) and meteorological factors such as wind speed, insolation, and temperature. PM2.5 and O<sup>3</sup> pollutants are known to cause the most damages to the human respiratory and cardiovascular system. A study for Terengganu State, Malaysia, showed that high levels of O<sup>3</sup> occurring under dry and warm conditions during the southwest monsoon, were higher in industrial areas, and were positively correlated with the maximum daily temperature [3].

The emission of NO<sup>x</sup> is primarily emitted from transportation and combustion process, while the emission of VOCs is primarily emitted from road traffic and the use of products containing organic solvents [4,5].

The emission of NO<sup>x</sup> and VOCs is responsible for the O<sup>3</sup> formation, in particular rural areas being NOx-sensitive while urban areas being VOC-sensitive. Nevertheless, the greater NO<sup>x</sup> emission reductions have contributed to a widespread shift in the O<sup>3</sup> production regime from NOx-saturated (high-NOx) to NOx-sensitive (low-NOx) in some urban areas, while O<sup>3</sup> production in rural areas is even more sensitive to NOx.

TSPs are primary contributors to premature death worldwide, with over four million premature deaths being recorded due to exposure to high levels of ambient PM2.5 [6–8]. PM2.5 can penetrate deep into the lungs when being inhaled, which leads to both acute and chronic health issues [1,6]. NO<sup>2</sup> and TSPs are responsible for 412,000 and 71,000 premature death per year, respectively, in the European Union [9,10]. Moreover, previous studies show a strong correlation between short-term exposure to NO<sup>2</sup> and both the number of hospital outpatients with eye and adnexa diseases (EADs) [11] and the number of hospital admission due to cardiovascular diseases (CVD) [12]. The Chinese National Ambient Air Quality Standard (NAAQS) has set the threshold of PM10, PM2.5, and O3 MAX concentration at 150 µg/m<sup>3</sup> , 75 µg/m<sup>3</sup> , and 160 µg/m<sup>3</sup> , respectively, while the WHO Air Quality Guideline has set the same thresholds at 50 µg/m<sup>3</sup> , 25 µg/m<sup>3</sup> , and 100 µg/m<sup>3</sup> , respectively. Compliance with the thresholds set by the WHO for PM2.5 could improve life expectancy in China by 0.14 years [13] and ambient air pollution has caused at least 3.7 million deaths, with more than 25% of deaths in Southeast Asia [14,15].

Air pollution forecasting models can provide important information for populations to adopt mitigation measures during high pollution days. To be useful, these models should be robust to deal with extreme variations in pollution levels, in particular during high-pollution peak days. Factors leading to extreme variation in pollution levels are diverse and include both human activities and meteorological factors.

In a study for Beijing, China, the reduction of traffic flow and vehicle emissions in downtown areas during the Chinese National Holiday, reduced air pollution, while, in contrast, fireworks during the Chinese New Year Holiday had the opposite effect [16]. When highway tolls were being waived for passenger vehicles during the Chinese National Holiday across the entire nation of China, air pollution increased by 20% and visibility decreased by 1 km, causing economic losses due to negative health impacts estimated at RMB 0.95 billion [17]. Nevertheless, the Chinese National Holiday is known to be a golden week of tourism, in which the Chinese tourist flock to different tourist destinations around the world to celebrate the national holiday. Due to the vibrant casinos and entertainment industry and close proximity to mainland China, Macao is also one of the favorite destinations for Chinese tourists, so the influx of tourist during the period of Chinese National Holiday may lead to an increase of emissions in Macao.

Likewise, the recent COVID-19 crisis has had an extreme impact in air pollution levels. The Wuhan Health Commission has first reported cases of pneumonia linked to the Wuhan wet market in Hubei Province, China, back in December 2019 [18]. Preventive measures were implemented soon after that abruptly reduced industrial activities and transportation. Nevertheless, the levels of air pollutants, in particular of PM2.5, remained severe in northern China throughout the end of January 2020 due to adverse meteorological conditions that have overwhelmed the benefits of emission reduction in transportation and industrial sectors [19].

Previous work showed that there is an increase in the level of O<sup>3</sup> concentrations and a decrease in the level of NO2, PM10, and PM2.5 concentration during the period of COVID-19 pandemic lockdown in several cities of China, due to the significant reduction of transportation and industrial activities [4,5,20,21].

Several methodologies have been developed and applied to forecast air quality across the world, including deterministic, statistical, and machine learning methods [22–26]. Some studies showed that statistical models are more accurate and efficient compared to deterministic models, particularly in regions with complexed terrain [27–30] Moreover, prediction of NO2, PM10, PM2.5, and O3 MAX concentrations based on MLR and CART models have been successfully implemented in Macao, Bangkok, Changsha City, Beijing, Bilbao, and Pakistan [26,31–35].

In this context, it is relevant to develop a reliable methodology to forecast the concentration of air pollutants, which is presented and tested for a high pollution episode (associated with the Chinese National Holiday) and a low pollution episode (during COVID-19 preventive measures).

#### **2. Materials and Methods**

The air quality and meteorological variables that were considered to build all of the air quality statistical models were obtained from Macao Meteorological and Geophysical Bureau (SMG). The air quality data was gathered from the air quality monitoring network, namely for: Macao Roadside, Macao Residential, Taipa Ambient, Taipa Residential, and Coloane Ambient stations, which have a suitable historic dataset of surface air quality measurements for the levels of NO2, PM10, PM2.5, and O<sup>3</sup> concentrations. These background stations (residential and ambient) can capture the regional contribution of PM<sup>10</sup> and PM2.5. There is a higher population and traffic density in Macao Roadside and Macao Residential, which are located in the main peninsula, in comparison to Taipa Ambient, Taipa Residential, and Coloane Ambient stations, which are located on the outlying islands.

Meteorological data was obtained from surface observations at SMG's Taipa Grande Meteorological Station, hourly observations from automatic weather stations, such as temperature, relative humidity, precipitation, average wind speed, and dew point temperature, as well as upper-air observations (from Hong Kong King's Park location) such as geopotential heights, thickness, stability, temperature, relative humidity, and dew point temperature at various altitudes. In the present work, statistical models such as multiple linear regression (MLR), and classification and regression tree (CART), are developed, based on historical measurements of meteorological and air quality variables. Table 1 presents all the variables considered as predictors in the MLR and CART forecast models, as shown in previous work [22]. The air quality variables considered included the levels of NO2, PM10, PM2.5, and O3 MAX concentration from 00:00 to 23:00 of the previous day, two days and three days ago, and from 16:00 of the previous day and 15:00 of today. The meteorological variables being considered included the upper-air observations from King's Park location, Hong Kong Observatory, surface observations and other variables from the monitoring network of Macao Meteorological and Geophysical Bureau (SMG). **Table 1.** Variables considered as predictors in the multiple linear regression (MLR) and classification and regression tree (CART) models in all of the air quality forecast models.


Meteorological variables: \* Daily sounding at 12H (GMT+8) at King's Park Meteorological Station—Hong Kong Observatory.

In this study, meteorological and air quality variables for 2013 to 2016, 2015 to 2018, and 2013 to 2018 were used to build three separate forecasting models. The 2013 to 2016 model was constructed for the initial evaluation for the application of the statistical model to forecast air quality in Macao, while the 2015 to 2018 models and the 2013 to 2018 models are a follow-up, to determine if any improvement could be made with two additional years of data. The comparison of extended data ranging from 5 to 6 years are considered to be adequate lengths to test if there is any significant difference between the time series. Simultaneously, it would not be ideal to trace back too far with the time series, because regional emissions are constantly changing, and therefore the level of pollutants concentration may also be changing. The dataset from 2019 was the most recent dataset, which would be used for the model validation for all the models. This study is an empirical approach and also region-specific, which may also be chemical-regime dependent.

The final selected variables to predict the levels of PM2.5 and O<sup>3</sup> concentration are common to different locations of Macao air quality monitoring stations. Some variables initially selected were rejected from the forecast models due to collinearity. The final objective is to obtain prediction models with the lowest number of variables, but with the maximum explained variance as translated by the coefficient of determination (R<sup>2</sup> ).

After selecting the best model, it was applied to forecast pollution levels during an extremely high pollution episode, and a low pollution period. The high and low pollution selected episodes were, respectively: (i) the period of Chinese National Holiday, a week before the Chinese National Holiday from September 23rd to 30th, 2019, and the week during the Chinese National Holiday from October 1st to 7th, and (ii) the preventive measures period of COVID-19, from February 5th to 20th, 2020.

The statistical model was built using IBM SPSS Statistics version 26 with MLR (stepwise) and CART methods [26,36]. SPSS is a statistical software that is applied to solve research problems through hypothesis testing and predictive analysis.

Model performance indicators were calculated, such as, coefficient of determination (R<sup>2</sup> ), root mean square error (RMSE), mean absolute error (MAE), and systematic error (BIAS).

#### **3. Results and Discussion**

#### *3.1. Air Quality Forecast Models*

The model performance indicators obtained for the 2013 to 2016 model and for the 2013 to 2018 model, validated with 2019 data, are listed in Tables 2 and 3, respectively. The models chosen to figure in Tables 2 and 3 are the ones that performed the worst and best 2019 validation results.

The results showed that the model for the 2013 to 2018 period was the one that performed best in predicting next-day concentrations levels in 2019, with high R<sup>2</sup> between predicted and observed daily average concentrations (between 0.78 and 0.89 for all pollutants) and low RMSE, MAE, and BIAS. The additional two years of data helped to improve the air quality forecasting model. Nevertheless, with the two other models (2013–2016 and 2015–2018) a significant R<sup>2</sup> (between 0.78 and 0.89 for all pollutants) was also obtained, but it translated into a less reliable air quality forecast.

Regarding model performance indicators obtained per pollutant and station, the majority of models show a good agreement and a similar R<sup>2</sup> range values (from 0.81 to 0.89), except for O3 MAX, which is more difficult to predict. MLR was used for all pollutants, while CART analysis was used in almost all the O3 MAX models (Tables 2 and 3). This CART analysis complement was an approach to obtain improved results, mainly regarding a better prediction of high pollutant levels.

Table 4 presents the final model equations obtained for each pollutant, per air quality monitoring station, in the 2013 to 2018 model. Additionally, the final equations used to predict the levels of NO2, PM10, PM2.5, and O3\_MAX concentrations are presented in Table 4.

#### *Int. J. Environ. Res. Public Health* **2020**, *17*, 5124


#### **Table 2.** Model performance indicators for the 2013 to 2016 model validation with 2019 data.


**CART** 

**CART** 

**CART** 

**CART** 

**CART** 

**CART** 

−

−

−

−

− 

−

−

−

−

−

−

−

−

−

−

−

−

−

−

− −

−

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 − 

 − 

 − 

 − 

 

− 

−

−

−

−

− 

−

−

−

− 

−

−

−

−

−

−

<sup>−</sup>

−

−

−

 

 

 

 −

 

 

 

 

 

 −

 

 

 

 

 

 

 

 

 

 

 

 

−

−

−

−

 

 

> 

 

 − 

 

 

 

−

−

−

−

−

−

−

−

−

−

−

−

−

−

<sup>−</sup>

−

<sup>−</sup>

−

−

−

− −

−

−

−

<sup>−</sup>

− −

−

−

−

−

−

−

−

−

−

− −

−

−

−

−

−

−

−−

−

−

−

−

−

−

−

<sup>−</sup>

−

<sup>−</sup>

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−


**Table 4.** Variables and model equations for each pollutant per air quality monitoring station in the 2013 to 2018 model.


#### *3.2. Air Quality During the High Pollution Episode*

Taipa Ambient is the representative background location for Macao, and was chosen to assess the background levels of PM2.5 and O<sup>3</sup> during the extreme pollution episode.

The influx of tourists coming to Macao, in light of the Chinese National Holiday, contributed to an high pollution episode that occurred during late September and early October 2019, with peak daily levels of PM2.5 concentration exceeding 55 µg/m<sup>3</sup> and O3 MAX levels exceeding 400 µg/m<sup>3</sup> , largely exceeding the threshold level recommended by the WHO.

*Int. J. Environ. Res. Public Health* **2020**, *17*, 5124

The levels of PM2.5 and O3 MAX concentrations for Taipa Ambient during the Chinese National Holiday in 2019 (from September to November) are presented in Figures 1 and 2.

**Figure 1.** PM2.5 concentrations for Taipa Ambient highlighting a pollution episode immediately before, and during, the Chinese National Holiday of 2018 and 2019 (September to November).

**Figure 2.** O3 MAX concentrations for Taipa Ambient highlighting a pollution episode immediately before, and during, the Chinese National Holiday of 2018 and 2019 (September to November).

μ μ

μ μ μ

Figures 1 and 2 showed the comparison of daily average PM2.5 and O3 MAX concentration during 2018 and 2019, from a month before in September and a month after in November of the Chinese National Holiday. The pollution episode of 2019 occurred just before and going well into the period of Chinese National Holiday (1 to 7 October).

As shown in Figures 1 and 2, the levels of PM2.5 and O3 MAX concentration peaked immediately before, and during, the Chinese National Holiday in late September and early October 2019. The monthly mean concentration of PM2.5 (from September to November) during the Chinese National Holiday in 2019 was 19 µg/m<sup>3</sup> , 24 µg/m<sup>3</sup> , and 28 µg/m<sup>3</sup> , respectively. In addition, the monthly mean concentration of O3 MAX (from September to November) during the Chinese National Holiday in 2019 was 181 µg/m<sup>3</sup> , 163 µg/m<sup>3</sup> , and 172 µg/m<sup>3</sup> , respectively.

The levels of O3 MAX concentrations reached its peak during the late September and early October due to meteorological factors including predominant winds from the north and east, from the Guangdong Province and Hong Kong, respectively. Temperatures were high in conjunction with low wind speed. The average daily temperature during the ozone peak episode that took place the two-weeks before the Chinese National Holiday (October 1st) was 28 ◦C, while the maximum daily average was 31 ◦C. Average wind speed was 2.5 m/s.

Due to the shutdown of nearby industrial sectors during the period of Chinese National Holiday, there were lower emissions of nitrogen oxides associated with the decreased load from the coal power plants in the northern region, usually supporting the operation of the factories. Therefore, this caused a decrease NOx, the precursor of O3. However, the increase in emissions of VOCs and NOx by vehicles, with chemical reactions in the presence of sunlight, may have caused the peak levels of ozone concentrations under these high temperature favorable conditions.

#### *3.3. Air Quality During the Low Pollution Episode*

In contrast, the COVID-19 pandemic has led to the Macao government's decision to temporarily suspend the operation of the casinos and entertainment industry and highly restrict cross border movements, as a preventive measure to reduce population mobility within the region of Macao. As a result, it has caused a low pollution episode during late January and early February 2020, with daily levels of PM2.5 concentration reaching a record low at 2 µg/m<sup>3</sup> and O3 MAX levels at 50 µg/m<sup>3</sup> . The reduction of population mobility, and consequently, of traffic emissions in Macao and its nearby Guangdong Province, lead to this lowest PM2.5 concentration levels.

As shown in Figure 3, the levels of PM2.5 concentrations remained low during the initial outbreak of COVID-19 pandemic in Macao (from January to February 2020), slowly recovering to pre-COVID-19 values in March 2020. As shown in Figure 4, the levels of O3 MAX concentration remained high during the initial outbreak of COVID-19 pandemic in Macao (from January to February 2020) and the high levels continued into March 2020. The higher levels of O3 MAX concentration were associated with lower NO<sup>X</sup> emissions, which led to a weakened O<sup>3</sup> titration by NO during the COVID-19 pandemic lockdown in the nearby Guangdong Province [4].

Despite industrial emission being a major contributor to the PM2.5 pollution in China prior to COVID-19 pandemic lockdown period, the residential emission contributed to 39% of total PM2.5 emissions in China, so the emissions of PM2.5 during the lockdown period may have originated from residential areas [5].

The comparison of PM2.5 and O3 MAX concentrations for Taipa Ambient during the previous year of 2019 and COVID-19 pandemic in 2020 (January to March) is presented in Figures 3 and 4.

**Figure 3.** Comparison of PM2.5 concentrations for Taipa Ambient during the previous year of 2019 and COVID-19 pandemic in 2020 (January to March).

**Figure 4.** Comparison of O3 MAX concentrations for Taipa Ambient during the previous year of 2019 and COVID-19 pandemic in 2020 (January to March).

μ μ μ

μ μ μ

As shown in Figure 5, the difference between monthly mean concentration (from January to March) of PM2.5 concentration in 2019 and 2020 was 16 µg/m<sup>3</sup> , 2 µg/m<sup>3</sup> , and 1 µg/m<sup>3</sup> , respectively. As shown in Figure 6, the difference between monthly mean concentration (from January to March) of O3 MAX concentration in 2019 and 2020 was 12 µg/m<sup>3</sup> , 21 µg/m<sup>3</sup> , and 9 µg/m<sup>3</sup> , respectively.

μ μ μ

μ μ μ

**Figure 5.** Monthly mean PM2.5 concentrations for Taipa Ambient during the previous year of 2019 and COVID-19 pandemic in 2020 (January to March).

**Figure 6.** Monthly mean O3 MAX concentrations for Taipa Ambient during the previous year of 2019 and COVID-19 pandemic in 2020 (January to March).

The monthly mean concentration of PM2.5 and O3 MAX concentration for Taipa Ambient during the previous year of 2019 and COVID-19 pandemic in 2020 (January to March) is presented in Figures 5 and 6. Overall, the preventive measures of COVID-19 pandemic may not have caused a significant difference in the levels of PM2.5 and O<sup>3</sup> concentration in Macao, as the levels from February to March 2020 were similar to that of the previous year, 2019.

#### *3.4. Air Quality Pollution Episodes Discussion*

The air quality of Macao, a territory with only 32.8 km<sup>2</sup> , is heavily influenced by external factors, in particular by human activities that occur in the much larger and neighboring Guangdong province. Our study shows the extent to which an increase in mobility associated with Chinese National Holiday, or a decrease in the same factors, associated with the COVID-19 preventive measures period, impacts air quality in Macao.

The levels of PM2.5 concentrations significantly reduced after the first confirmed case of COVID-19 pandemic in Macao on January 22nd, 2020, which caused panic and anxiety in the local population, and continued by the announcement of casino closures by the Macao government as part of the preventive measures for COVID-19 from February 5th to 20th, 2020. As some of the preventive measures, in particular, the 15 days mandatory casino closure have been lifted, the fear and tension of the local residents has eased, which has promoted population mobility. Although the levels of PM2.5 concentrations in Macao improved significantly during late January and early February 2020, the levels of PM2.5 concentrations gradually returned to normal in March 2020 after some of the preventive measures began to be lifted in Macao and its nearby Guangdong Province.

#### *3.5. Air Quality Pollution Episodes Forecast*

Regarding the model behavior in predicting PM2.5 and O3 MAX during the high pollution episode (Chinese National Holiday), observed and predicted PM2.5 and O3 MAX concentrations are presented in Figures 7 and 8.

As shown in Figures 7 and 8, the levels of PM2.5 and O3 MAX concentration peaked during late September and early October 2019. The PM2.5 predicted levels followed the primary trend of the measured concentrations and followed the concentration peak represented in Figure 7. The model for O3 MAX also followed the primary trend, but it was more difficult to represent the concentration peak. The forecast model for PM2.5 has a higher R<sup>2</sup> in comparison to the model of O3 MAX, because the maximum hourly concentration of O3 MAX is more challenging to predict in comparison to the 24 h average of PM2.5, as there is influence from the regional precursors sources and also its complex chemistry with solar radiation for O<sup>3</sup> formation, which led to a higher degree of variability.

Due to the different nature of PM2.5 and O3 MAX, the forecast model performed better in the prediction of PM2.5 in comparison to O3 MAX. This can be demonstrated in the higher R<sup>2</sup> values in the PM2.5 forecast model. The observed and predicted PM2.5 and O3 MAX concentrations, during the low pollution episode (implementation of COVID-19 preventive measures), are presented in Figures 9 and 10.

The 2013 to 2018 model successfully predicted both the high and low pollution episodes, for PM2.5 and O3 MAX, obtaining a significant R<sup>2</sup> of 0.88 and 0.83, respectively, for the high pollution period (from September to November 2019), and an R<sup>2</sup> of 0.82 and 0.75, respectively, for the low pollution period (from January to March 2020). The R<sup>2</sup> obtained for the entire year of 2019 was 0.86 for both PM2.5 and O3 MAX. The statistical forecast model has been shown to be capable to predict, with a high coefficient of determination, the next 24 h.

**Figure 7.** Observed and predicted PM2.5 concentrations for Taipa Ambient during Chinese National Holiday (from September to November 2019).

**Figure 8.** Observed and predicted O3 MAX concentrations for Taipa Ambient during Chinese National Holiday (from September to November 2019).

**Figure 9.** Observed and predicted PM2.5 concentrations for Taipa Ambient during preventive measures of COVID-19 pandemic (from January to March 2020).

**Figure 10.** Observed and predicted O3 MAX concentrations for Taipa Ambient during preventive measures of COVID-19 pandemic (from January to March 2020).

#### **4. Conclusions**

As expected, the 2013 to 2018 model performed best with the highest R<sup>2</sup> and lowest RMSE, MAE, and BIAS as compared with the 2013 to 2016 model and the 2015 to 2018 model. The additional two years of data helped to improve the accuracy and stability of the forecast of the 2013–2018 model.

The 2013–2018 model was able to successfully predict the high pollution episode during the Chinese National Holiday in late September and early October 2019 and the low pollution episode during the preventive measures period of COVID-19 pandemic in late January and early February 2020. This shows that this model can be reliably applied to forecast next-day pollutants concentrations across different magnitude levels of air pollution, being a useful tool for mitigation of air pollution impacts.

In addition, this shows that an improvement of global air quality in the territory is possible but it is tightly linked to the implementation of air pollution control measures in the industry and mobility sectors in Macao, in particular, in Guangdong Province. As previously studied, the air pollution

problem associated with PM2.5 and O3 MAX is a regional problem that is not only limited to Macao, but also in the nearby regions of Hong Kong and Guangdong Province.

**Author Contributions:** Data curation, M.T.L., J.M., and L.M.; funding acquisition, D.G. and F.F.; methodology, M.T.L.; software, L.M.; supervision, D.G. and F.F.; validation, M.T.L., J.M., and L.M.; writing—original draft, M.T.L.; and writing—review and editing, M.T.L., J.M., D.G., and F.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [Fundação para a Ciência e Tecnologia, I.P., Portugal] grant number [UID/AMB/04085/2019] and the APC was funded by [CENSE].

**Acknowledgments:** The work developed was supported by The Macao Meteorological and Geophysical Bureau (SMG).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Five Year Trends of Particulate Matter Concentrations in Korean Regions (2015–2019): When to Ventilate?**

**Dohyeong Kim 1,**† **, Hee-Eun Choi 2,**† **, Won-Mo Gal <sup>2</sup> and SungChul Seo 2,\***


Received: 6 June 2020; Accepted: 2 August 2020; Published: 10 August 2020

**Abstract:** Indoor air quality becomes more critical as people stay indoors longer, particularly children and the elderly who are vulnerable to air pollution. Natural ventilation has been recognized as the most economical and effective means of improving indoor air quality, but its benefit is questionable when the external air quality is unacceptable. Such risk-risk tradeoffs would require evidence-based guidelines for households and policymakers, but there is a lack of research that examines spatiotemporal long-term air quality trends, leaving us unclear on when to ventilate. This study aims to suggest the appropriate time for ventilation by analyzing the hourly and quarterly concentrations of particulate matter (PM)10 and PM2.5 in seven metropolitan cities and Jeju island in South Korea from January 2015 to September 2019. Both areas' PM levels decreased until 2018 and rebounded in 2019 but are consistently higher in spring and winter. Overall, the average concentrations of PM10 and PM2.5 peaked in the morning, declined in the afternoon, and rebounded in the evening, but the second peak was more pronounced for PM2.5. This study may suggest ventilation in the afternoon (2–6pm) instead of the morning or late evening, but substantial differences across the regions by season encourage intervention strategies tailored to regional characteristics.

**Keywords:** particulate matter; natural ventilation; indoor air quality; regional variation

#### **1. Introduction**

The quality of the air we breathe every day, both outdoors and indoors, is a matter of great concern since it contains various polluting substances which have negative effects on human health and on the environment [1]. According to the 2017 Global Burden of Disease Study, exposure to outdoor air pollution is one of the leading risk factors for premature death, accounting for 3.4 million deaths each year [2]. One of the most critical air pollutants is particulate matters (PM) which are considered as one of major reasons for increased prevalence or exacerbation of respiratory diseases [3,4], cardiovascular diseases [5], and diabetes [6]. The International Agency for Research on Cancer (IARC) has designated the atmospheric PM as a carcinogen of the same class as asbestos and those found in tobacco smoke [7]. However, a growing number of studies have reported that indoor PM could be more harmful in comparison to outdoor air quality because potential air pollutants are built up in confined environments [8,9]. Pollutants are endlessly generated in indoor environments due to human activities and inflow of external air and the purification of contaminated indoor air quality is difficult [10]. Due to a modern lifestyle in which people stay indoors for longer in most industrialized countries [11], the health impacts of indoor PM have been reported to outweigh those of outdoor air pollution, particularly for long-term exposure even at low concentration levels [12,13]. This is

particularly relevant for susceptible groups such as children and elderly people who spend most of their time indoors [14,15].

Policymakers and professionals have developed strategies to manage indoor air quality such as improved ventilation and filtration, regulations and standard setting, routine monitoring and inspection, etc. [16]. Natural ventilation has been considered as the most economical and effective means of improving indoor air quality [17], which leads to improved productivity [18], mitigates allergic diseases [19], and even reduces the risk of airborne contagion [20]. Although natural ventilation may be less effective than mechanical ventilation in some settings [21], simply opening the windows for 10 min was found to be effective in achieving satisfactory indoor comfort conditions and air quality [22]. While natural ventilation is widely recommended for indoor air quality, its effectiveness depends on the quality of outdoor air [23,24]. Various concerns have been reported regarding its application, such as ambient particle concentration, indoor source intensity, indoor-outdoor air exchange rate, and circulation of outdoor air within the building, but the largest concern centers on the entry of polluted air from outdoors [25]. Due to the high filtration rate of particles from outdoors, the linkage between indoor and outdoor particle concentrations is more pronounced for a very small particle such as PM2.5 [26]. In particular, residential housing units are highly influenced by the quality of outdoor air since they tend to depend more on natural ventilation even if mechanical ventilation units are installed [27].

Considering a strong connectivity and dependency between indoor and outdoor air quality, natural ventilation could generate negative health impacts where the outdoor air quality is not acceptable for ventilating a building [28]. Despite the importance of risk-risk tradeoff between the exposure to indoor and outdoor air pollutants [29], except for a few governmental reports [17] there is no research that directly provides guidelines for the appropriate time and duration for natural ventilation based on the detailed assessment of temporal trends of ambient air pollutants in multiple communities. As hourly patterns of PM 2.5 and PM 10 could vary due to local characteristics such as major sources of contamination and seasonal weather conditions, an empirical investigation of long-term, site-specific air quality data is critical for developing evidence-based guidelines.

Despite the efforts of the Korean government to improve air quality during the past decades based on several laws and regulations such as the Clean Air Conservation Act of 1991 and the Special Act on the Improvement of Air Quality in Seoul Metropolitan Area of 2003, the problem of air pollution is far from being solved in South Korea. The public's concerns about the health risks associated with particulate matter exposure have substantially increased recently in South Korea, particularly after experiencing the record levels of outdoor PM 2.5 concentrations in many parts of the country in February and March 2019 [30]. Soon after these events, the national assembly of South Korea passed a series of bills which declared PM air pollution as a "social disaster" and provided emergency measures to tackle the problem, such as mandatory installation of air purifiers in classrooms, distribution of masks to vulnerable groups, and the promotion of low-emission vehicles [31]. It has been reported that Koreans spend more than 20 h a day indoors, particularly the most vulnerable groups such as patients, the elderly, and infants [32]. Although the Korea Ministry of Environment (KME) has provided several guidelines to citizens in response to indoor PM exposure, including ventilation, indoor water cleaning, and indoor air quality monitoring [33], they still lack details concerning how to actually implement the measures.

Only a few research articles have reported the seasonal patterns of particulate matter in South Korea, showing that PM10 and PM2.5 concentrations peak in spring and winter respectively, but both were lowest in summer [34]. However, to our knowledge there is no recent study in South Korea looking at hourly patterns of PM10 and PM2.5 concentrations by season in multiple regions based on long-term data. In South Korea, most people do not ventilate frequently due to Asian dust and air pollution problems attributed to industrial activities in China [35,36]. For this reason, indoor air quality of dwellings in South Korea often becomes worse due to anthropogenic activity including cooking which could lead to increases in the level of PM2.5 as well as NO<sup>2</sup> [37]. In this regard, KME suggested

ventilation for over 30 min three times a day during daytime to improve indoor air quality [33], but this recommendation does not indicate a specific time and lacks tangible evidence of dynamic fluctuation patterns of each pollutant. Therefore, this study aims to analyze the historical records of atmospheric PM10 and PM2.5 concentrations in major cities and regions over the past five years (2015–2019) to assess hourly trends of each pollutant for each season and discuss potential reasons behind the patterns. The findings of this study should provide guidance on the most desirable and undesirable time during the day for not only natural ventilation but also outdoor activities, which could reduce the level of exposure to particulate matter and minimize the adverse effects of air pollution. Moreover, the improved understanding of dynamic trends of PM concentrations and their variations across multiple regions and seasons could help policymakers design a more effective strategy in identifying and monitoring region-specific sources of atmospheric air pollution. –

The structure of this paper is as follows: the next section gives an overview of the eight study regions along with data collection and analytic processes. The following sections show the results of the analyses focusing on hourly trends of outdoor PM concentrations by year, season and region, followed by the discussion and conclusion sections that summarize the findings of this study and the implications and limitations of these findings.

#### **2. Materials and Methods**

The hourly PM10 and PM2.5 concentration data supplied by Air Korea (www.airkorea.or.kr) for KME were obtained for seven metropolitan cities (Seoul, Incheon, Daejeon, Daegu, Ulsan, Gwangju, Busan) and Jeju Island, for the past five years from January 2015 to September 2019. Figure 1 shows the mapped locations of the eight regions in South Korea and Table 1 summarizes the sociodemographic, geographic, meteorological, traffic and other relevant information for each region.

**Figure 1.** Location of the study regions in South Korea.



The observations from multiple monitoring stations in each region were integrated to calculate the average; very few observations were missing due to impediment or equipment failures during the period, and these were excluded from the analysis. Statistical analysis software (SAS for Windows version 9.1; SAS Institute Inc., Cary, NC, USA) was used to calculate the geometric means of PM10 and PM2.5 concentrations (µg/m<sup>3</sup> ) for each hour of a day (24 time slots) in each region, which were then aggregated by year and season to investigate the annual and seasonal changes and variations. Geometric means were chosen as the best summary statistic for the PM data since their distributions were highly skewed. Four distinct seasons in Korea were classified as: spring (March–May), summer (June–August), fall (September–November), and winter (December–February). A series of hourly time graphs were created to illustrate the trends of PM concentrations over a day for each year, season and region, in order to identify good and bad time slots for natural ventilation. Some additional analysis and discussion followed to explain the potential sources of the patterns, including hourly traffic volume in each region.

#### **3. Results**

#### *3.1. Hourly Trends of Outdoor PM Concentrations by Year*

Figure 2 shows the hourly patterns of PM10 and PM2.5 concentrations as an annual average for each year between 2015 and 2019. For both pollutants, there has been a gradual reduction from 2016 until 2018, but they rebounded in 2019. The PM10 concentrations were mostly below the KME's annual average standard of 50 µg/m<sup>3</sup> throughout a single day, while the PM2.5 concentrations were way above the recently-strengthened annual average standard of 15 µg/m<sup>3</sup> during the past five years. Despite some annual fluctuations, the overall hourly patterns for both PM10 and PM2.5 show some similarity throughout the five years; they peaked at 8–11am, declined in the afternoon and rose again in the evening. The second peak in the evening until dawn appears more conspicuous in PM2.5 concentrations, which is due to insufficient air circulation via inversion layer caused by the temperature drop on the earth's surface during these time ranges [38]. tly below the KME's – earth's surface during the

**Figure 2.** *Cont*.

–

–

earth's surface during the

tly below the KME's

**Figure 2.** Diurnal variation of particulate matter (PM) concentrations (µg/m<sup>3</sup> ) by year (all regions).

#### *3.2. Hourly Trends of Outdoor PM Concentrations by Season*

KME's annual average standard of 50

– Figure 3 compares the hourly patterns of both pollutants across the four seasons, showing a noticeable seasonal pattern: high in spring and winter and low in summer and fall. Similar to the previous study [34], PM10 concentration was highest in spring while PM2.5 concentration was highest in winter. Regardless of season, both PM10 and PM2.5 peaked between 9am–noon and decayed afterwards. The slope of decline in the afternoon was much steeper for PM2.5 than PM10, except for summer, when both concentration levels appear relatively stagnant throughout a day. The PM10 concentrations during spring and winter were above the KME's annual average standard of 50 µg/m<sup>3</sup> particularly between 8am and 2pm, while these were under the standard during summer and fall. PM2.5 concentrations were, however, much higher than the annual average standard of 15 µg/m<sup>3</sup> during all four seasons. the KME's annual average standard of 50

**Figure 3.** *Cont*.

–

– –

the KME's annual average standard of 50

– **Figure 3.** Diurnal variation of PM concentrations (µg/m<sup>3</sup> ) by season (2015–2019).

#### *3.3. Hourly Trends of Outdoor PM Concentrations by Season: Regional Variation*

KME's annual average standard of 50 – – Figures 4 and 5 show regional variations of hourly and seasonal trends of PM10 and PM2.5 respectively, including Seoul (capital city of South Korea), six metropolitan cities (Incheon, Daejeon, Daegu, Ulsan, Gwangju, Busan) and Jeju Island. As seen in Figure 4, the seasonal and hourly patterns of PM10 are relatively similar across all regions, exhibiting high concentrations in spring and winter and a peak between 7am and 1pm, but the height of the peak and the declining patterns after the peak seem to vary by region. PM10 concentrations in Busan, Incheon and Ulsan were above the KME's annual average standard of 50 µg/m<sup>3</sup> in the daytime period during spring, but its duration was much longer in Incheon (9am to 6pm) than Busan (9am–1pm) and Ulsan (10–11am). During winter, however, only Incheon and Daegu show PM10 concentration above the standard for a relatively short period (10am to noon). As for PM2.5, illustrated in Figure 5, the regional variation appears more noticeable. The PM2.5 concentrations are relatively similar between spring and winter in Seoul, Busan, Incheon, Gwangju, Ulsan and Jeju Island, while the winter concentrations were substantially larger than the spring concentrations in Daejeon and Daegu. For all eight regions, the spring and winter concentrations of PM2.5 were above the annual average standard of 15 µg/m<sup>3</sup> throughout a given day. However, during the summer and fall seasons, the concentrations were above the standard only during the peak time periods in some regions. Unlike PM10, all regions except for Daegu show two peaks: one in the morning and the other in the evening or night.

The distinct patterns of PM10 and PM2.5 concentrations across the regions is mainly due to different size of particles. Relatively speaking, PM10 is associated with dust on roads, while PM2.5, a very tiny air particulate matter, even smaller than PM10, and relates to emission from vehicles such as aerosol and nitrogen oxide (NOx) [39]. For instance, the hourly trends of PM2.5 concentrations in Ulsan and Jeju Island look quite different from those of other regions, particularly in the afternoon and evening, possibly because of its unique traffic patterns. Ulsan is an industrial city known for factories and plumes leading to male-dominated demographics and workforce characteristics [40]. Thus, considering relatively less traffic congestion during commuting time in the city, air quality might be related more to industrial emissions than local traffic emissions. Jeju Island is a famous tourist destination located at the southern end of South Korea, exhibiting unique demographic, traffic and environmental characteristics. To further investigate the impact of traffic volume on PM2.5 concentration, the traffic volume data were obtained for the same eight regions for three years between 2016 and 2018 and the average amount of traffic was calculates for each hour. Figure 6 shows somewhat

similar hourly patterns to the PM2.5 patterns found in Figure 5, beginning to increase at around 6am, fluctuating during the day, and declining at 6pm. This degree of similarity could indicate the level of contribution of traffic emissions on PM2.5 concentration trends, which should vary across the regions depending on site-specific characteristics [31,41]. As mentioned above, the traffic volumes in Ulsan and Jeju Island are substantially lower and flatter than the other regions due to their unique sociodemographic and mobility patterns.

**Figure 4.** Diurnal variation of PM10 concentrations (µg/m – 3 ) by season in each region (2015–2019).

**Figure 5.** Diurnal variation of PM2.5 concentrations (µg/m<sup>3</sup> ) by season in each region (2015–2019). –

**Figure 6.** Diurnal variation of traffic volume in each region during 2016 and 2018.

#### **4. Discussion**

This study examined hourly and seasonal patterns of PM concentrations for multiple regions in South Korea based on the recent five-year data in order to find desirable and undesirable times for natural ventilation due to high levels of outdoor concentrations. Overall, the average concentrations of both PM10 and PM2.5 peaked in the morning, declined in the afternoon, and rebounded in the evening (particularly for PM2.5). The highest concentrations were generally found in spring (March to May) and winter (December to February), but PM2.5 concentrations were more concerning since they were above the KME's annual average standard of 15 µg/m<sup>3</sup> during peak times in all four seasons. Thus, it is generally advisable that natural ventilation is recommended during the afternoon (2–6pm) but should be avoided in the morning (9am to noon) or in the late evening (8–11pm). However, substantial

differences were observed across the eight regions across the four seasons, suggesting intervention strategies tailored to regional climate and emission patterns. Moreover, the actual decision on whether ventilation is needed at a specific time should be based on a more comprehensive consideration of both indoor and outdoor air quality conditions.

Our results were consistent with previous studies in other countries showing a bimodal pattern with peaks during morning and evening rush hours [42]. This pattern was prominent in many other metropolitan cities such as New York, Los Angeles, Beijing, and London, implying that it was mostly attributed to anthropogenic activity, in particular motor vehicle traffic patterns [43–45]. Likewise, the morning peak noticed in this study could be mostly due to enhanced anthropogenic activity during commuting hours, while the afternoon valley is mainly due to a higher atmospheric mixing layer, which is beneficial for air pollution diffusion [46]. The morning peak was relatively more rigorous than the evening risk because, in South Korea, work starting time is relatively standard (around 8 or 8:30 am) but finishing time is more widely distributed due to variability in work schedule. This could partly explain why the PM levels in the evening were somewhat lower than those in the morning in this study.

The seasonal pattern of PM concentrations in South Korea is partly due to its unique geographic and climatic characteristics. The average PM10 concentration was highest in spring mainly because of yellow dust transported from northern China and the deserts of the Mongolian plateau by the prevailing westerlies [47]. The low PM levels during summer time are attributed to rainout and washout processes due to the rainy period as well as frequent typhoons, and rapid air circulation in the fall helps in reducing PM concentrations [24,48]. Generally, PM2.5 can stay longer in the air compared to PM10 due to a smaller size [49]. Local meteorological conditions (very low temperature, low wind speeds, surface layer inversions) and weak wind circulation during winter make it difficult to remove PM2.5 [50]. An increased usage of heating fuels during the winter raises PM levels, and they stay in the air for a longer period of time due to the cold surface of the earth and a low mixture rate [51,52]. This could explain why PM2.5 concentrations, unlike PM10, were higher in winter than spring since finer particles tend to be generated widely by man-made sources of emission such as solid fuel heating in winter [38].

In general, the level of PMs within a big metropolitan city is affected by traffic volume, point source pollution around the city (e.g., power plant), and weather conditions (especially wind direction) by season [53]. Increased fuel consumption for heating during winter could contribute to increasing the level of PM2.5, but other factors may also apply. The high winter concentrations of PM2.5 in Daejeon and Daegu (shown in Figure 4) are found to be associated with region-specific seasonal factors, such as the seasonal rise of bio-incineration via fuel use for household heating and cooking and emissions from petroleum-related industries [54,55] and traffic emission and gas-form pollutants [34]. The studies showed that the high level of PM2.5 during winter season in Daejeon area could be attributed to emission pollutants from two big coal-fired power plants located in the northwest of the city under the influence of the main wind direction during the winter season [56]. Similar patterns of PM2.5 concentrations over winter in Daegu could also be due to its geographic and climatic conditions. During winter, wind blowing from the northwest causes pollutants from the industrial complex located in the northwestern side of Daegu to enter the downtown area, and the particles remain stagnant in the area due to the lack of air circulation in the basin-shaped city [57]. Of course, a number of other geographic, climatic and socioeconomic factors specific to each region may contribute to PM trends as well [23].

#### **5. Conclusions**

Controlling air pollution is a daunting task in all industrialized countries because so many factors are involved and some of those even cross a country's border. Eliminating or reducing the sources of emission could bring other types of social conflict or dilemma which would require a long-term effort and investment to be resolved. Natural ventilation could be considered as a simple and low-cost solution to diminish the level of exposure to indoor air pollutants, but only when the outdoor air quality is acceptable. This study emphasizes the importance of thorough monitoring of ambient air quality to be used in providing guidelines for when indoor environments should be ventilated in different region. Despite substantial distinctions across the regions by season, this study provides some general suggestions that the time between 2–6pm is most suitable for natural ventilation but it is not desirable either in the morning or late evening, particularly during spring and winter. It also highlights PM2.5 concentrations surpassing the KME's annual average standard of 15 µg/m<sup>3</sup> during peak times in all four seasons and rebounding in size since 2019. Although further studies would be needed to confirm this suggestion and implement solutions in practice, the results of this study can be used as basic information for designing a comprehensive environmental health policy in consideration of dynamic exposure to air pollution.

Of course, it would be ideal to consider simultaneously both indoor and outdoor PM values in order to confirm the suggestions of this study. However, we believe that this study is valid in itself because indoor PM data are not readily available in most of the households and, even if available, the patterns would vary greatly vary according to building structure and lifestyle. This study aims to provide a broad guideline on appropriate times for natural ventilation based on outdoor air quality, particularly where indoor PM levels are unknown, but specific implementation should be carried out by considering indoor-outdoor dynamics pertaining to each building or household. Future study should be directed towards the confirmation of the trends found in this study based on actual experimental data on indoor and outdoor PM concentrations measured in sample buildings in multiple regions, instead of the public aggregated data. In addition, gathering more data on covariates indicating geographic, sociodemographic, climatic, industrial and traffic patterns in each region would enable multivariate analysis and modeling for the determinants of spatiotemporal and seasonal changes in PM concentrations. In spite of limitations and the inability to explain some patterns, this study can inform the public regarding the desirable or undesirable times for natural ventilation as well as outdoor activities in each season and region. Moreover, it can encourage policymakers to design season-specific environmental management strategies and practices tailored to regional climate and emission factors such as regional traffic monitoring networks and air quality alert systems for ventilation, which may promote a cost saving, by specialization of intervention, and enhance policy outcomes.

**Author Contributions:** Conceptualization, D.K. and S.S.; Data curation, H.-E.C.; Formal analysis, H.-E.C. and S.S.; Funding acquisition, S.S.; Methodology, D.K. and S.S.; Project administration, S.S. and W.-M.G.; Supervision, S.S. and W.-M.G.; Visualization, H.-E.C.; Writing—original draft, H.-E.C. and S.S.; Writing—review & editing, D.K. All authors read and approved the final manuscript.

**Funding:** This work was supported by Korea Environment Industry & Technology Institute (KEITI) through the Environmental Health Action Program, funded by Korea Ministry of Environment (MOE) (grant number: 2018001350005). This support is greatly appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Evolution of Urban Haze in Greater Bangkok and Association with Local Meteorological and Synoptic Characteristics during Two Recent Haze Episodes**

**Nishit Aman 1,2 , Kasemsan Manomaiphiboon 1,2,\*, Natchanok Pala-En <sup>3</sup> , Eakkachai Kokkaew 1,2, Tassana Boonyoo <sup>4</sup> , Suchart Pattaramunikul <sup>4</sup> , Bikash Devkota 1,2 and Chakrit Chotomosak <sup>5</sup>**


Received: 18 November 2020; Accepted: 15 December 2020; Published: 18 December 2020 -

**Abstract:** This present work investigates several local and synoptic meteorological aspects associated with two wintertime haze episodes in Greater Bangkok using observational data, covering synoptic patterns evolution, day-to-day and diurnal variation, dynamic stability, temperature inversion, and back-trajectories. The episodes include an elevated haze event of 16 days (14–29 January 2015) for the first episode and 8 days (19–26 December 2017) for the second episode, together with some days before and after the haze event. Daily PM2.5 was found to be 50 µg m−<sup>3</sup> or higher over most of the days during both haze events. These haze events commonly have cold surges as the background synoptic feature to initiate or trigger haze evolution. A cold surge reached the study area before the start of each haze event, causing temperature and relative humidity to drop abruptly initially but then gradually increased as the cold surge weakened or dissipated. Wind speed was relatively high when the cold surge was active. Global radiation was generally modulated by cloud cover, which turns relatively high during each haze event because cold surge induces less cloud. Daytime dynamic stability was generally unstable along the course of each haze event, except being stable at the ending of the second haze event due to a tropical depression. In each haze event, low-level temperature inversion existed, with multiple layers seen in the beginning, effectively suppressing atmospheric dilution. Large-scale subsidence inversion aloft was also persistently present. In both episodes, PM2.5 showed stronger diurnality during the time of elevated haze, as compared to the pre- and post-haze periods. During the first episode, an apparent contrast of PM2.5 diurnality was seen between the first and second parts of the haze event with relatively low afternoon PM2.5 over its first part, but relatively high afternoon PM2.5 over its second part, possibly due to the role of secondary aerosols. PM2.5/PM<sup>10</sup> ratio was relatively lower in the first episode because of more impact of biomass burning, which was in general agreement with back-trajectories and active fire hotspots. The second haze event, with little biomass burning in the region, was likely to be caused mainly by local anthropogenic emissions. These findings suggest a need for haze-related policymaking with an integrated approach that accounts for all important emission sectors for both particulate and gaseous precursors of secondary aerosols. Given that cold surges induce an abrupt change in local

meteorology, the time window to apply control measures for haze is limited, emphasizing the need for readiness in mitigation responses and early public warning.

**Keywords:** urban haze; temperature inversion; Obukhov length; HYSPLIT; biomass burning; cold surge, emission

#### **1. Introduction**

Particulate matter (PM) with a size less than or equal to 2.5 µm (PM2.5) is an environmental concern worldwide. Suspension of such particulate matter in the atmosphere (known as aerosols) affects human health, atmospheric visibility and also impacts weather and climate both directly and indirectly [1]. The aerosols are either emitted directly in the atmosphere (known as primary sources) from combustion, wind-borne dust, sea spray, volcanic emission, and biogenic aerosol or formed in the atmosphere by conversion of primary precursor gases to secondary particles through nucleation and complex multiphase chemical reactions. High PM pollution in any place always depends upon the complex interplay between local emissions, secondary particle formation, along with local and synoptic meteorology [1,2]. Around 58% of the world population lives in areas with PM2.5 > 35 µg m−<sup>3</sup> (in terms of daily average according to the WHO Interim Target 1) [3], many of which are highly urbanized areas or large cities. Although PM pollution episodes (sometimes known as haze episodes) are limited to few days to weeks, exposure to high PM levels even for shorter times can have an effect on human health and ecosystems [3].

Haze pollution has been studied globally to understand its formation and evolution mechanism [4–7], potentials sources contribution [5,8,9], mitigation strategies [10,11], and early warning or forecasting [12]. A haze episode may be induced by one or a number of factors combined, which encompasses emissions, secondary aerosols [13–17], and atmospheric transport [8,17], with unfavorable weather conditions acting as an accelerating factor [18–25]. The effects of the atmospheric boundary layer (ABL) structure, near-surface atmospheric stability, and synoptic conditions have received the attention given that they strongly dictate how haze and its associated thermal and dynamical processes evolve with time in the lower part of the troposphere [23–25]. Quan et al. (2013) [19] and Petaja et al. (2016) [20] suggested a positive feedback cycle for heavy air pollution where heat flux decrease significantly due to decreased solar radiation blocked by the haze layer, which in turn further decreases ABL height and trap air pollutant within. Tie et al. (2017) [21], and Liu et al. (2018) [22] reported that decreased ABL height leads to increased relative humidity, which enhances secondary aerosol formation and also promotes the hygroscopic growth of aerosols and light scattering. The latter effect reduces incident solar radiation more. Temperature inversion at different levels of the troposphere may be caused by radiation, advection, large-scale subsidence, etc. Low-level temperature inversion layers effectively reduce the atmospheric volume available for diluting air pollutants while weak winds poorly ventilate them out of a polluted area [23–25]. Turbulence is another important factor as the capability of mixing airborne constituents vertically in the ABL. Low turbulence, typically under stable conditions, suppresses vertical mixing, which allows pollutants and/or precursors to accumulate at near-surface levels [23,26]. Synoptic weather is another relevant aspect because it sets favorable or unfavorable background conditions to haze formation and evolution and influences local or urban-scale weather processes as well [23–27].

Bangkok, the capital of Thailand, and its five neighboring provinces (Samut Prakan, Samut Sakhon, Nonthaburi, Nakhon Pathom, and Pathum Thani) are collectively known as Greater Bangkok (GBK) or Bangkok Metropolitan. It is one of the largest urban agglomerations in Southeast Asia. It has experienced haze pollution typically found in the dry season, posing great concern to the general public and challenges to the local and central governments for mitigation and prevention. PM2.5 exceeds the daily (i.e., 24-h average) national ambient air quality standard (NAAQS) of 50 µg m−<sup>3</sup> several times

per year, according to the Pollution Control Department (PCD) [28], particularly in the dry season. Several PM-related studies in GBK have been conducted e.g., [29–38], see Table S1 in Supplementary Materials for details, ranging from source apportionment, chemical characterization, emission inventory, and human health. On-road vehicles (i.e., traffics) and biomass burning were identified as the major PM2.5 sources [29–32]. In the dry season, agricultural burning to clear crop residues on lands within GBK and its vicinity and forest fires in the northern region are generally intensified, and smoke can also be dispersed or transported to GBK [31–34]. However, biomass burning contributes little to air pollution in GBK during the summer due to a shift in the prevailing winds [35]. Pham et al. (2008) [36] estimated gaseous and particulate emissions from industrial and power plant sectors, finding the central and central regions (among all regions in Thailand) to have the largest intensities and shares. Secondary aerosols were also reported as a nonnegligible contributor [30,32]. Effects of increased PM2.5 on health have been emphasized and quantified [37,38]. Although these past studies provide useful information of the PM sources and PM effects for the study area in question, elevated haze and its associated meteorological dependence are still lacking or little addressed. Accordingly, this study aims to fill this knowledge gap by investigating how urban haze is influenced by meteorology at both local and synoptic scales. Specifically, we seek to understand how events of elevated haze (in terms of PM2.5) temporally evolve with meteorological factors using several observational datasets combined, with particular attention to cold surge as a synoptic disturbance relevant to the region and here suspected to induce favorable conditions for haze to elevate. Here, an intensive observational analysis was performed and discussed for two recent wintertime haze episodes in GBK using data from various sources to examine the association of haze with local and synoptic weather conditions.

#### **2. Data and Methods**

#### *2.1. Study Area*

Greater Bangkok (GBK) is located in the lower part of Central Thailand. It currently has a registered population of 11 million and an area of 7762 km<sup>2</sup> [39]. It is the largest national hub of the economy, accounting for 46% of its total gross domestic product (GDP) [40]. It has a complex mixed (built and natural) landscape with the co-existence of commercial, residential, agricultural, and industrial areas [41]. The overall terrain of GBK is generally flat with limited heights (<10 m above mean sea level or MSL). Its general climate is tropical and humid and governed mainly by the northeast monsoon (November–February as the winter) and the southwest monsoon (May–October as the wet season) [42]. The former monsoon brings cool, dry air from continental mid-latitudes over which persistent strong high-pressure systems are present. The latter monsoon brings moist air from the Indian Ocean and the Gulf of Thailand, causing abundant rain in most parts of Thailand. The monsoon trough or intertropical convergence zone (ITCZ), which moves along the north-south direction, and tropical cyclones developed in the North Indian and Western Pacific Ocean Basins can modulate rain at a sub-seasonal scale. Importantly, a cold surge is a synoptic phenomenon, characterized by a transient southward propagation or extension of a high-pressure system from mid-latitudes (i.e., mainland China and East Asia) to the Indochina Peninsula and the equatorial South China Seas [43–45]. During the winter, cold surges occur episodically more often, strengthening the northeast monsoon. The arrival of a cold surge typically brings strong winds due to high-pressure gradients with an abrupt drop in temperature. Once it weakens, the cold surge recedes back or dissipates. The transitional period between the two monsoons (March–April) has relatively warm conditions, corresponding to the summer season. Here, the winter and summer combined are called the dry season [42].

#### *2.2. Data*

The PCD is the main government agency that administers air quality monitoring stations across Thailand. Here, hourly PM2.5 and PM<sup>10</sup> data for the years 2015–2017 at three air quality stations (P27, P59, and P61) in the study area were requested and obtained (Table 1 and Figure 1c). P27 is located near a busy major highway (by about 70 m). P59 is in a semi-general area but not far from a busy local street and a major expressway to the west by 0.6 km. P61 is in a general area with no major road nearby within 1 km, representatively selected as the key station to support several parts of the analysis. It is noted that the PM2.5 data at P59 are available only from April 2015. Given missing data and discontinuity in meteorological measurements at these three stations, the 100-m tower (M6) of the PCD at Techno Thani in Pathum Thani was the main source of meteorological data for use. The tower measures air temperature (*T*) at 2 m, 50 m, 75 m, and 100 m, wind speed (*WS*) and wind direction (*WD*) at 10 m, 50 m, and 100 m, and other near-surface variables, i.e., relative humidity (*RH*), rain (*RN*) and global radiation (*GR*). PM2.5 and PM<sup>10</sup> are real-time monitored using the standard beta-ray method at 3 m above ground level (AGL). Data are intensively quality controlled/assured by the PCD internally before distribution. The detectable limits or probable ranges of data are PM2.5 and PM<sup>10</sup> (3 to 1000 µg m−<sup>3</sup> ), temperature (–5 to 50 ◦C), relative humidity (0 to 100%), wind speed (0 to 50 ms −1 ), wind direction (0 ◦ to 360 ◦ ), rain (0 to 1000 mm h −1 ) and global radiation (0 to 1000 W m−<sup>2</sup> ) [46,47]. We applied these ranges for data screening and also removed suspicious or erratic values if found visually. Relative humidity data at the M6 tower are absent throughout 2017. Simple linear-regression extrapolation using the data at P59 was employed to gap-fill them (details not shown). The tower is in a suburban/rural area, with most of the land near or surrounding the tower being water and paddy fields. Built-up areas are also well present over its northeastern quadrant within 2 km. Since many buildings taller than 10 m are also in their proximity of 100 m, 10-m wind data were discarded. Upper-air sounding data at the Bang Na weather station of the Thai Meteorological Department (TMD) was obtained (available at http://weather.uwyo.edu/upperair/sounding.html). However, radiosonde soundings are routinely operated only at 7 local time (LT) (i.e., 0 UTC), with pilot balloons alone at 1, 13, and 19 LT. The upper-air sounding and data of the TMD typically follow the standard quality assurance of the World Meteorological Organization (WMO). Here, radiosonde-based upper-air data up to a height of 4 km were considered and extracted. To convert hourly data to a daily scale for air quality and meteorological variables, 24-h (1–24 LT) averaging was typically applied. For global radiation, a period of 11–16 LT was used to represent late-morning to mid-afternoon hours. All computations and statistical tests were performed using software R (R development core team, 2019) [48]. In any statistical calculation, 50% of data as valid/ non-missing was necessarily required as the minimum threshold. <sup>−</sup> – − − − – –

**Figure 1.** (**a**) Thailand; (**b**) Greater Bangkok and its provinces; (**c**) monitoring stations. In (**b**), the urban built-up areas (urban residential, industrial, and commercial classes combined) are shown in shaded

gray, based on LDD (Land Development Department) (2016) [41]. In (**c**) P27, P59, and P61 are the air quality stations at Samut Sakhon Witthayalai School, Public Relation Department, and Bodindecha School, respectively, M6 is the 100-m tower at Techno Thani, and the Thai Meteorological Department (TMD) is the standard weather station (WMO No. 48453) at Bang Na District.


**Table 1.** Monitoring stations considered in this study.

Note: PM2.5 and PM<sup>10</sup> (µg m−<sup>3</sup> ): particulate matter with size not larger than 2.5 µm and 10 µm, respectively; *T2* and *T50* ( ◦C): temperature at 2 m and 50 m, respectively; *WS50* and *WS100* (m s −1 ): wind speed at 50 m and 100 m, respectively; *WD50* and *WD100* (degrees from the north): wind direction at 50 m and 100 m, respectively; *RH* (%): relative humidity; *RN* (mm): rain; *GR* (W m−<sup>2</sup> ): global radiation; the stations have the terrain elevations of 2–4 m MSL. − − −

–

#### *2.3. Selection of Haze Episodes*

A haze day is here defined as the day with PM2.5 exceedance (i.e., daily PM2.5 level ≥ 50 µg m−<sup>3</sup> ) registered at one station at least. Then, a haze event is simply defined as the period of consecutive haze days. As seen from Figure 2a (also see Figure S1 in Supplementary Materials), December–February is typically the period when haze intensifies. Given emphasis to extended (e.g., haze days over a week or longer) severe episodes and the amount of the valid surface and upper-air data, two haze episodes (EP1 and EP2) were representatively chosen (Figure 2b,c). In EP1, it comprises the following three periods in sequence: pre-HZ1, HZ1, and post-HZ1. The haze event (namely, HZ1) spans 16 days (14–29 January 2015). Pre-HZ1 and post-HZ1 are the periods before (7 days) and after (4 days) HZ1, respectively, when PM2.5 was at relatively low levels. Similarly, EP2 covers pre-HZ2, HZ2, and post-HZ2. The haze event (namely, HZ2) spans 8 days (19–26 December 2017), with 6-day pre-HZ2 and 5-day post-HZ2. ≥ − – – –

**Figure 2.** *Cont*.

– – **Figure 2.** (**a**) Monthly particulate matter (PM) with a size less than or equal to 2.5 µm (PM2.5 ) and number of haze days observed at P61 over 2015–2017; (**b**) daily PM2.5 during haze episode 1 (EP1); (**c**) daily PM2.5 during haze episode 2 (EP2). In (**a**), the filled and unfilled circles are the averages, the vertical bars are the standard deviations, and the *x*-axis labels (J, F, M, . . . , N, and D) denote the months of the year. In (**b**,**c**), the suffixes N, D, J, and F on the *x*-axis labels correspond to November–February, respectively. In (**b**,**c**), the dark-gray shading marks the HZ periods, and the light-gray shading marks the pre-haze event (HZ) and post-HZ periods for each haze episode.

#### *2.4. Synoptic Patterns*

The 6-hourly 0.5 ◦ -resolution Climate Forecast System Version 2 (CFSv2) reanalysis data [49] were used to provide sea-level pressure at 13 LT for all days in EP1 and EP2 to help construct daily synoptic surface charts. Moreover, satellite infrared images (at the same local time) captured by geostationary Himawari-8 (available at http://weather.is.kochi-u.ac.jp/sat/gms.sea/) were used to examine the presence of low and high clouds over the central region of Thailand, the ITCZ, and tropical cyclones. It is noted that the constructed CFSv2-based surface charts were also compared with the corresponding TMD surface charts (available at https://www.tmd.go.th/en/weather\_map.php), both of which were found to be similar. A total of four simple synoptic patterns (numbered as 0, 1, 2, and 3) were proposed here to help support the analysis, using our visual examination of the constructed charts. The representative surface charts for EP1 and EP2 are shown in Figures S2 and S3 (Supplementary Materials), respectively. The four synoptic patterns classified are as follows:


–

#### *2.5. Temperature Inversion and Obukhov Length*

A temperature inversion refers to an atmospheric condition when the air temperature increases with height, and its presence can restrict the volumetric dilution of air pollutants [2,50]. Using the radiosonde-based upper-air data (available only at 7 LT), temperature inversion layers were identified. Those only found over the heights of 100–1500 m are of interest here. Those below 100 m were excluded since they are typically induced by continuous radiative surface cooling over the nighttime and early morning hours. The upper limit of 1500 m conservatively marks the ABL thickness. A single inversion layer is given as all successive vertical levels from the sounding with temperature monotonically increasing with height. Inversion intensity (*IV*) is a simple parameter used to indicate the extent of difficulty or blockage to which air pollutants penetrate an inversion layer, which was here computed

as the rate of temperature change (◦C per 100 m depth) from the bottom to the top of the inversion layer, similar to Dai et al. (2020) [7]. Only inversion layers with *IV* > 0.1 ◦C per 100 m were considered. Obukhov length (*L*) is an important measure of near-surface dynamic stability [2,18,23], suggesting the capability of vertical mixing for air pollutants between the surface and higher levels. Its positive/negative values correspond to stable/unstable conditions. The smaller magnitude of *L*, the larger degree of stability/instability. A very large *L* in magnitude corresponds to the neutral condition. The term "dynamic" implies that both mechanical and thermal turbulence production processes are taken into account. Following the Monin–Obukhov similarity theory [18,23], *L* is computed by

$$L = \frac{1}{2} \frac{(T\_{v1} + T\_{v2})u\_\*^2}{kg\theta\_\*} \,\text{.}\tag{1}$$

$$\mu\_{\ast} = k \, \mathcal{U}(z\_{\mathcal{U}}) \Big[ \ln \left( \frac{z\_{\mathcal{U}}}{z\_0} \right) - \Psi\_M \left( \frac{z\_{\mathcal{U}}}{L} \right) + \Psi\_M \left( \frac{z\_0}{L} \right) \Big]^{-1}, \text{ and} \tag{2}$$

$$\theta\_\* = k \left( \theta\_{v2} - \theta\_{v1} \right) \left[ \ln \left( \frac{z\_{\theta\_{v2}}}{z\_{\theta\_{v1}}} \right) - \Psi\_H \left( \frac{z\_{\theta\_{v2}}}{L} \right) + \Psi\_H \left( \frac{z\_{\theta\_{v1}}}{L} \right) \right]^{-1} \tag{3}$$

where *u*∗ is the frictional velocity, θ<sup>∗</sup> is the temperature scale, *k* is the von Karman constant (0.4), *g* is the acceleration due to gravity, *z<sup>0</sup>* is the roughness length, *z<sup>u</sup>* is the single measurement height of wind (here, 50 m or 100 m separately), *U*(zu) is the hourly wind speed, *z*θ*v*<sup>1</sup> and *z*θ*v*<sup>2</sup> are the two measurement heights of temperature (here, 2 m and 50 m, respectively), θ*v1* and θ*v2* are the virtual potential temperature at those two heights, *Tv1* and *Tv2* are the virtual temperatures at those two heights, respectively, and Ψ*<sup>M</sup>* and Ψ*<sup>H</sup>* are the Businger stability correction functions for wind and temperature, respectively. As in Kamma et al. (2020) [51], the surrounding area of the M6 tower within a 2 km radius was assessed using Google Earth (https://www.google.com/earth/) (Google, Mountain View, CA, USA) and visually examined, finding ponds and paddy fields being dominant with built-up areas present in its northeast quadrant. Using the classification by Stewart and Oke (2012) [52], the approximate local climate zone is "sparsely built" with terrain roughness (or Davenport) class 5, whose roughness length (*z0*) equals 0.25. The concept of virtual (potential) temperature is necessary for humidity correction, which was implemented using the "aiRthermo" package in R [53]. In doing so, surface pressure data were required, extracted from 3-hourly 0.25◦ -resolution global land data assimilation system (GLDAS) [54] (available at https://giovanni.gsfc.nasa.gov/giovanni/), and then linearly interpolated to an hourly scale. If pressure at any higher levels was needed, the hydrostatic adjustment was applied. We attempted to compute hourly *L* with *z<sup>u</sup>* = 50 m and *z<sup>u</sup>* = 100 m separately and found a very high correlation (0.99) between the results from the two cases. Thus, the average *L* values over both cases were used.

#### *2.6. Back-Trajectories*

To investigate the potential transport of air pollutants from nearby and far areas [23,24], daily kinematic back-trajectories were simulated for all days in both episodes (EP1 and EP2) by the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT) of the National Oceanic and Atmospheric Administration (NOAA) [55]. Each back-trajectory starts at 13 LT (as a typical midday time with developed ABL) and at 500 m AGL (as a typical mid-ABL height) and migrates backward in time for 48 h. HYSPLIT was run online (at https://www.ready.noaa.gov/HYSPLIT.php) using hourly 0.5◦ -resolution global data assimilation system (GDAS) data for driving wind fields. Given that biomass burning (agricultural burning and forest fires) in Upper Southeast Asia is well present in the dry season [32–34], daily 1-km active fire hotspots detected by the MODIS (Moderate Resolution Imaging Spectroradiometer) sensors onboard of both Terra and Aqua satellites (MCD14ML Collection 6) [56] (available at https://firms.modaps.eosdis.nasa.gov/download/) were downloaded for both episodes. The fire hotspots were then summed and gridded to 0.5◦ according to the pre-HZ, HZ, and post-HZ periods for each episode.

#### **3. Results and Discussion**

The urban haze in GBK is generally associated with multiple factors and their interdependence or interplays, ranging from emissions from local sources and biomass burning, mid- and long-range transport, secondary aerosols, and meteorological conditions at both local and synoptic scales. Here, the last factor is our main focus, for which general and distinct meteorological features and how they are coupled with the urban haze evolving during the selected two haze episodes are described.

#### *3.1. Haze Episode EP1*

−

This episode (EP1) occurred mostly in January 2015, and the haze event (HZ1) spans14–29 January. Figure 3 displays the day-to-day variation of PM2.5, PM2.5/PM<sup>10</sup> ratio, temperature, relative humidity, wind speed, global radiation, rain, and synoptic pattern in the episode. Every variable was of 24-h average, except for global radiation (11–16 LT). For PM2.5/PM10, its daily values were of the 24-h average of the ratio of hourly PM2.5 to hourly PM10. As seen from the figure, PM2.5 was relatively high in HZ1 but low during both pre-HZ1 and post-HZ1. It showed two peaks, 80.7 µg m−<sup>3</sup> on 23 January and 68.3 µg m−<sup>3</sup> on 27 January. PM2.5/PM<sup>10</sup> did not appear to vary much (ranging between 0.45 and 0.66), being slightly higher in HZ1 than pre-HZ1 and post-HZ1, suggesting fine PM mode to increase when the haze was more developed. Temperature and relative humidity were 28.2 ◦C and 62.9%, at the start of EP1, respectively. Both decreased to 22.0 ◦C and 47.6%, respectively, at the start of HZ1, as caused by the synoptic change from Pattern 1 to Pattern 2 (i.e., cold surge reaching GBK with cool, dry air), but later climbed up continuously until HZ1 ended, corresponding synoptically to the weakening cold surge or its eventual dissipation. Winds appeared to follow the synoptic patterns, which were relatively strong during the cold surge arrival (i.e., Pattern 2), as seen on 11–14 January and became weaker on the other days, particularly during HZ1, supporting the buildup of haze.

**Figure 3.** Daily PM2.5, PM2.5/PM<sup>10</sup> (at station P61) and other meteorological variables during EP1.

−

−

Global radiation was relatively low in pre-HZ1 (with the minimum of 237.1 W m−<sup>2</sup> ) but turns relatively large in HZ1 (with the maximum of 849.4 W m−<sup>2</sup> ), which could be attributed partly to fewer clouds induced by the high-pressure cold surge (i.e., shifted synoptic patterns). No rain was observed over the entire episode, except on one single day in pre-HZ1 with a light amount (4.4 mm).

In view of diurnal variation (Figure S4 in Supplementary Materials), during pre-HZ1, PM2.5 was generally low (<30 µg m−<sup>3</sup> for most of the hours). Once the haze sets in, PM2.5 showed diurnality more clearly, i.e., relatively low in the afternoon but high in the nighttime and morning. Furthermore, we noticed that diurnal variation in the first (14–23 January) and second (24–29 January) parts of HZ1 showed a sharp contrast and PM2.5 became relatively high in the afternoon in second part as compared to the afternoon PM2.5 in the first part. We looked to daytime Obukhov length but did not find any dramatic change in instability over the days at all (Figure S4 in Supplementary Materials). Hence, it was not possible to explain it directly, and we then suspected that secondary aerosols were enriched in the afternoon for the later part of HZ1, given that this haze event extended as long as 16 days with increasing temperature and humidity and ample global radiation.

Figure 4 shows daily vertical temperature profiles with inversion layers identified. Low-level inversion occurred for 4 days (out of 7) during pre-HZ1 but persistently appeared almost every day (14 out of 15 days with non-missing data) in HZ1. The latter highlights the limited diluting volume for air pollutants and facilitate PM2.5 elevation. Inversion intensity varied day-to-day during HZ1, with a minimum of 0.13 ◦C/100 m on 22 January and a maximum of 2.39 ◦C/100 m on 25 January. As mentioned previously, the temperature was relatively low after the cold surge arrives, and before it receded or dissipated, the low-level inversion was more easily induced. Even though global radiation and near-surface dynamic instability were present, it still took more heat or longer time to warm the surface to break up the inversion, based on the concept of the bulk model of daytime mixing height [47]. With a sequence of inversion-breakup failures, multiple low-level inversion layers could form, which was seen twice in HZ1. Atmospheric aerosols were known for radiative effects, e.g., black carbon or soot to absorb heat and sulfate to scatter radiation, and they could play a role in modifying temperature profiles. However, this subject is beyond the scope of the current study. Upper-level inversion also existed before the cold surge arrival and maintained over most of EP1, suggesting the presence of large-scale subsidence inversion aloft.

Lastly, the transport of air pollutants was investigated using back-trajectories (Figure 5 and Figure S5 in Supplementary Materials). Over the course of EP1, most of the back-trajectories moved from the eastern and northeastern directions, allowing traveling air masses to absorb and carry air pollutants or emissions to the study area. Coincidentally, biomass burning was present in Laos, Cambodia and the central and northeastern regions of Thailand and became relatively intense during HZ1. Only two air masses (out of 7) pass through fire areas in pre-HZ1, while all air masses move over such areas during HZ1 (Figure S6 in Supplementary Materials). Slow-moving and low-level back-trajectories found in the second part of HZ1 (24–29 January) possibly worsened the PM2.5 situation in the study area (Figure S5 in Supplementary Materials). During post-HZ1, half of the back-trajectories (2 out of 4) were maritime (i.e., originating in or passing over the Gulf of Thailand) and thus relatively clean. It is, thus, fair to say that, besides local and synoptic meteorology, the long-range transport potentially impacted the urban haze in GBK to a certain extent in this episode.

– – **Figure 4.** Daily vertical temperature profiles at 7 local time (LT) during EP1. The blue rectangles are the low-level inversion layers (found at 0.1–1.5 km), while the gray rectangles are the upper-level ones (found at 1.5–4 km). In the figure, *IV* represents the inversion intensity in ◦C/100 m.

68

**Figure 5.** Daily 48-h back-trajectories: (**a**) EP1 in 2015; (**b**) EP2 in 2017.

#### *3.2. Haze Episode EP2*

− − − − This episode (EP2) occurred in December 2017, with pre-HZ2 on 13–18 December, the haze event (HZ2) on 19–26 December 2017 and post-HZ2 on 27–31 December. Daily PM2.5, PM2.5/PM<sup>10</sup> ratio, temperature, relative humidity, wind speed, global radiation, rain, and synoptic pattern for EP2 are given in Figure 6. As seen, the average PM2.5 was less than 50 µg m−<sup>3</sup> for all days during pre-HZ2. At the start of HZ2 on 19 December, PM2.5 increased to 54.4 µg m−<sup>3</sup> and remains consistently higher than 50 µg m−<sup>3</sup> until 26 December (except on 25 December). PM2.5 dropped sharply after 26 December due likely to wet scavenging by rain reported for three consecutive days (26–28 December) caused by a tropical depression in the Gulf of Thailand [57]. The tropical depression was the final phase of Typhoon Tembin maturely developed in the South China Sea. PM2.5/PM<sup>10</sup> appeared to be higher in EP2 (0.69–0.94) than that in EP1 (0.45–0.66). The lower ratio in EP1 was due partly to significant contributions from agricultural burning and forest fires. Based on a global emission database by Klimont et al. (2017) [58], it was found that the PM2.5/PM<sup>10</sup> ratio for emissions from agricultural burning and forest fires was lower as compared to those for other anthropogenic emissions. Thus, biomass burning observed in EP1 may have decreased the PM2.5/PM<sup>10</sup> ratio. Haze episode EP2 did not show effects of biomass burning, which is explained in the latter part of this section. Temperature and relative humidity on 13 December was 30 ◦C and 74.1%, respectively and did not change much until 16 December. The synoptic change from Pattern 1 to Pattern 2 on 16–17 December, which caused an abrupt decrease in both variables until Pattern 2 prevailed, i.e., 20 December. Temperature and relative humidity dropped to 20.5 ◦C and 45.1%, respectively, on 20 December and gradually climbed up to 28.8 ◦C and 56.4%, respectively, on 24 December due to the weakening of cold surge leading to the return of the warm moist air. Once the cold surge recedes, Pattern 0 and Pattern 1 were dominant for most of the days. Global radiation increased from 393.8 W m−<sup>2</sup> on 13 December to 594.7 W m−<sup>2</sup> on 19 December (i.e., the start of HZ2) due to less cloud cover caused by the cold surge. Over 24–27 December, both temperature and global radiation declined and relative humidity increased, as influenced by the tropical depression. The tropical depression and rainfall observed during 26–28 December likely increased relative humidity during 24–27 December. The temperature

rose after 27 December until the end of EP2 due to an increase in global radiation, which in turn decreased relative humidity. Global radiation was consistently higher during the haze event until 24 December with less variation and decreased later between 24 and 27 December and then again increased until the end of EP2 during 27–31 December. The wind followed the synoptic patterns quite closely. It was relatively strong during the active phase of cold surge during pre-HZ2 (particularly, 19–22 December) and starting four days of HZ2 while for the other days it was relatively weak.

− −

**Figure 6.** Daily PM2.5, PM2.5/PM10, and other meteorological variables during EP2.

The diurnality during HZ2 was clearer as compared to that during the pre-HZ2 and post-HZ2 (Figure S4 in Supplementary Materials). An obvious diurnal pattern with higher PM2.5 during the night and early morning and drop during the afternoon were observed during HZ2. Obukhov length indicated stability conditions for some hours during a few days. However, no clear change in stability was observed throughout the haze event before the tropical depression.

The vertical temperature profiles with inversion layers for all days in this episode are shown in Figure 7. Three out of six days during pre-HZ2 do not show any temperature inversion while all non-missing seven days during HZ2 between 19 and 26 December show temperature inversions of varying intensity strength with a minimum of 0.25 ◦C/100 m on 26 December and a maximum of 2.0 ◦C/100 m on 24 December (Figure 7). Multiple low-level inversion layers were found during 21–22 December. 2 out of 5 days during post-HZ2 did not show any inversion. High-level inversions existed for most of the days during EP2 due to the presence of large-scale subsidence inversion aloft in the winter.

**Figure 7.** Same as Figure 4, but for haze episode EP2. *IV* represents the inversion intensity in ◦C/100 m.

– Lastly, we used back-trajectories from HYSPLIT to investigate the effect of long-range transport of air pollutants (Figure 5 and Figure S5 in Supplementary Materials). Except for 13–14 December, all the days had longer trajectories and coming from the eastern and northeastern directions. Trajectories during pre-HZ2 and post-HZ2 were mostly low level, while for HZ2, it was at a high level for many days (Figure S5 in Supplementary Materials). Fewer fire hot spots for all three periods: pre-HZ2, HZ2, and post-HZ2 were found in nearby regions suggesting that this haze episode HZ2 was mainly caused by the synergetic effect of local emission and local/synoptic meteorology (Figure S6 in Supplementary Materials). Comparing the two episodes, the EP2 back-trajectories tended to be longer (i.e., move faster). Fast-moving back-trajectories generally had short residential times to absorb atmospheric constituents. These imply limited contribution from biomass burning to haze, and elevated haze in this episode may have been more driven by local emissions.

It is noted that the contrasting results of the potential impact of biomass burning on the two episodes suggest that the influence of biomass burning may have varied by episode. Following our literature survey for Bangkok or Greater Bangkok (see Table S1 in Supplementary Materials), Phairuang et al. (2019) [33] found agricultural burning in the central region of Thailand and forest fires in remote areas (e.g., the northern region) to contribute to particulate matter in Bangkok, but the impact is more evident during the winter (November–January) when air masses are generally continental and

–

favorably transporting haze to the study area, as compared to the summer (March) when air masses are relatively clean due to originating in or passing over the Gulf of Thailand. Similarly, Kayee et al. (2020) [35] looked into a summertime haze episode and found no significant contribution from biomass burning in the northern region due likely to lack of continental air masses. Dejchanchaiwong et al. (2020) [34] reported biomass burning in Thailand and Cambodia with favorable air masses as a major contributor to a wintertime haze episode.

#### **4. Conclusions**

Several local and synoptic meteorological aspects associated with the selected two wintertime haze episodes in Greater Bangkok were examined using observational data, including classified four synoptic patterns, day-to-day and diurnal variation, dynamic stability, temperature inversion, and back-trajectories. The first episode included an elevated haze event of 16 days (14–29 January 2015), together with some days before and after the haze event. However, the second episode had an elevated haze event of only 8 days (19–26 December 2017). Daily PM2.5 was found to be 50 µg m−<sup>3</sup> or higher over most of each haze event. The two haze events commonly had cold surges as the background synoptic feature to initiate or trigger haze evolution. The cold surge reached the study area before the start of each haze event, causing temperature and relative humidity to drop abruptly initially but then to gradually increase as the cold surge weakened or dissipated. Wind speed was relatively large when the cold surge was active. Global radiation was generally modulated by cloud cover, which turns relatively high during each haze event because the cold surge induced fewer clouds. Daytime dynamic stability was generally unstable along the course of each haze event, except being stable at the ending of the second haze event due to a tropical depression. In each haze event, low-level temperature inversion existed well, with double layers seen in the beginning, effectively suppressing atmospheric dilution. Large-scale subsidence inversion aloft was also persistently present. In the first haze event, relatively low PM2.5 was observed in the afternoon over its first part, but a higher afternoon PM2.5 occurred over its second part, due possibly to the role of secondary aerosols. Comparing the two episodes, the PM2.5/PM<sup>10</sup> ratio was relatively low in the first episode because of more impact of biomass burning, which was in general agreement with back-trajectories and active fire hotspots. The second haze event, with little biomass burning in the region, was likely to be caused mainly by local anthropogenic emissions.

According to the above findings, certain policy-related implications are given as follows:


Future works on the urban haze problem in the study area may be extended to the following aspects: urban heat island, local recirculation, and systematic synoptic pattern classification [59], hourly mixing height variation and application of high-resolution light detection and ranging (LiDAR) if available [60], low-visibility events, and episode-specific source–receptor studies using source apportionment, chemical characterization (e.g., carbonaceous components in PM2.5), and air quality modeling.

#### *Int. J. Environ. Res. Public Health* **2020**, *17*, 9499

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1660-4601/17/24/9499/s1. Table S1. Selected air pollution studies in Greater Bangkok. Figure S1. Daily PM2.5 variation at the three air quality stations by year: (a) 2015; (b) 2016; (c) 2017. Figure S2. Example surface charts in EP1 corresponding to the four synoptic patterns classified: (a) Pattern 0; (b) Pattern 1; (c) Pattern 2; (d) Pattern 3. Figure S3. Example surface charts in EP2 corresponding to the four synoptic patterns classified: (a) Pattern 0; (b) Pattern 1; (c) Pattern 2; (d) Pattern 3. Figure S4. Diurnal variation of PM2.5 and Obukhov length (*L*) during: (a) EP1 in 2015; (b) EP2 in 2017. Figure S5. Vertical migration of each individual daily 48-h back-trajectory: (a) EP1 in 2015; (b) EP2 in 2017. Figure S6. Average daily active fire hotspots per pixel (0.5◦ by 0.5◦ ): (a) EP1 in 2015; (b) EP2 in 2017.

**Author Contributions:** K.M. initiated, conceptualized, and administrated the research. N.A. performed the entire analysis and drafted the manuscript. Both K.M. and N.A. jointly discussed the results. N.P.-E., E.K., T.B., S.P., and C.C. technically supported the research. B.D. contributed to the data processing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded mainly by the Energy Conservation and Promotion Fund (ENCONFUND) of the Ministry of Energy, Thailand, and partly by the Joint Graduate School of Energy and Environment (JGSEE) and the Postgraduate Education and Research Development Office (PERDO).

**Acknowledgments:** The authors sincerely thank the Pollution Control Department (PCD) and the Thai Meteorological Department (TMD) for the observational data used in the study. Seksan Sangkow (PCD) and Thiantawan Chulathipyachat (PCD) shared useful perspectives on air quality issues in Greater Bangkok and the use of air quality data. We also thank the two anonymous reviewers for the comments and suggestions that helped improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Meteorological Variables and Synoptic Patterns Associated with Air Pollutions in Eastern China during 2013–2018**

**Zhujun Dai 1,2, Duanyang Liu 1,3,4,\* , Kun Yu <sup>2</sup> , Lu Cao <sup>5</sup> and Youshan Jiang <sup>2</sup>**


Received: 24 February 2020; Accepted: 4 April 2020; Published: 7 April 2020

**Abstract:** Steady meteorological conditions are important external factors affecting air pollution. In order to analyze how adverse meteorological variables affect air pollution, surface synoptic situation patterns and meteorological conditions during heavy pollution episodes are discussed. The results showed that there were 78 RPHPDs (regional PM2.5 pollution days) in Jiangsu, with a decreasing trend year by year. Winter had the most stable meteorological conditions, thus most RPHPDs appeared in winter, followed by autumn and summer, with the least days in spring. RPHPDs were classified into three patterns, respectively, as equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT) according to the SLP (sea level pressure). RPHPDs under EQP were the most (51%), followed by ACF (37%); INT was the minimum (12%). Using statistical methods and meteorological condition data on RPHPDs from 2013 to 2017 to deduce the thresholds and 2018 as an independent dataset to validate the proposed thresholds, the threshold values of meteorological elements are summarized as follows. The probability of RPHPDs without rain was above 92% with the daily and hourly precipitation of all RPHPDs below 2.1 mm and 0.8 mm. Wind speed, RHs, inversion intensity(ITI), height difference in the temperature inversion(ITK), the lower height of temperature inversion (LHTI) and mixed-layer height (MLH) in terms of 25%–75% high probability range were respectively within 0.5–3.6 m s−<sup>1</sup> , 55%–92%, 0.7–4.0 ◦C 100 m <sup>−</sup><sup>1</sup> , 42–576 m, 3–570 m, 200–1200 m. Two conditions should be considered: whether the pattern was EQP, ACF or INT and whether the eight meteorological elements are within the thresholds. If both criteria are met, PM2.5 particles tend to accumulate and air pollution diffusion conditions are poor. Unfavorable meteorological conditions are the necessary, but not sufficient condition for RPHPDs.

**Keywords:** air pollution; synoptic situation pattern; meteorological variables; threshold values

#### **1. Introduction**

PM2.5 is the air particulate matter with an aerodynamically equivalent diameter of less than 2.5 µm [1]. Rising PM2.5 levels can lead to worsening air quality, seriously affecting human health. Epidemiological studies have shown associations between ambient air pollution and changes in heart rate variability (HRV) [2–5]. PM2.5 is characterized by the small size and lightweight. It can be transported to some faraway places from its source region, causing a wide range of heavy air pollution [6,7], leading to a wide range of human health hazards [8,9]. As a typical pollutant in heavy

pollution weather, PM2.5 has become a hot spot of research recently [10–13]. Beijing–Tianjin–Hebei region [14–17], the Yangtze River Delta region [18–22] and the Pearl River Delta region [23–25] are the three highly air-polluted regions in China. As a result, the incidence of air pollution-related diseases in major cities in these regions has been increasing year by year [8,9,26,27].

Although air pollution is usually linked with human activities, natural processes may also determine noticeable concentrations of hazardous substances in the low atmosphere. Many studies have suggested that continuous air polluting processes are jointly affected by pollutant emissions and the weather conditions that favor the accumulation of pollutants [22,23,28]. Comprehensive analyses of meteorological factors on air pollution were conducted by Zhou [29], who established the statistical model of PM2.5 concentrations and meteorological element field. The levels of pollutants may be reduced when emissions can be controlled. However, the impact of meteorological variables on concentrations measured may be marked, and these variables cannot be controlled. If pollutant emission is the internal factor of heavy air pollution, then the meteorological condition should be the external factor [30,31]. Deducing the threshold values of meteorological elements that are conducive to the pollution accumulation is very necessary.

The Yangtze River Delta where Jiangsu Province is located is one of China's most economically developed regions—and one of the most polluted regions. In this region, many questions have been studied by scientists, such as how meteorological conditions including wind speed affect the concentration of the atmospheric pollutants, how the air pollutants transported from other heavily polluted areas by the synoptic situation patterns and meteorological variables and how changes of atmospheric mixed-layer and temperature inversions affect the air quality. Since most of the past studies are cases studies, they have not provided systematic and comprehensive answers. Moreover, the thresholds of various meteorological elements are not given. Therefore, in order to in-depth knowledge and systematic study the different synoptic situation patterns and meteorological variables on air pollution, heavy air pollution from 2013 to 2018 in Jiangsu Province, China was analyzed in this article. According to the sea level pressure (SLP), the characteristics of the boundary meteorological element distribution are discussed in different seasons, different synoptic patterns and meteorological variables.

The remainder of this paper is organized as follows: Data and methods are described in Section 2. We analyzed synoptic patterns and meteorological conditions associated with heavy PM2.5 air pollutions in Section 3. The conclusions are given in Section 4.

#### **2. Materials and Methods**

#### *2.1. Study Area and Data Description*

The geographical location of Jiangsu Province, China (UTC/GMT+08:00), the distribution of meteorological stations and state-controlled environmental protection stations (SCEPSs) in 13 prefecture-level cities in Jiangsu were shown in Figure 1. In regards to the issue of mismatch between SCEPSs and meteorological stations, the data of SCEPSs nearest to 13 meteorological stations was adopted. According to diurnal variation characteristics and regional differences of heavy PM2.5 pollution, 13 cities in Jiangsu Province were divided into four regions: South Jiangsu (Suzhou, Wuxi and Changzhou), Coastal Jiangsu (Lianyungang, Yancheng and Nantong), Southwest Jiangsu (Zhenjiang, Yangzhou, Taizhou and Nanjing) and North Jiangsu (Xuzhou, Huaian and Suqian). Although SCEPSs were distributed in various environments (cities, suburbs, roadsides, parks, etc.), this paper focused on the regional PM2.5 heavy pollution, therefore the environmental differences were ignored.

The PM2.5 data of Jiangsu Province begins from 2013, therefore the data from 2013–2018 were used in this paper. Meteorological characteristics were discussed by the data from 2013 to 2017 and verified by the data of 2018. The dataset consists of four parts: (1) routine meteorological observation data including the hourly wind speed, wind direction, surface relative humidity (RHs) and precipitation at 13 meteorological stations distributed separately in prefecture-level cities in Jiangsu from 2013 to 2018; (2) the vertical distribution of temperature was obtained from the observational data from four radiosonde stations(Xuzhou, Sheyang, Nanjing and Baoshan) at 8 o'clock. Jiangsu just had three radiosondes, which were respectively Xuzhou, Sheyang, Nanjing, located in the North, Coastal, Southwest Jiangsu. Baoshan, the nearest radiosonde station to South Jiangsu, was located in Shanghai. Hence the data from the four radiosondes may be used to analyze the sounding situation of four regions in Jiangsu. (3) daily reanalysis data with spatial resolution of 2.5◦ × 2.5◦ from National Centers for Environmental Prediction (NCEP); (4) hourly PM2.5 mass concentration data provided by Jiangsu Environmental Monitoring Center from 2013 to 2018.

**Figure 1.** (**a**) The geographical location of Jiangsu province, China (UTC/GMT+08:00); (**b**) The distribution of 13 meteorological, 13 environmental monitoring and 4 radiosonde stations. (Lianyungang: in short LYG; Xuzhou: XZ; Nanjing: NJ; Nantong: NT; Suqian: SQ; Zhenjiang: ZJ; Taizhou: TZ; Huaian: HA; Yancheng: YC; Changzhou: CZ; Suzhou: SZ; Yangzhou: YZ; Baoshan: BS; Sheyang: SY).

Wind speed referred to the horizontal distance of air movement per unit time. Hourly wind speed took 1 second(s) as the time step to calculate the arithmetic average value of the 2minute(min) before the point hour. The sampling frequency of wind speed was 4 times/s.

Wind direction referred to the direction of wind. Hourly wind direction took 1 s as the time step to calculated the vector average value of 2 min before the point hour. The sampling frequency of wind direction was 1 times/s.

Relative humidity referred to the percentage of the actual partial pressure of water vapor in the air and the partial pressure of saturated water vapor under the same conditions. Hourly RHs were calculated as the arithmetic average value at 1 min before the hour. The sampling frequency of RHs was 30 times/min.

Hourly precipitation calculated the cumulative amount of 60 min before the hour. The sampling frequency of precipitation was 1 time/min.

There were four seasons in a year, namely spring (March, April, May), summer (June, July, August), autumn (September, October, November) and winter (January, February, December).

#### *2.2. Definition of Heavy Pollution Day*

The calculation method for daily average PM2.5 concentration was based on "Ambient Air Quality Standards" (GB3095-2012) (Ministry of Environmental Protection of the PRC, 2012b). If the daily average PM2.5 concentration above 150 µg m−<sup>3</sup> was observed in no less than 3 cities of Jiangsu Province in one day, it was defined as a PM2.5 HPDs in Jiangsu. If the daily average PM2.5 concentration above 150 µg m−<sup>3</sup> appeared in no less than three cities in a region, it was recorded as a regional PM2.5 HPDs (RPHPDs). The pollutant concentrations of the four regions in Jiangsu Province varied greatly, so the definition of RPHPDs was very necessary.

#### *2.3. Classification of Synoptic Patterns*

We analyzed the weather situation, especially SLP, to classified synoptic patterns on RPHPDs (as well as on non-RPHPDs) by the subjective approach, the classification method was the same as Peng, et al. [21] and Huth et al. [32–34]. The results showed that, surface conditions on RPHPDs were only dominated by three patterns, respectively were equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT) (Figure 2).

**Figure 2.** The synoptic situation for each type: equalized pressure (EQP: (**a**,**d**,**g**)), advancing edge of a cold front (ACF: (**b**,**e**,**h**)), inverted trough of low pressure (INT: (**c**,**f**,**i**)); 500 hPa (**a–c**), 850 hPa (**d–f**) and SLP (**g–i**). Among them, contours were pressure field (unit: hPa) and vectors were wind field (unit: m s−<sup>1</sup> ).

The three patterns were the necessary external condition for RPHPDs. There were any other synoptic conditions which frequently affect the flow over Jiangsu Province, but RPHPDs only occurred in the context of the three synoptic conditions. Nevertheless, non-RPHPDs were also observed under the three patterns, this was because that the synoptic pattern was the external factors causing RPHPDs [35].

(i) EQP (Figure 2g): when the cold air was blocked in the north, the domain was controlled by equalized pressure;


Upper and lower altitudes under three synoptic flow patterns were described as follows. At 500 hPa, the northwest westerly wind prevailed under EQP; the northwest wind prevailed under ACF; the straight west wind prevailed under INT. That is, in the meridional direction: ACF >EQP>INT (Figure 2a–c). At 850hPa, the west wind prevailed under EQP; the northwest wind prevailed under ACF; the southwest wind prevailed under INT (Figure 2d–f). In the sea-level-pressure field, Jiangsu was on the rear of the high pressure moving towards the sea under EQP; under ACF, Jiangsu was in the wedge at the front of cold high; under INT, Jiangsu was in the uniform pressure field at the front of low pressure. All the wind speeds were very low, and wind directions varied in different regions (Figure 2g–i).

Due to great differences of physical quantities in different seasons, RPHPDs were further classified (Table 1). In terms of four seasons and surface synoptic flow patterns, eight types were obtained: (a) spring INT, (b) summer EQP, (c) autumn EQP, (d) autumn ACF, (e) autumn INT, (f) winter EQP, (g) winter ACF, (h) winter INT.


**Table 1.** RPHPDs (regional PM2.5 pollution days) during 2013–2017 <sup>1</sup> . (equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT)).

<sup>1</sup> The date of RPHPDs,\: No RPHPDs in this type.

Considering the contingency of the data when there were less than three RPHPDs, characteristics of meteorological elements were studied only when RPHPDs were above three days in four regions and under different synoptic flow patterns (Table 1). The heavy pollution episodes which met the condition only happen in winter. The requirement was met in all four regions under winter EQP and winter ACF. While for winter INT, it was satisfied in only two regions—Southwest and North Jiangsu.

#### *2.4. Method for Calculating MLH*

The dry adiabatic curve method was used to calculate the mixing layer height (MLH). The dry adiabatic curve method was proposed by Holzworth [36,37] when studying the average maximum height of the mixing layer in some areas of the United States. As the required data were easy to obtain, the estimation method was widely used in China. The results showed that the dry adiabatic method was more practical than the methods established in the national standard GB/T 3840—91 [38] or Nozaki method [39,40] in calculating the MLH, and it was the most representative method. And its calculated results can denote the actual characteristics of the atmosphere to some extent. The daily maximum height of the mixing layer can be defined as the height below the intersection point of the dry adiabatic line and the temperature profile in the afternoon based on the daily radiosonde data at 0800 LST (Local Standard Time) and the maximum surface temperature data. This method was convenient for areas with radiosonde data and can be used for air pollution potential forecasting. Then, MLHs of Xuzhou, Sheyang, Nanjing and Baoshan radiosonde stations were calculated. Among them, Xuzhou, Sheyang and Nanjing were in North Jiangsu, Coastal Jiangsu and Southwest Jiangsu respectively; Baoshan station was a radiosonde station in Shanghai nearest to South Jiangsu. Thus, these four radiosonde stations were selected to represent the four regions (North Jiangsu, Coastal Jiangsu, Southwest Jiangsu and South Jiangsu) to analyzed the characteristics of MLH in each area.

#### *2.5. Method for Calculating Temperature Inversion*

Temperature inversion meant that the lower-layer temperature was lower than the higher-layer temperature. The lower height of a temperature inversion (LHTI) was the height (unit: m) nearest the surface when the temperature inversion layer develops. The upper height of a temperature inversion (UHTI) was the height (unit: m) at the top of the temperature inversion layer. The temperature difference in the temperature inversion was the upper temperature of the temperature inversion layer (UTTI) minus the low-level temperature of the temperature inversion layer (LTTI). The height difference in the temperature inversion (ITK) was the height of the LHTI minus the height of the UHTI. The intensity of the temperature inversion (ITI) was the temperature difference divided by the height difference multiplied by 100(unit: ◦C/100m): ITI = (UTTI − LTTI)/(LHTI − UHTI) \* 100.

#### **3. Results**

#### *3.1. Distribution Characteristics of RPHPDs*

Statistical characteristics were obtained by the data from 2013 to 2017 and verified by the data of 2018. PM2.5 HPDs in 13 prefecture-level cities of Jiangsu from 2013 to 2017 were shown in Figure 3. The largest number of PM2.5 HPDs appeared in Xuzhou (109 days), followed by Huaian (83 days), Taizhou (77 days), Suqian(70 days), Nanjing (64 days), Nantong (59 days), Lianyungang (57 days), Changzhou (57 days), Wuxi (53 days), Suzhou (52 days) and Zhenjiang (51 days); the number in Yancheng was the minimum (46 days).

During 2013–2017, there were a total of 78 RPHPDs in Jiangsu (Table 1), including 33 days in 2013, 17 days in 2014, 17 days in 2015, 7 days in 2016 and 4 days in 2017, that is to say, the number of RPHPDs decreased year by year. Source emissions as a major factor in pollution declined over the years, but were not the focus of this study. This paper analyzed the role of meteorological conditions as objective factors in pollution days, so the threshold value was one of the criteria of pollution intensity. We focused on the local pollution in Jiangsu Province, the transport from other provinces were not taken into account.

According to seasons, the largest number of RPHPDs was observed in winter (65 days), followed by autumn (7 days) and summer (5 days) and the least in spring (1 days).

Considering different synoptic flow patterns, from 2013 to 2017 (Table 1), RPHPDs under EQP was the most (40 days, 51% of the cases), followed by ACF (29 days, 37%); the number under INT was the minimum (9 days, 12%).

**Figure 3.** Heavy pollution days of 13 cities in Jiangsu Province, China (unit: days). (Lianyungang: in short LYG; Xuzhou: XZ; Nanjing: NJ; Nantong: NT; Suqian: SQ; Zhenjiang: ZJ; Taizhou: TZ; Huaian: HA; Yancheng: YC; Changzhou: CZ; Suzhou: SZ; Yangzhou: YZ; Baoshan: BS; Sheyang: SY).

#### *3.2. Frequencies of RPHPDs under Di*ff*erent Synoptic Flow Patterns*

The number of RPHPDs varied in different regions (Table 1). RPHPDs were 37, 26, 43 and 39 in South Jiangsu, Coastal Jiangsu, Southwest Jiangsu and North Jiangsu, respectively. The number of RPHPDs was the largest in Southwest Jiangsu, followed by North Jiangsu; the minimum number was in the Coastal Jiangsu.

The frequency of RPHPDs also varied in different seasons (Table 1). In spring, only one RPHPDs occurred in Southwest Jiangsu in May under INT. In summer, three RPHPDs and two RPHPDs respectively occurred in Southwest Jiangsu and North Jiangsu in June, which may be related to straw burning. In autumn, RPHPDs can occurred in all three patterns, but the probability was low and it only appeared in November. And the occurrence probability under EQP was higher than under ACF and INT. RPHPDs in South Jiangsu, Coastal Jiangsu and Southwest Jiangsu were one day, one day and three days respectively under EQP in autumn, while for ACF and INT in autumn, only one RPHPDs was observed in North Jiangsu and Southwest Jiangsu, respectively. RPHPDs appeared the most frequently in winter, mainly in December and January, followed by February. With regards to the three synoptic flow patterns, RPHPDs occurred more frequently under EQP and ACF and less frequently under INT. Yet no RPHPDs appeared under EQP in February, and no RPHPDs occurred in South Jiangsu under INT in winter.

The frequencies of RPHPDs under the three patterns of EQP, ACF and INT were also different (Table 1). Under EQP, RPHPDs were observed in summer, autumn and winter; under ACF, RPHPDs only occurred in autumn and winter, mainly in winter; under INT, RPHPDs occurred in spring, autumn and winter. The frequency of RPHPDs was high under patterns of g and g2 and low under INT. Moreover, the occurrence of RPHPDs was closely related to the circulation background. Under ACF, the high pressure was in the north and the cold air was active, thus RPHPDs only occurred in autumn and winter. Besides local pollution, large amounts of transportation and rapid accumulation of pollutants were also accompanied under ACF. Therefore, under EQP, the accumulation concentration of local pollutants was relatively high, and the continuous high-pressure control led to lower wind speed and stable stratification, and radiation inversions frequently happened in the night of sunny days. All these conditions finally resulted in the persistent accumulation of pollutants and a high frequency of RPHPDs. As for INT, the uniform pressure field were dominating with low wind speed and stable stratification. Located in the front of low pressure with the prevailing southeast wind, the surface convergence in Jiangsu was strong, which can easily led to the pollutant accumulation, formation

of inversion weather, rise of RHs and hygroscopic growth, increasing the occurrence probability of RPHPDs.

Furthermore, frequencies of RPHPDs under eight types varied in four regions (Table 1). For South Jiangsu, the number of RPHPDs for winter ACF were the highest (17 days), followed by winter EQP (15 days), and the minimum number was for autumn EQP (1 day). In Coastal Jiangsu, the number for winter EQP (12 days) were slightly larger than that for winter ACF (11 days), followed by winter INT (2 days) and autumn EQP (1 day). In Southwest Jiangsu, the number for winter ACF (16 days) was larger than that for winter EQP (14 days), followed by winter INT (5 days), summer EQP (3 days) and autumn EQP (3 days), and the minimum number was for spring INT(1 day) and autumn INT (1 day). In North Jiangsu, the number for winter EQP (17 days) was larger than that for winter ACF (14 days), followed by winter INT (5 days) and summer EQP (2 days), and the minimum number was for autumn ACF (1 day). In all four regions, the number for winter EQP or winter ACF was significantly larger than that for other types.

Considering the contingency of the data when there were less than three RPHPDs, winter were not the only season to be considered (Table 1). Winter EQP, winter ACF and winter INT were hereinafter referred to as EQP, ACF and INT for short. Thereby, RPHPDs can be categorized into ten types: EQP in North Jiangsu (EQP\_nth), EQP in South Jiangsu (EQP\_sth), EQP in Southwest Jiangsu (EQP\_sw), EQP in Coastal Jiangsu (EQP\_cst), ACF in North Jiangsu (ACF \_nth), ACF in South Jiangsu (ACF \_sth), ACF in Southwest Jiangsu (ACF \_sw), ACF in Coastal Jiangsu (ACF \_cst), INT in North Jiangsu (INT\_nth) and INT in Southwest Jiangsu (INT\_sw).

#### *3.3. Variation Characteristics of PM2.5 and Meteorological Elements for Ten Types*

#### 3.3.1. Concentration of PM2.5

Among the ten types, ACF\_nth had the highest values in PM2.5 concentration distribution, the high-value range of 25%–75% and the average value (Figure 4). In terms of 25%–75% high-value range and average value, for EQP, the concentration were not the highest in EQP\_cst and lowest in EQP\_sth; for ACF, the concentration were not the highest in ACF\_nth and lowest in ACF\_sth. For INT, the concentration in INT\_nth was higher than that in INT\_sw and the PM2.5 concentration in ACF\_sth was the lowest among the ten types. In the following section, characteristics of meteorological elements under the ten types were discussed.

**Figure 4.** PM2.5 concentration distributions of RPHPDs under the ten types. Rectangle: 25%–75% (quartile); rectangle midline: the median; end of the line: 5%–95%; ×: 1% and 99%; : average. (equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT).EQP in North Jiangsu (EQP\_nth), EQP in South Jiangsu (EQP\_sth), EQP in Southwest Jiangsu (EQP\_sw), EQP in Coastal Jiangsu (EQP\_cst), ACF in North Jiangsu (ACF \_nth), ACF in South Jiangsu (ACF \_sth), ACF in Southwest Jiangsu (ACF \_sw), ACF in Coastal Jiangsu (ACF \_cst), INT in North Jiangsu (INT\_nth) and INT in Southwest Jiangsu (INT\_sw)).

#### 3.3.2. Wind Direction and Speed

The prevailing wind direction was different under ten different types (Figure 5). The prevailing wind in EQP\_nth was the southeast direction (Figure 5a), the northwest direction in EQP\_sth (Figure 5b), the southeast direction in EQP\_sw and the southwest direction in EQP\_cst, both followed by north direction (Figure 5c,d). In ACF\_nth, the most frequent wind direction was north (Figure 5e). For ACF\_sth, the direction was northwest (Figure 5f). In ACF\_sw, the predominant wind direction was northwest, followed by northeast direction (Figure 5g). In ACF\_cst, the wind was northwest (Figure 5h). In INT\_nth, the most frequent wind direction was northwest and southeast (Figure 5i). Southeast wind prevails in INT\_sw (Figure 5j). Overall, wind directions under EQP in different regions were more dispersed, whereas, under ACF and INT, they were more concentrated. Wind speed under EQP was mostly below 3 m/s. Wind speed under ACF and INT were larger, mainly distributed in 0–5 m/s and 2–5 m/s, respectively.

**Figure 5.** Wind-rose map of flow vector wind in ten types ((**a**) EQP\_nth, (**b**) EQP\_sth, (**c**) EQP\_sw, (**d**) EQP\_cst, (**e**) ACF\_nth, (**f**) ACF\_sth, (**g**) ACF\_sw, (**h**) ACF\_cst, (**i**) INT\_nth, (**j**) INT\_sw). The color scale represents the average wind speed during the two minutes before the integral hour point. (equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT).EQP in North Jiangsu (EQP\_nth), EQP in South Jiangsu (EQP\_sth), EQP in Southwest Jiangsu (EQP\_sw), EQP in Coastal Jiangsu (EQP\_cst), ACF in North Jiangsu (ACF \_nth), ACF in South Jiangsu (ACF \_sth), ACF in Southwest Jiangsu (ACF \_sw), ACF in Coastal Jiangsu (ACF \_cst), INT in North Jiangsu (INT\_nth) and INT in Southwest Jiangsu (INT\_sw)).

#### 3.3.3. Diurnal Variation of RHs and Wind Speed

Higher RHs was more likely to cause heavy pollution. It may be caused by the accelerated reaction and hygroscopic growth of particulate matters under high RHs conditions [41], The rise of RH altered the particle size and thus visibility. The wind was an important dynamic factor affecting the diffusion of pollutants in the boundary layer. The wind direction determined the transport direction of pollutants in the atmosphere. Wind speed determined the diffusion and dilution speed of pollutants in the atmosphere, and variations of low-level wind direction and wind speed directly affected the convergence, divergence and concentration distribution of pollutants. Therefore, RHs and wind were important meteorological elements affecting the formation and development of heavy pollution.

The average diurnal variations of RHs and wind speed were calculated for ten types. Diurnal variations of RPHPDs in different regions were like each other, but there were still some differences (Figure 6). Except for an insignificant negative correlation observed between wind speed and RHs under INT\_sw (Figure 6e,f), negative correlation coefficients under the other 9 types were all in the range of 0.77–0.92, passing the 99.9% significance test.

**Figure 6.** Diurnal variations of wind speed and RHs for EQP(**a**,**b**), ACF(**c**,**d**) and INT(**e**,**f**); red: North Jiangsu, black: South Jiangsu, green: Southwest Jiangsu, blue: Coastal Jiangsu, X-axis was time (unit: h); Rectangle: 25%-75% (quartile); rectangle midline: the median; end of the line: 5%-95%; ×: 1% and 99%; : average. (equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT)).

Diurnal variations of RHs and wind speed under EQP and ACF were more significant than INT (Figure 6a–d). The RHs decreased greatly and the wind speed increased notably in the afternoon, reaching a minimum values during 1700–1900 LST. Diurnal variation ranges of RHs and wind speed also varied in different regions. It can be seen that ranges of RHs diurnal variation were larger in Coastal Jiangsu and South Jiangsu than those in North Jiangsu and Southwest Jiangsu. Wind speeds in Coastal Jiangsu and South Jiangsu were slightly larger than that in North Jiangsu for ACF, and the range of RHs diurnal variation in Coastal Jiangsu was the smallest. The wind speed for ACF began to fluctuated and increased at 1600 LST, reaching its peak at 2000 LST and remained at a high level at night.

Compared with EQP and ACF, the diurnal variation of INT was insignificant, and diurnal variation ranges of RHs and wind speed were smaller, but wind speed was larger (Figure 6e,f). RHs in Southwest Jiangsu were significantly larger than North Jiangsu for INT.RHs decreased slightly at 1600 LST and reached the bottom value between 1900 LST and 2000 LST. Wind speed in INT\_sw was within 1.8–2.7 m s <sup>−</sup><sup>1</sup> and RHs was within 76%–84% and those for INT\_nth were 1.6–2.9 m s−<sup>1</sup> and 65%–73%, respectively.

#### *3.4. Threshold Values of Meteorological Elements Causing RPHPDs*

Using statistical methods and meteorological condition data with RPHPDs from 2013 to 2017 to deduce the thresholds, meteorological data of 2018 were used as an independent dataset to validate the proposed thresholds. The threshold values of meteorological elements for RPHPDs were analyzed as follows.

#### 3.4.1. Precipitation

Precipitation rarely occurred in winter on RPHPDs in Jiangsu. There were only three days, two days, one day and three days of precipitation, respectively in South Jiangsu, Coastal Jiangsu, Southwest Jiangsu and North Jiangsu on RPHPDs during 2013–2017. Probabilities of no precipitation on RPHPDs were as high as 92%, 92%, 98% and 92%, respectively. That is to say, the probability of RPHPDs without rain was above 92% with the daily and hourly precipitation below 2.1 mm and 0.8 mm. This was because the wet deposition brought by heavy precipitation can reduce the PM2.5 concentration. When the light precipitation occurred on RPHPDs in Jiangsu, RH values increased to 83%–100%, high RH promoted the hygroscopic growth of pollutants [42], which can be explained by the Köhler's theory [43]. The rise of RH and hygroscopic growth altered the particle size and aggravated the pollutant accumulation [44–46].

#### 3.4.2. Wind Speed and RHs

In this section, the threshold values of wind speed and RHs on RPHPDs were investigated. First, time points of RPHPDs were selected. Then average, maximum and minimum values of wind speed and RHs were calculated. Finally, the threshold values of wind speed and RHs were counted in terms of 25–75% high probability range when RPHPDs occurred. The results were shown in Figure 7. It can be seen from Figure 7a that wind speeds in ten types were very low, with an average of 1.3–2.6 m s−<sup>1</sup> . Wind speeds can reached more than 5 m s−<sup>1</sup> at some time points in all ten types. In terms of 25%–75% high-value range of RPHPDs, the wind speed varied within 0.5–3.6 m s−<sup>1</sup> . The result was in good agreement with the conclusion of Dai et al. [47]. Specifically, the maximum wind speed was observed in INT\_nth, followed by ACF\_cst and INT\_sw. Wind speeds in INT\_nth were mainly within 1.5–3.5 m s <sup>−</sup><sup>1</sup> and the maximum value reaches 7.5 m s−<sup>1</sup> . Average wind speeds in INT\_nth, INT\_sw and ACF\_cst were all above 2 m s−<sup>1</sup> .

The wind speed under EQP was minimum with the average wind speed within 1.3–1.5 m s−<sup>1</sup> and the wind speed within 0.5–2.3 m s−<sup>1</sup> in terms of 25%–75% high probability range. Wind speeds in different regions were basically the same. The average wind speed was 1.7–2.4 m s−<sup>1</sup> under ACF, which slightly higher in South Jiangsu. The wind speed in terms of 25–75% high probability range was within 0.8–2.7 m s−<sup>1</sup> . In a word, the rank of patterns based on the wind speed value on the RPHPDs was INT > ACF > EQP.

The distribution of RHs in RPHPDs was shown in Figure 7b. RHs for ten types varies greatly (21%–100%) with the average RHs within 68%–82%, the maximum was in INT\_sw, followed by INT\_nth; the minimum was in ACF\_sth. Comparing EQP and ACF, INT had the maximum average RHs. It was found that the average RHs under EQP was higher than ACF except in Southwest Jiangsu. RHs on RPHPDs in terms of 25%–75% high probability range was within 55%–92%, with the maximum under INT (65%–92%) and RHs for EQP was basically equivalent to ACF.

**Figure 7.** Wind speed (**a**) and RHs (**b**) on RPHPDs (regional PM2.5 pollution days) under the ten types. Rectangle: 25%–75% (quartile); rectangle midline: the median; end of the line: 5%–95%; ×: 1% and 99%; : average.

To sum up, for INT\_nth and INT\_sw, the light wind was more likely to cause heavy pollution than the still wind. This was because, under INT, Jiangsu was in the uniform pressure field in the front of low pressure (Figure 2i) and the southeast wind prevails (Figure 5i,j). The light southeast wind was more likely to form a weak wind field convergence near the surface than the still wind. Humidities in INT\_nth and INT\_sw were also higher than those in other types due to the high RHs of southeast wind. Under EQP, the entire Jiangsu Province was in the center of high pressure (Figure 2g), and wind speeds in four regions were very low. RPHPDs under ACF often occurred before the arrival of cold air moving along the middle or west path (Figure 2h). Dry cold air led to a very low RHs under ACF whose wind speed was slightly higher than that under EQP. Coastal Jiangsu were in the front of the cold front, coupled with the land-sea thermal difference, the wind speed in Coastal Jiangsu was slightly larger than that in other areas.

For the ten patterns, light wind was more likely to cause RPHPDs than still wind. The mechanism was still wind urge pollutants to accumulate in the same place, while light wind prompted pollutants to accumulate in a small area, resulting in weaker vertical shear and vertical mixing of pollutants, which will increased the scope of air pollution. Ultimately, it affected the health of more and more people. The strong wind increased the horizontal vertical shear, and the vertical mixing of pollutants was stronger, which was conducive to the diffusion of pollutants to the upper air. The explanation was consistent with the conclusion of Zhang, et al. [48].

#### 3.4.3. ITI, ITK, LHTI and MLH

Under thermal inversion, especially low-level inversion, the vertical movement and turbulent exchange in the atmosphere were restrained, preventing the vertical transportation and horizontal diffusion of water vapor, smoke and other pollutants in the atmosphere. Thus, it led to persistent heavy pollution. The influence of inversion on the pollution diffusion was related to its ITI, ITK and LHTI [49–52].

Four radiosonde stations, Xuzhou, Sheyang, Nanjing and Baoshan, were selected as representative stations for North Jiangsu, Coastal Jiangsu, Southwest Jiangsu and South Jiangsu, respectively. The radiosonde data at 0800 LST was used to obtain the characteristics of ITI, ITK and LHTI at the lowest level in ten types of RPHPDs (Figure 8). Statistical results showed that except for type EQP\_sth

on January 18, 2014, inversion were observed in all RPHPDs. Temperature inversion and heavy pollution were bidirectional feedback mechanisms. A capping inversion was not conducive to the diffusion and convection of pollutants and was easy to form a static and stable weather condition that lead to accumulation of pollution. The occurrence of heavy pollution was often accompanied by temperature inversion. When pollution accumulates to a certain level, aerosol pollution can lead to temperature inversion (Zhong, 2019) [53]. ITI was 0–6.5 ◦C 100 m <sup>−</sup><sup>1</sup> , with an average of 1.6–3.1 ◦C 100 m <sup>−</sup><sup>1</sup> . ITK was within 0–651 m, with an average value 142–354 m. LHTI was within 0–1015 m, with an average value 71–407 m. When RPHPDs occurred, ITI, ITK and LHTI were respectively within 0.7–4.0 ◦C 100 m <sup>−</sup><sup>1</sup> , 42–576 m and 3–570 m in terms of 25%–75% high probability range.

**Figure 8.** Inversion intensity(ITI) (**a**), height difference in the temperature inversion(ITK) (**b**), the lower height of temperature inversion (LHTI) (**c**) and mixed-layer height (MLH) (**d**) on RPHPDs (regional PM2.5 pollution days) under the ten types. Rectangle: 25%–75% (quartile); rectangle midline: the median; end of the line: 5%–95%; ×: 1% and 99%; : average.

Specifically, the average ITI in EQP\_cst was the strongest and that in ACF\_sth was the weakest (Figure 8a). Under EQP and ACF, ITI was the strongest in Coastal Jiangsu, and the weakest in Southwest Jiangsu. However, under INT, ITI in Southwest Jiangsu was slightly higher than that in North Jiangsu. For the same region, the average ITI under ACF was lower than that under EQP and INT, and the average ITK under ACF was lower than that under the other two patterns (Figure 8b). Under EQP and ACF, the largest ITK was in North Jiangsu, followed by Coastal Jiangsu, Southwest Jiangsu and South Jiangsu in turn. Under INT, ITK was larger in Southwest Jiangsu than that in North Jiangsu. Only in ACF\_sth, ACF\_cst and INT\_nth, LHTI on a few RPHPDs was relatively higher. LHTI in other types were all below 300 m in terms of 25%–75% high probability range (Figure 8c).

The mixing layer referred to the atmospheric layer with sufficient turbulence below the discontinuous-turbulence interface [54–56]. The MLH affected the vertical diffusion of pollutants. Figure 8d showed that the MLH was relatively low (83–2284 m) when RPHPDs occurred. And the average value was less than 1 km (364–993 m). The MLH was within 200–1200 m in terms of 25%–75%

high probability range when regional RPHPDs occurred. The MLH was low and the atmospheric dispersion and dilution capability were weak, which was conducive to the accumulation of pollutants. Specifically, the mixing layer was the lowest under INT (INT\_nth and INT\_sw), and heights were basically the same under EQP and ACF, and MLH were in ranges of 246–740 m, 700–1200 m, 500–1100 m and 202–1100 m in terms of 25%–75% high probability range when RPHPDs occurred in North Jiangsu, South Jiangsu, Coastal Jiangsu and Southwest Jiangsu, respectively.

The ITI, ITK, LHTI and MLH were within 0.7–4.0 ◦C 100 m−<sup>1</sup> , 42–576 m, 3–570 m, 200–1200 m, respectively. The threshold values of inversion and MLH was coincide with the results of Zhong J.T. [56], Miao Y.C. et al. [57], Huang X. et al. [58] and Ding A. J. et al. [59]. According to Ding A. J. et al. [59], the threshold values ranged extremely, the lower MLH or the greater ITI, the poorer air pollution diffusion condition. Unfavorable meteorological conditions were the necessary but not sufficient condition for RPHPDs.

#### *3.5. Reliability Test of Threshold Values*

The threshold values were deduced from the winter of 2013–2017, so the heavy pollution data in of 2018 as an independent dataset was used to conduct the reliability tests. In the winter of 2018 consists of January, February and December, with a total of 90 days, there were 29 days which meteorological conditions met both criteria: (1) the pattern was EQP, ACF or INT, (2) the eight meteorological elements were within the thresholds. There were only 11 RPHPDs that occurred in January among the 29 days, accounted for 38%. The proportion was not high, because the meteorological condition was one of factors leading to RPHPDs.

The threshold values of meteorological elements for RPHPDs in Jiangsu in winters of 2013–2017 were tested in Section 3.3. The results (Table 2) show that in the winter of 2018, RPHPDs only occurred in January, and the corresponding synoptic patterns were g and j. The frequency of RPHPDs under EQP was the highest, while only one RPHPDs under INT occurred. The number of RPHPDs in North Jiangsu was the largest, followed by Southwest Jiangsu and South Jiangsu, and no RPHPDs occurred in Coastal Jiangsu. The daily and hourly precipitation of RPHPDs in 2018 was less than 2.1 mm and 0.8 mm. The wind speed, wind direction, RHs, ITI, ITK and LHTI were all distributed within the aforementioned threshold ranges. That is to say, ITI is within 0.7–4.0 ◦C 100 m <sup>−</sup><sup>1</sup> ; ITK was within 42–576 m; LHTI was within 3–570 m; MLHs in North Jiangsu, South Jiangsu and Southwest Jiangsu were also within the range.

**Table 2.** Dates of RPHPDs (regional PM2.5 pollution days) and variation range of eight meteorological elements in Jiangsu in the winter of 2018. (equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT).EQP in North Jiangsu (EQP\_nth), EQP in South Jiangsu (EQP\_sth), EQP in Southwest Jiangsu (EQP\_sw), EQP in Coastal Jiangsu (EQP\_cst), ACF in North Jiangsu (ACF \_nth), ACF in South Jiangsu (ACF \_sth), ACF in Southwest Jiangsu (ACF \_sw), ACF in Coastal Jiangsu (ACF \_cst), INT in North Jiangsu (INT\_nth) and INT in Southwest Jiangsu (INT\_sw)).


After testing, it was found that when the synoptic flow pattern were EQP, ACF or INT and eight meteorological elements were within the threshold ranges, PM2.5 particles tend to accumulate, and air pollution diffusion conditions were poor. Unfavorable meteorological conditions were the necessary but not sufficient condition for RPHPDs. For example, in January of 2018, the probabilities of RPHPDs in South Jiangsu, North Jiangsu and Southwest Jiangsu were 2/3, 2/3 and 1/2, respectively, but RPHPDs was not necessarily to occurred as the meteorological condition was only one of factors leading to heavy pollution.

#### **4. Conclusions**

Continuous air polluting processes were jointly affected by pollutant emissions and the weather conditions that favor the accumulation of pollutants. Steady meteorological conditions were important external factors affecting air pollution. Deducing the threshold values of meteorological elements that were conducive to the pollution accumulation was very necessary to achieve a better control of air pollution. This study provided an in-depth, systematic and comprehensive research about the different synoptic situation patterns and meteorological variables on air pollution, the heavy air pollution from 2013 to 2018 in Jiangsu Province, China.

Surface conditions on RPHPDs were only dominated by three patterns, respectively were equalized pressure (EQP), advancing edge of a cold front (ACF) and inverted trough of low pressure (INT). Most RPHPDs appeared under EQP, followed by ACF, and the last days were under INT. The three patterns were the necessary condition for RPHPDs. There were any other synoptic conditions which frequently affect the flow over Jiangsu Province, but RPHPDs only occurred in the context of the three synoptic conditions. Nevertheless, non-RPHPDs were also observed under the three patterns, since the synoptic pattern was the external factors causing RPHPDs.

The threshold values of meteorological elements when RPHPDs occurred in terms of 25%–75% high probability range were as follows: The daily precipitation was below 2.1 mm; the wind speed was within 0.5–3.6 m s−1; the RHs was within 55%–92%; the thermal ITI was within 0.7–4.0 ◦C 100 m−<sup>1</sup> ; ITK was within 42–576 m; LHTI was within 3–570 m; MLH was within 200–1200 m. The MLH in terms of 25%–75% high probability range when heavy pollution occurred in South Jiangsu, Coastal Jiangsu and Southwest Jiangsu was within ranges of 246–1000 m, 163–1300 m, 162–1200 m and 200–1200 m, respectively. After validation with one-year data, threshold values were proved to be relatively reliable.

We focused on the discussion of meteorological conditions for air pollution diffusion in this paper. the source emissions as a major factor in pollution should be the main solution for air pollution control.

The influence of regional transport was not discussed. In future, if the regional PM2.5 concentrations continue to decrease, the threshold values would remain applicable. When the threshold values were met, it represents PM2.5 particles tend to accumulate, and air pollution diffusion conditions were poor. Unfavorable meteorological conditions always were the necessary but not sufficient condition for RPHPDs.

We took five-year data to determine threshold values, and only one-year data to finished the reliability test. RPHPDs cases of Jiangsu in 2018 were few and more cases were needed to verify the reliability of threshold values, the combined indicators need to be further constructed. In addition, this threshold needs localization correction when it extended to other regions in China or other regions around the globe. The aim was to provide an easy and quick approach to carry out heavy pollution forecasting and early warning services.

**Author Contributions:** Conceptualization, D.L.; validation, K.Y.; formal analysis, Z.D. and L.C.; project administration, D.L.; resources, Y.J.; writing—original draft preparation, Z.D.; writing—review and editing, D.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was jointly supported by the National Key Project of MOST (2016YFC0203303 and 2016YFC0201903), the National Natural Science Foundation of China (91544231 and 41575135), the 333 Project of Jiangsu Province (BRA2016565), the Jiangsu Meteorological Bureau Fund (KM202010), the Jiangsu Meteorological Bureau Youth Fund (KQ201913, KZ201902).

**Acknowledgments:** The authors would like to thank Jiansu Wei from the Meteorological Observatory of Jiangsu Province, China, Cong Li and Rong Lu from the Nanjing Meteorological Bureau of Jiangsu Province, China and Hongbin Wang from the Jiangsu Institute of Meteorological Sciences, China for their advice on text revision. Thanks to Nanjing Environmental Monitoring Center of Jiangsu Province for providing the PM2.5 data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Combining Cluster Analysis of Air Pollution and Meteorological Data with Receptor Model Results for Ambient PM2.5 and PM<sup>10</sup>**

#### **Héctor Jorquera 1,2,\* and Ana María Villalobos <sup>3</sup>**


Received: 30 September 2020; Accepted: 12 November 2020; Published: 15 November 2020 -

**Abstract:** Air pollution regulation requires knowing major sources on any given zone, setting specific controls, and assessing how health risks evolve in response to those controls. Receptor models (RM) can identify major sources: transport, industry, residential, etc. However, RM results are typically available for short term periods, and there is a paucity of RM results for developing countries. We propose to combine a cluster analysis (CA) of air pollution and meteorological measurements with a short-term RM analysis to estimate a long-term, hourly source apportionment of ambient PM2.5 and PM10. We have developed a proof of the concept for this proposed methodology in three case studies: a large metropolitan zone, a city with dominant residential wood burning (RWB) emissions, and a city in the middle of a desert region. We have found it feasible to identify the major sources in the CA results and obtain hourly time series of their contributions, effectively extending short-term RM results to the whole ambient monitoring period. This methodology adds value to existing ambient data. The hourly time series results would allow researchers to apportion health benefits associated with specific air pollution regulations, estimate source-specific trends, improve emission inventories, and conduct environmental justice studies, among several potential applications.

**Keywords:** air pollution; source apportionment; cluster analysis; health risks; residential wood burning; sustainable urban development

#### **1. Introduction**

Ambient air pollution is a major environmental risk worldwide. Current estimates of premature mortality brought by population exposure to ambient PM2.5 (particulate matter with aerodynamic diameter below 2.5 µm) vary between five million [1] to ten million [2] annual deaths. This burden of disease has prompted a worldwide effort to regulate ambient PM2.5 concentrations. A customary set of public policies targeting that overarching goal includes: ambient air pollution monitoring, emission inventories, air pollution modeling, health effects studies, and economic valuation of regulation measures (emission standards, ambient air quality standards, market-based instruments, etc.).

Ambient PM2.5 is directly emitted by traffic, industrial, commercial, and residential sources in urban zones, and by agriculture, mining, industrial, and natural sources elsewhere. However, the resulting ambient PM2.5 concentrations are far more complex to characterize. Some gaseous pollutants (sulfur and nitrogen oxides) undergo atmospheric oxidation, leading to sulfuric and nitric acids, which, in turn, react with ammonia to produce ammonium sulfate and nitrate particles, forming secondary PM2.5 [3]. Furthermore, volatile (VOC) and semi-volatile (SVOC) organic compounds are

also oxidized in the atmosphere to become higher molecular weight compounds that condense to the particulate phase, generating secondary organic PM2.5 [4]. The large number of organic compounds released, along with their wide range of gas/particle phase partitioning, add further complexity to ambient PM2.5 [5]. Meteorology also plays a role in the previously mentioned chemical and physical processes, adding seasonality to ambient PM2.5 at different time scales: wind advection moves pollutant away from their sources, temperature enhances atmospheric chemical reactions and VOC emissions, solar radiation modulates photochemical processes, relative humidity modifies fugitive PM emissions, atmospheric stability controls the degree of vertical mixing of air pollutants, rainfall brings particles and gases down to the ground, etc.

Hence, a major hurdle in quantifying population exposure to ambient PM2.5 lies on the wide range of sources and processes that drive ambient PM2.5 concentrations, and their variability on different time scales: diurnal, weekly, seasonal, long-term. Although ambient PM2.5 monitoring is becoming widespread worldwide [6], a quantitative knowledge of the major sources accounting for those observed concentrations is still scarce. To obtain such knowledge, there are two customary modeling tools to conduct such a task: dispersion models and receptor models, which are described next.

Dispersion models (DM) use atmospheric emission inventories to model ambient PM2.5 at different spatial scales: urban, regional, continental, and global [7,8]. They use atmospheric modeling to simulate how those emission sources disperse, react, deposit on the surface as they are advected by the wind field, and mixed by atmospheric turbulence. In principle, these models provide a detailed spatial and temporal representation of ambient PM2.5 concentrations to evaluate regulatory policies at a national level. However, there are several sources of uncertainty in the emission inventories themselves and also in the parametrization of physical processes occurring in the lower atmosphere. Those two uncertainties propagate in the simulated ambient concentrations. For instance, urban zones in complex terrain pose a great challenge to current meteorological modeling tools [9], and the low-wind, stable meteorological conditions that drive air pollution episodes are difficult to be represented in current models [10]. In a review of photochemical and PM2.5 dispersion modeling applications in the US and Canada [11], the percentage of observed variance explained by the dispersion model (r<sup>2</sup> ) varies between 0.2 and 0.6 for hourly ozone and between 0.1 and 0.6 for PM2.5 and its major species. These are interquartile ranges and, in some cases, the model performance is better. These results represent the typical uncertainty expected for this DM approach.

Receptor models (RM) use chemical speciation of ambient PM2.5 samples, taken on filters (daily or longer time-accumulated samples), to estimate the major sources contributing to ambient PM2.5 mass concentrations [12–16]. They are based on the mass conservation principle applied to a set of tracer species (elements, isotopes, ions, carbonaceous aerosol, organic molecules). The required chemical speciation analyses are expensive, so most RM studies have been conducted in developed countries and have short-term results with a small proportion of long-term (multiyear) studies. A challenge for RM studies is to cover most urban zones in any given country to obtain a national level assessment of pollution sources to evaluate a regulatory policy. Regarding the uncertainty associated in RM results, one way of estimating it is to conduct a round robin exercise in which the very same data set is analyzed by different research groups. In one of these exercises conducted with a synthetic database [12], it was found that 87% of 344 RM source contribution estimates (SCE), generated by 47 participant groups, were within 50% of the correct result. Larger sources (SCE > 10% of total PM mass) were better quantified by RM, while smaller sources (SCE near 1% of total PM mass) were a challenge for most RM. In another European study (with an actual chemical speciation database), 38 RM were applied by 33 participant groups [15]. Additionally, 72% of the estimated time series of SCE were within 1-sigma uncertainty and 91% of the estimated SCE passed the z-score test. Therefore, receptor model results for SCE have a typical upper bound uncertainty of 50%.

A more recent generation of continuous ambient PM2.5 instruments has been developed. They provide continuous information on chemical speciation [17–19] and optical properties [20–24], which also leads to source apportionment, even at a finer temporal scale (minutes, hours) than

filter-based RM applications, which typically correspond to daily or weekly samples. That kind of instrumentation is more expensive than traditional ambient monitors for regulated pollutants, so it is currently less widespread than RM applications.

The preceding paragraphs show that a quantitative apportionment of major sources of ambient PM2.5 in urban areas worldwide is yet to be accomplished, particularly in developing countries. This is the current knowledge gap.

With an ever-increasing global ambient air pollution monitoring by regulatory authorities, the possibility of extracting information out of these databases has been on the rise, estimating temporal trends in air pollution [25–29]. Furthermore, more information is expected to be collected due to the rise of low-cost ambient monitoring [30–37]. More insight into the likely sources of ambient PM2.5 and PM<sup>10</sup> can be achieved by combining meteorological and air pollution measurements. Bivariate plots of ambient concentration as a function of wind speed and direction provide information on sources contributing to ambient concentrations [38,39]. For instance, tall stack and area sources contribute most when wind speeds are higher and lower, respectively [3]. One step further consists in performing a cluster analysis (CA) of this type of bivariate representation, so that major clusters contributing to ambient concentrations may be visualized and explored [40]. One limitation of this CA approach is that there is no simple recipe to know how many clusters should be chosen in a given analysis [41]. In other words, more information is needed to decide how many clusters may be resolved in any given set of ambient air pollution data.

In a review of CA applied to air pollution analysis between 1980 and 2019 [41], it is shown that most CA applications to surface observations (i.e., air quality monitoring) consisted in exploring spatial associations among monitor sites within the same city, region, or country. The temporal evolution of those clusters was focused on their seasonal variability. In one publication [42], chemical speciation of ambient PM2.5 was included in the cluster analysis to classify daily samples according to differences in PM2.5 chemical composition. Then back trajectory analyses were developed to estimate air mass origins associated with each resolved cluster.

Our goal is to show that, by combining short-term RM analyses with results gotten from long-term (multiyear) CA, for the same location, it is possible to identify the major sources (clusters) in ambient air pollution data therein. With this approach, the RM results provide additional information to decide how many clusters should be considered in the CA, so both methods complement each other. We show a proof of the concept for this combined approach in three case studies of urban zones with widely different dominant sources and meteorological conditions. As far as we are aware of, this is the first time this joint analysis has been proposed for ambient air pollution data.

This combined analysis is facilitated with a set of simple rules to identify sources: PM2.5/PM<sup>10</sup> ratios, dependence of ambient concentrations with temperature, wind speed and relative humidity, and the location of the monitoring site with respect to obvious sources: highways, industrial zones, etc. All these rules come from the literature of RM studies worldwide—see [16,43,44] for references —along with well-known results for the dispersion of stack emissions under different meteorological conditions [3] and the main features of motor vehicle emissions [45–52] as well as other combustion emissions [53–58].

The outcome of this proposed computational process is a long-term (multiyear), hourly time series of source contributions to ambient PM2.5 and PM<sup>10</sup> (and other measured air pollutants as well) suitable for different research purposes: exploring association of health effects with specific sources (example: traffic), tracking source trends, evaluate effectiveness of sector regulations, constrain sector emissions through dispersion modeling, conduct environmental justice analysis, further analysis of long-range sources using backward wind trajectories, etc. A major implication of our work is that ambient air quality databases are a rich source for further air pollution analysis worldwide, especially for source apportionment of PM2.5 and PM10.

#### **2. Materials and Methods**

#### *2.1. Computational Methodology*

We use the open access R software environment, including the openair package for air pollution analysis [59]. Openair has several dedicated functions for air pollution analysis, including the polarCluster function that groups ambient pollutant concentrations using bivariate plots. Briefly, polarCluster performs a k-means cluster analysis on the vectors [u, v, c] where u, v are the zonal and meridional wind components (sometimes referred to as the east and north wind components) and c as the pollutant concentration under analysis (PM2.5, PM10, etc.). Since all these variables have different scales, they are all standardized before the k-means algorithm proceeds [40]. This k-means clustering algorithm is well-known and has been extensively used in air pollution analyses [41], as mentioned above.

Regarding the receptor model methodology, we use already published RM results for two Chilean cities: Santiago and Temuco [60,61]. For Santiago, USEPA Positive Matrix Factorization (PMF, version 5.0) was the RM used [62], while for Temuco, the USEPA Chemical Mass Balance (Version 8.2) was the RM used [63]. Both models solve the following mass balance equation for p sources [13].

$$X\_{ij} = \sum\_{k=1}^{p} g\_{ik} f\_{kj} + e\_{ij} \tag{1}$$

where *Xij*, *gik*, and *fkj* are matrices whose entries are the j-th species mass concentration measured in the i-th sample, the mass concentration from the k-th source contributing to the i-th sample, and the j-th species mass fraction in the k-th source, respectively. *p* is the total number of resolved sources. The residuals *eij* are assumed random and normally distributed.

The Chemical Mass Balance (CMB8.2) model solves Equation (1) for the case when the source profiles {*fkj*} have been experimentally measured for the k-th source. Hence, Equation (1) is solved for {*gik*} using an effective variance least squares approach [64]. The Chemical Mass Balance model assumes that all major sources in the study area have been identified and included in its input. This model is usually run for different combinations of sources until a satisfactory solution is achieved in terms of the percentage of observed variance explained by the model.

The Positive Matrix Factorization (PMF5) model solves Equation (1) without making any assumption regarding source profiles {*fkj*}. However, in such a case, there are more unknowns than equations in Equation (1). This means additional information must be supplied into the model [65]. Usually, source compositions {*fkj*} and contributions {*gik*} are required to be non-negative. This is the case of software PMF5, which minimizes the following function.

$$Q = \sum\_{i=1}^{n} \sum\_{j=1}^{m} \left[ \left( X\_{ij} - \sum\_{k=1}^{p} g\_{ik} f\_{kj} \right) / \sigma\_{ij} \right]^2 \tag{2}$$

where σ*ij* is the estimated uncertainty in the j-th species at i-th PM sample. In this approach, all observations are individually weighted by their respective uncertainties. The above minimization is carried out including the previously mentioned non-negative constrains upon compositions {*fkj*} and contributions {*gik*}.

More details on the information used for each RM application (sampling period, chemical analyses, etc.) are provided in the respective references above [60,61].

#### *2.2. Case Studies*

We have selected, as case studies to test the proposed methodology, three urban zones in Chile with widely different meteorological conditions and dominant sources in which RM results are known. Calama (22◦27'S, 68◦55'W, population 2017: 158,600 [66], elevation: 2500 m) is a city located within the Atacama Desert (Koppen-Geiger classification BWk [67]) and the city is close to a large mining district. The major pollution problem is ambient PM10, given the desert conditions and rather strong winds that promote road dust resuspension by traffic and aeolian dust blown from the desert environment. Santiago (33◦28'S, 70◦40'W, population 2017: 6.85 million [66], elevation: 600 m) is the capital city, and extends along a valley surrounded everywhere by ranges and the Andes cordillera to the east. The climate is a warm temperate with dry summer (Koppen-Geiger classification CSb). Air quality regulations in Santiago were the first to be enacted in Chile and they have been successful in curbing down ambient PM2.5 [68–71]. Nonetheless, ambient PM2.5 in Santiago currently exceeds World Health Organization (WHO) guidelines and Chilean ambient air quality standards (AAQS), and major PM2.5 sources are traffic, industry, commercial, and residential sources. The third case study is the city of Temuco (38◦44'S, 72◦35'W, population 2017: 308,600 [66], elevation: 350 m) with a warm temperate fully humid climate (Koppen-Geiger classification Cfb). In this city, residential wood burning emissions rise in colder months, leading to severe ambient PM2.5 concentrations [60] that increase indoor concentrations as well [72].

Ambient air pollution (PM10, PM2.5, CO, SO2, NOx) and surface meteorological measurements were downloaded from Chile's Air Quality Information System (https://sinca.mma.gob.cl/). Table 1 shows the summary of information used to conduct the analyses, where one monitoring station per city was chosen. Figure 1 shows the locations of the three monitoring sites used in our analyses. Data were screened for obvious outliers and these were removed before the computational analyses.



<sup>1</sup> For the Santiago, Las Condes (LAC) site, meteorology is only available from 01/2004 through 06/2012.

**Figure 1.** Geographical locations of the three cities analyzed.

#### *2.3. Simple Comparison Rules*

To map the outcome of the CA results with the RM results for the same site, we propose to use a set of rules. Some of these rules have been traditionally used in identifying RM results, such as looking at the weekly seasonality of source contributions to identify traffic or industrial sources [61,70,71,73]. The former is lower over weekends while the latter is not. Other rules use the hourly resolution of CA results to generate scatter plots and time-variability plots of ambient PM2.5, PM10, and gases (CO, NOx, SO2), stratified by cluster, to check for specific patterns. For instance, large industrial stacks contribute more to ambient concentrations when wind speeds are higher [40] and ship emissions contribute more when ambient temperature rises [74]. A key result from scatter (X-Y) plots of air pollutants is the graphical interpretation of source compositions as 'limiting edge lines' [13,75], which define an upper (lower) edge line with the highest (lowest) Y/X ratio. Furthermore, the ratio PM2.5/PM<sup>10</sup> is different for combustion particles than for mechanically-generated particles [3]. The proposed rules and their rationale are detailed next.


in resolving intermittent sources with low contributions to ambient PM2.5 [12]. Nonetheless, they might indicate long-range, regional sources arriving to the monitoring site. Thus, these may be analyzed on their own using backward wind trajectories [77,78] to confirm their identity (natural or anthropogenic).


#### **3. Results**

#### *3.1. Results for Calama*

In this city in an arid environment, annual ambient PM2.5 concentrations are below 10 µg/m<sup>3</sup> , that is, ambient PM2.5 satisfies current WHO guidelines [79]. Thus, we focus on ambient PM<sup>10</sup> as the pollutant of concern, considering that fugitive sources, such as those windblown from desert surroundings and road dust, should be relevant. Both sources are also hard to estimate, given the intermittence of physical processes that lead to high wind speeds (gustiness conditions) in the former source, and the heterogenous processes responsible for surface dust loading on the city's streets in the latter source [51].

For this city, no RM results have been published. Nonetheless, available urban emission inventories for Calama [80] estimate that road dust (or traffic non-exhaust emission) is the dominant PM<sup>10</sup> emission source, being 88% of urban emissions (Supplementary Table S1). Therefore, traffic emissions dominate PM<sup>10</sup> emissions. However, no estimation of windblown dust is available for this study zone. We conduct the analysis for the monitoring station 'Centro' in Calama (see Table 1 for further details).

Figure 2 shows the outcome of the polarCluster routine when applied to ambient PM<sup>10</sup> data gathered from October 2012 through August 2020. The clusters associated with windblown dust emerge with high wind speeds and wind directions between 270◦ and 360◦ , and stay the same for solutions of CA with 5–9 clusters. The other clusters include low wind speed conditions and they split into several wind direction sectors as the number of clusters increases. For simplicity, we choose a solution with five clusters and explore the results using the simple rules presented in Section 2.3 to analyze the results.

Figure 3 shows the source contributions associated with the five-cluster solution for PM10. Clusters 1 and 2 are the dominant ones, followed by clusters 4, 3, and 5. The first analysis is performed looking at the time variability of this five-cluster solution. Figure S1 shows that clusters 3 and 5 present higher contributions between noon and 6 PM when wind speeds are higher in the region. Figure S2 shows scatter plots of ambient PM<sup>10</sup> versus wind speed. PM<sup>10</sup> concentrations from clusters 3 and 5 clearly increase with wind speed, which is distinctive of Aeolian emissions. By contrast, other clusters' contributions have little increase or decrease with wind speed (see R<sup>2</sup> values for instance). Cluster 3 comes from W and WNW directions and cluster 5 from NNW winds. Cluster 3 is characteristic of upslope, anabatic winds that develop all year long given the strong solar irradiation and extreme soil dryness [81]. Cluster 5 is characteristic of mountain synoptic conditions that happen in the austral fall

and winter seasons, when the Pacific subtropical high weakens, favoring a higher frequency of north winds. Nonetheless, both windborne dust sources are minor contributors to a long-term average PM10, as can be seen in Figure 3.

**Figure 2.** Cluster analysis source apportionment results for ambient PM<sup>10</sup> at Calama, 2012/10–2020/08, for 2–10 clusters solutions. The numbers stand for wind speed levels (m/s).

Figure 4 shows the PM<sup>10</sup> time variability for clusters 1, 2, and 4. It is clear from this figure that these three clusters correspond to traffic sources: their ambient contributions have distinctive peaks during morning and evening traffic rush hours, and they significantly decrease over the weekends. The same behavior can be seen in Figure 5, where the NOx time variability is plotted for clusters 1, 2, and 4. Figure 6 shows scatter plots of PM2.5 and PM<sup>10</sup> by cluster. In clusters 1, 2, and 4, the points lie between two limiting 'edge lines': one upper edge with high PM2.5/PM<sup>10</sup> ratio (traffic exhaust) and a lower edge with a low PM2.5/PM<sup>10</sup> ratio (lower that 0.1), characteristic of non-exhaust traffic emissions [46,51]. Thus, any point in the plot for those three clusters can be regarded as an air mass that arrives to the monitor site with a mixture of those two traffic emissions, as it is customarily presented in the receptor modeling literature [13,75].

**Figure 3.** Cluster analysis source apportionment results for ambient PM<sup>10</sup> at Calama, 2012/10–2020/08, for a 5-cluster solution.

**Figure 4.** Time variation plot for cluster analysis source apportionment results for ambient PM<sup>10</sup> at Calama, 2012/10–2020/08, for a five-cluster solution. Only clusters 1, 2, and 4 are plotted for clarity's sake.

**Figure 5.** Time variation plot for cluster analysis source apportionment results for ambient NOx at Calama, 2012/10–2020/08, for a five-cluster solution. Only clusters 1, 2, and 4 are plotted for clarity's sake.

**Figure 6.** PM-2.5–PM<sup>10</sup> scatter plot by cluster for a five-cluster solution for Calama.

#### *3.2. Results for Temuco*

In this city located in a wet temperate climate, RWB is the major source of ambient PM2.5. RM results show that, in the winter season, RWB contributions to ambient PM2.5 vary between 70% and 100% [60]. These RM results were found for a monitoring site located 700 m NW of the air quality monitoring site that we analyze here (denoted as LE by its Spanish name 'Las Encinas'). The data cover

the period from January 2009 through August 2020. In this case study, we use temperature instead of wind speed for conducting the CA. Had we chosen the traditional CA, we would have found difficulties in interpreting the resulting clusters because of an overlapping of high PM2.5 contributions at low wind speeds, which mixes traffic and RWB contributions (results not shown). By using temperature as an input variable, we are able to extract the RWB contribution and separate it from the traffic contribution. This approach works well precisely because the highest RWB impacts occur at nearly midnight in colder months when traffic sources are minimal and RWB emissions are the highest. Furthermore, since under those stable atmospheric conditions, wind direction is highly variable because of wind meandering effects, RWB contributions should come from all wind directions, defining a single cluster enclosing the origin of coordinates in the bivariate plot. This geometric feature helps in identifying RWB sources, even when they are not the dominant ones (see Section 3.3).

Figure 7 shows the cluster analysis results for two to eight clusters in Temuco, Las Encinas (LE) site. It can be clearly seen that a 'central cluster' develops for solutions with three or more clusters. For seven or more clusters, this central cluster splits itself into a pair of low-temperature and high-temperature clusters. They both correspond to RWB (results not shown here). For simplicity, we choose a three-cluster solution to identify the sources.

**Figure 7.** Cluster analysis source apportionment results for ambient PM2.5 at Temuco, 2009/01–2020/08, for 2–8 cluster solution.


Figure 8 shows the monthly source contributions of the three-cluster solution found in Temuco. Source 3 is the dominant one, followed by sources 1 and 2. Figure 9 shows the time variability of those three clusters. Clusters 1 and 2 show lower seasonality than cluster 3, which has dominant contributions near midnight. Figure S3 shows the pollution rose by the cluster. Cluster 3 has contributions for low wind speeds and from all wind directions, while clusters 1 and 2 are each associated with a set of narrow wind directions. Figure 10 displays PM2.5-PM<sup>10</sup> scatter plots. Clusters 1 and 2 have two clear limiting edge lines—characteristic of traffic sources—while cluster 3 points lie along a straight line with a high slope near 1. Figure S4 shows CO time variability by cluster. The cluster patterns mimic the ones already shown in Figure 9 for PM2.5. Figure 11 show scatter plots of PM2.5 against the temperature. Cluster 3 has a different behavior and higher PM2.5 concentrations than in clusters 1 and 2. Therefore, the CA results indicate that residential wood burning is the dominant source of ambient PM2.5, which is followed by traffic sources (clusters 1 and 2) with lower contributions.

We now turn to a quantitative comparison of the above CA results with RM results obtained for the period July–September 2014 at a monitoring site 700 m NW of the LE site. Table 2 shows a comparison of the average sources identified with CA and those resolved using the chemical mass balance RM (CMB8.2) with organic molecular markers [60]. The molecular markers measured in ambient PM2.5 samples were used to apportion organic carbon (OC) among different combustion sources (gasoline, diesel, wood, coal, natural gas) using Equation (1) with known source profiles {*fkj*}. Afterward, organic source contributions to PM2.5 mass were calculated from those source contributions to OC and specific OC/PM2.5 mass ratios for each source [82–85].

**Figure 8.** Cluster analysis source apportionment results for ambient PM2.5 at Temuco, 2009/01–2020/08, for a three-cluster solution.

**Figure 9.** Time variation plot for cluster analysis source apportionment results for ambient PM2.5 at Temuco, 2009/01–2020/08, for a three-cluster solution.

**Figure 10.** PM2.5–PM<sup>10</sup> scatter plot by cluster, for a three-cluster solution for Temuco.

**Figure 11.** PM2.5-temperature scatter plot by cluster for a three-cluster solution for Temuco.


**Table 2.** Comparison for cluster analysis (CA) and receptor model (RM) source apportionment estimates (±standard error) for Temuco, July–September 2014 <sup>1</sup> in (µg/m<sup>3</sup>

Average 42.5 45.2 2.2 2.2 <sup>1</sup> Supplementary material in Reference [60] provides further details and data. <sup>2</sup> Residential wood burning sources.

Week 7 23.1 ± 4.2 23.5 ± 5.0 5.1 ± 1.1 1.5 ± 0.1 Week 8 57.6 ± 9.0 48.6 ± 7.8 1.1 ± 0.3 2.8 ± 0.2

For the CA traffic sources, we compare them with the sum of diesel exhaust emissions, vegetative detritus (both resolved by CMB8.2), and dust (sum of Al, Si, Fe, Ca, and Ti oxides) to consider the non-exhaust contributions included in the CA (i.e., road dust). For the residential wood burning source, the comparison is more elaborated. We have found that, in Temuco, coal combustion is also used for space heating (identified using picene as organic tracer), so we have added RWB contributions (identified using levoglucosan as an organic tracer) to those from residential coal combustion since both processes happen under the very same environmental conditions. Furthermore, there is a substantial contribution of secondary organic aerosol (SOA) in Temuco, which is ascribed to the inefficient burning of wood. This combustion process is known to release semi-volatile organic compounds [53], which quickly oxidize and, thus, generate secondary organic aerosols [86]. The CMB8.2 RM cannot resolve secondary sources, so the unresolved organic carbon (OC) is denoted as 'Other OC'. This 'Other OC' is identified as SOA by its high correlation with water soluble organic carbon (WSOC), indicating a high degree of molecular oxygenation. Therefore, we added the corresponding SOA contribution to the RWB (and coal combustion contribution) to get the 'RWB RM' entry in Table 2.

We show in Table 2 the previously mentioned comparisons for the 8-week ambient monitoring campaign reported in Reference [60]. For every weekly sample, we computed the average CA contributions for both sources (Traffic: clusters 1 and 2, RWB: cluster 3), considering the same five sampling days per week as in the ambient measurement campaign. Standard errors for all estimates (from CA and RM results) were computed from error propagation. For both major sources that can be resolved with CA, the agreement is good with very similar winter average results. The zero value for one week in the traffic CA contribution is a result of trying to identify a weak signal in a data set dominated by a single source (RWB in this case). This is the usual outcome both in CA and RM analyses.

Therefore, the comparison results in Table 2 show that the proposed CA analysis is able to capture the major sources contributing to ambient PM2.5 in this case study of a zone dominated by residential wood burning emissions. For the (smaller) traffic contributions, the methodology has difficulties in capturing the temporal variation, even though, on average, the results are similar to the RM results.

#### *3.3. Results for Santiago*

In this case study of a large metropolitan area, we focus the analysis on a monitoring site located near the east edge of the city at a higher elevation in Santiago's basin (henceforth, denoted as LAC, which is short for its Spanish name 'Las Condes'). For that site, a previous RM study [61] has shown that traffic sources and residential wood burning are the dominant ones, although regional sources also contribute to ambient PM2.5. That study was conducted with another RM, Positive Matrix Factorization (PMF), using elemental concentrations in ambient PM2.5 daily samples as input data.

Figure 12 shows the results of CA for 2 to 10 candidate clusters, using temperature instead of wind speed as input variable, to resolve RWB contributions. A central cluster with the lowest ambient temperatures shows up for solutions with seven or more clusters, suggesting this is the RWB source. We choose an eight-cluster solution and we present these results here. Figure 13 shows the source contributions brought by this eight-cluster solution. Clusters 3 and 4 are the ones with the largest source contributions, followed by cluster 7 (that groups the lowest ambient temperatures) and clusters 1 and 6. The rest of the clusters are of minor relevance. Figure S5 shows that cluster 3 has predominant WSW directions while cluster 4 includes ENE and E directions. Figure 14 shows the PM2.5 time variability but only for the five major contributing clusters. It can be seen that cluster 4 has morning and evening peaks, coincident with rush hour traffic conditions. In cluster 3, its contribution increases from morning to early afternoon, indicating the arrival of traffic contributions from the city, brought by anabatic winds. This rise is followed by a decline of contributions later in the evening. This daylight behavior of anabatic winds and pollution transport toward the east side of the city has been recently measured and modeled for black carbon particles in Santiago [87]. Cluster 7 is the only one that does not decrease over weekends, unlike the other clusters shown in Figure 14, which decrease over weekends. This temporal pattern supports the identification of cluster 7 as the RWB sources.

Figure 15 shows PM2.5-PM<sup>10</sup> scatter plots by cluster, showing that, in cluster 7, the data have the highest slope of all, suggesting that this is the RWB source. For clusters 1–5 and 8, the scatter plots show that the respective data points have upper and lower limiting edge lines, as expected for traffic sources.

To better identify this eight-cluster solution, we present in Table 3 the monthly average PM2.5 contribution by cluster, for year 2004 for which we have RM results in this same monitoring site.


**Table 3.** Monthly source contributions to ambient PM2.5 by cluster for Santiago, year 2004 (µg/m<sup>3</sup> ).

**Figure 12.** Cluster analysis source apportionment results for ambient PM2.5 at Santiago, 2004/01–2012/06, for 2–10 cluster solutions.

**Figure 13.** Cluster analysis source apportionment results for ambient PM2.5 at Santiago, 2004/01–2012/06, for an eight-cluster solution.

**Figure 14.** Time variation plot for the cluster analysis source apportionment results for ambient PM2.5 at Santiago, 2004/01–2012/06, for an eight-cluster solution. For clarity's sake, only the five major clusters are included.

**Figure 15.** PM2.5–PM<sup>10</sup> scatter plot by cluster, for an eight-cluster solution for Santiago.

From Table 3, we can see that clusters 4, 7, and 8 increase their contributions during the fall and winter season, whereas clusters 1, 2, and 3 show minimum contributions in those colder seasons. The reason for this different behavior has to do with meteorological conditions. In the fall and winter seasons, subsidence conditions promote a low level of thermal inversion layers over Santiago's valley, leading to shallow planetary boundary layers (PBL) [88] and blocking transport of emissions from Santiago's lower valley (dominant wind direction for clusters 1, 2, and 3 as seen on Figure S5). This also explains why local sources (clusters 4, 7, and 8) increase their contributions during the fall and winter. This topography-induced effect has been reported before for total ambient PM concentrations [89] and, more recently, in the simulation of black carbon transport from Santiago toward the Andes mountains east [87].

Regarding cluster 6, it has a different time variability as compared with the rest of the resolved clusters with contributions peaking in the afternoon and increasing in the fall and winter seasons. Therefore, they do not come from Santiago's lower valley. RM results [61] indicate that regional anthropogenic sources contribute to ambient PM2.5 at the monitoring site, diagnosed by the presence of arsenic and sulfur in ambient PM2.5 samples. To check whether cluster 6 could represent those regional sources, Figure S6 shows the source contributions to ambient SO2. It can be seen that cluster 6 SO<sup>2</sup> contributions show up all year long, so they cannot come from Santiago's lower valley. Besides, Figure S5 shows that cluster 6 has W and WSW wind directions, suggesting that those regional sources are located west of Santiago. This result agrees with the location of regional sources of arsenic and sulfates identified in Reference [61] using backward trajectory analyses. Therefore, we conclude that cluster 6 can be identified as representative of regional sources of PM2.5.

In order to make comparisons between the CA results and the RM results, we need to consider the limitations of both analyses. The RM results were obtained using Positive Matrix Factorization software and elements as tracers. Organic tracers were not included and this is a limitation. For instance, more recent results [90,91] have shown that secondary organic aerosols are relevant in the spring and summer seasons in Santiago. This secondary PM2.5 was not resolved by the PMF solution including only elements. Another limitation of the PMF solution is the use of potassium as a tracer of wood burning. Potassium may be a reasonable tracer in fall and winter seasons, but not so specific in the spring and summer when soil dust becomes relevant in Santiago's semi-arid climate [70]. Because of

these issues, we have decided to compare CA and RM results only for the months from May through August 2004. Table 4 shows such a comparison. We have grouped clusters 3 and 4 as traffic sources, cluster 7 as the RWB source, and cluster 6 as regional sources. The smaller contributions from other clusters have not been considered, given the limitations of both CA and RM to resolve smaller (or intermittent) sources. Since the filter-based RM results consider only a subset of days per month, these specific days have been extracted from the CA results to compute comparable averages.


**Table 4.** Ambient PM2.5 source apportionment (±standard error), Santiago, May–August 2004 (µg/m<sup>3</sup> ).

From the results in Table 4, it can be seen that CA tends to overestimate the RM result for traffic contributions (which includes traffic and soil dust contributions) but CA results for RWB and regional sources contributions that are below those estimated by the RM in the very same monitoring site. In fact, the sum of CA estimates in Table 3 has an average of 33.3 (µg/m<sup>3</sup> ), while the corresponding average of RM results is 41.3 (µg/m<sup>3</sup> ). The reason for this discrepancy is ascribed to the different PM2.5 measurement techniques: filter samples were taken in low-volume dichotomous samplers (Andersen Instrument, Inc, Smyrna, GA, USA, 15 L/min) while continuous measurements were made with a Tapered Element Oscillating Microbalance equipment (TEOM, Rupprecht & Patashnick, MA, USA). The latter instrument is known to present negative artifacts (i.e., underestimation of PM2.5 concentrations) due to partial volatilization of sampled PM2.5 in the TEOM's heating inlet used to dry the samples [92]. We think this instrument artifact explains why filter-based RM results are higher than the CA results reported here. This effect explains the lower contributions found for RWB in the CA because RWB sources have a larger fraction of volatile compounds emitted, as compared with other combustion sources [53]. We do not have such an issue in the case of Temuco because, in that monitoring site, a beta-attenuation monitor (BAM, MET-ONE 1020, Met One Instruments Inc., Grants Pass, OR, USA) has been used to measure PM2.5 and PM10.

#### **4. Conclusions**

We have proposed that a CA approach, applied to ambient air pollution and meteorological data, provides a source apportionment of ambient PM2.5 and PM<sup>10</sup> on an hourly basis. In order to achieve this result, we need to compare the outcome of the CA with RM results for the same monitoring site—or using an emission inventory in case no RM result is available—to identify the major sources (clusters) at play on a given zone. We have shown that our rule-based CA works for three different case studies in Chile: a city in a warm, desert region (Calama), another one dominated by RWB pollution in a cold, wet region (Temuco), and a large metropolitan area in a semi-arid region (Santiago).

In the case of Calama, the CA for ambient PM<sup>10</sup> is able to resolve local sources (traffic) and windblown dust coming from the nearby desert environment. Both are fugitive sources, which are difficult to estimate, because of the amount of information required on each case, such as particle size distributions. CA results indicate that traffic sources dominate ambient PM<sup>10</sup> concentrations, as suggested by the emission inventory of PM<sup>10</sup> sources for that city. CA results show that, over the long-term, outbursts of windblown dust are not significant within the city. This is relevant for policy purposes. The evolution of traffic contributions in Calama (Figure 3) suggests that the ongoing street sweeping program has been a successful one with a clear decreasing trend. This is an example of how a sector-specific regulation in the city may be evaluated. Furthermore, within the traffic source, ambient PM<sup>10</sup> data can be regarded as a combination of exhaust and non-exhaust emissions, providing

additional insights such as whether new exhaust emission standards for motor vehicles have curbed down exhaust emissions.

In the case of Temuco, with dominant RWB contributions to ambient PM2.5, the use of ambient temperature instead of wind speed improves the apportionment of all major sources of ambient PM2.5. RWB contributions show up in the bivariate polar plots as a 'central' cluster because, under stable atmospheric conditions (with lowest ambient temperatures and wind speeds), wind direction is highly variable. Hence, the monitor site will sample air masses from all wind directions. CA source contributions for RWB sources are statistically comparable with RM results obtained for RWB sources in a short-term campaign in 2014.

For the large metropolitan area of Santiago, the CA methodology resolved RWB, traffic, and regional sources of ambient PM2.5. The CA analysis again required using ambient temperature to resolve the RWB contribution. The seasonality of the resolved clusters showed a distinctive effect brought by topography: Santiago's plume does not reach the eastern side of the city in fall and winter seasons due to low PBL depths in colder months. This geographic condition provided an additional criterium to identify local and non-local traffic sources. However, the RM resolved source contributions (short-term campaign in 2004) were higher than the CA resolved counterparts. This is a result of an underestimation of ambient PM2.5 monitoring brought by volatilization losses in the continuous TEOM PM2.5 monitor.

The rule-based CA presented here generates added value to existing ambient air pollution databases. Regarding RM results to complement the CA analysis, methods that use specific source tracers (like organic molecular markers) are preferable.

The rule-based CA results may be applied to the following analyses:


The proposed rule-based CA analysis has limitations, though. The methodology works well for resolving the larger sources at play in a city, but smaller sources remain a challenge.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1660-4601/17/22/8455/s1. Table S1: Emission inventory for Calama, year 2016 [ton/year]. Figure S1: PM<sup>10</sup> Time variability results for a five-cluster solution for Calama. Figure S2: PM10-wind speed scatter plot by cluster for a 5-cluster solution for Calama. Figure S3: Pollution rose results by cluster for a 3-cluster solution for Temuco. Figure S4: Time variability of CO by cluster for a 3-cluster solution for Temuco. Figure S5: Pollution rose by cluster for an 8-cluster solution for Santiago. Figure S6: Source apportionment for SO<sup>2</sup> for an 8-cluster solution for Santiago.

**Author Contributions:** Conceptualization, H.J. Methodology, H.J. Validation, A.M.V. Writing—original draft preparation, H.J. Writing—review and editing, H.J. and A.M.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Chile's ANID, grant number FONDAP/15110020 (CEDEUS).

**Acknowledgments:** We acknowledge Gabriela Mancheno for her help in drawing some figures. We also thank two anonymous reviewers for their helpful comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Analysis of Alpha Activity Levels and Dependence on Meteorological Factors over a Complex Terrain in Northern Iberian Peninsula (2014–2018)**

#### **Miguel Ángel Hernández-Ceballos <sup>1</sup> , Fernando Legarda <sup>2</sup> and Natalia Alegría 2,\***


Received: 18 September 2020; Accepted: 28 October 2020; Published: 29 October 2020

**Abstract:** Alpha ambient concentrations in ground-level air were measured weekly in Bilbao (northern Spain) by collecting aerosols in filters between 2014 and 2018. Over this period, the alpha activity concentrations in the aerosol's samples range from 13.9 µBq/m<sup>3</sup> to 246.5 µBq/m<sup>3</sup> , with a mean of 66.49 ± 39.33 µBq/m<sup>3</sup> . The inter-annual and intra-annual (seasonal and monthly) variations are analyzed, with the highest activity in autumn months and the lowest one in winter months. Special attention has been paid to alpha peak concentrations (weekly concentrations above the 90th percentile) and its relationship with regional meteorological scenarios by means of air mass trajectories and local meteorological parameters. The meteorological analysis of these high alpha concentrations has revealed two airflow patterns-one from the south with land origin and one from the north with maritime origin-mainly associated with these alpha peak concentrations. Surface winds during representative periods of both airflow patterns are also analyzed in combination with <sup>222</sup>Rn concentrations, which demonstrated the different daily evolution associated with each airflow pattern. The present results are relevant in understanding trends and meteorological factors affecting alpha activity concentrations in this area, and hence, to control potential atmospheric environmental releases and ensure the environmental and public health.

**Keywords:** gross alpha activity; northern Iberian Peninsula; radon; airflow patterns; surface winds

### **1. Introduction**

Radioactivity is naturally present in the air due to cosmogenic radionuclides (14C, 7Be, and 3H), cosmic radiation, and natural radionuclides present in all rocks and soils (238U and progenies, <sup>232</sup>Th and progenies, and <sup>40</sup>K), and its natural background is increased by anthropogenic sources, such as nuclear/radiological accident, nuclear weapon tests, or the production and application of phosphorus fertilizers [1].

The impact of radioactivity on the environment and population health makes the control of radioactivity levels necessary. In this sense, alpha activity concentration is usually used as a control indicator of the radioactivity levels in air, and it is a common measurement of detecting activity peaks associated with radioactive released to the atmosphere [2]. This activity concentration can be measured counting the number of alpha emissions produced per unit of time and air volume (gross alpha), which are determined from the sampling of atmospheric aerosols. The determination and formation of radioactive aerosols in the atmosphere air are analyzed in detail in Papastefanou [3]. In case of alpha activity concentration in aerosols, their origin can be natural due to radon daughters [3].

The long-term monitoring of gross alpha activity concentrations enables us to study trends and periodic variations and, in combination with meteorological parameters, to identify the main meteorological causes influencing the spatial and temporal variability of the environmental radioactivity levels [4]. The seasonal variation and its link with meteorological parameters such as rainfall, temperature, relative humidity, pressure, visibility, wind speed, and wind direction, which influence the concentration of aerosols radiotracers in the atmosphere, has been analyzed in several sites in Spain. Dueñas [5,6] and Piñero-García [7] in the south, García-Talavera [8] in the west, and Saéz-Muñoz [9] in the east provide information about the spatial distribution and the meteorological conditions driven gross alpha concentrations in surface air in Spain. However, in the northern part of Spain, there are completely different meteorological and climatic conditions, and to the authors' knowledge, there are no studies analyzing alpha concentrations.

The direct relevance of alpha activity on the effective dose causes the need of analyzing and characterizing these concentrations. Depending on the aerosol's size, in the respiration process, some aerosol settling occurs in the alveolar spaces [10], and if there are some radionuclides on them, internal irradiation is produced [11]. Due to the inhalation of the short-lived radon decay products, this radionuclide delivers the largest proportion of the natural radiation dose to humans [3]. The alpha energy released by Rn progeny is close to 35 GeV/Bq.

In this paper, we analyze the temporal variation of the radioactive aerosol concentrations by measuring weekly the gross alpha activity over five years (2014–2018) at Bilbao (northern Spain), and due to the potential health impact of alpha peak concentrations, we investigated the meteorological factors controlling them. To this end, the link between local/regional meteorological condition and alpha peak concentrations are examined by means of the calculation of backward trajectories (regional scale) and the temporal variability of surface winds and <sup>222</sup>Rn measured (local scale). This last analysis is based on the fact that, in normal situations, gross alpha origin in air is mainly explained due to the presence of long-lived daughters of gaseous <sup>222</sup>Rn that are attached to aerosols after cluster formation [3].

The following research issues are addressed in the present study:


The structure of this paper is the following: Section 2 describes the measurement place, the radioactivity measurements, and the meteorological parameters and data used for the study. Section 3 presents, in a first step, the statistical analysis of the alpha activity concentration, while in a second step, the correlation of high alpha concentrations with regional and local meteorological parameters is investigated. To finalize this results section, the dependence of alpha activity concentrations, by means of <sup>222</sup>Rn values, on the surface wind conditions is analyzed. In Section 4, the conclusions are shown.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Bilbao City (43.26 N, −2.94 W) is placed in the Basque Country region (northern Spain) (Figure 1). Bilbao is about 16 km away from the sea, in a mountainous coastal area at the Gulf of Biscay, and in the narrow valley of the Nervion river. The mountains surrounding Bilbao have low altitude (300–800 m), although 40 km away there are mountains with altitudes of 1500 m.

The climate in Bilbao is warm and temperate, with a significant amount of rainfall during the year. Air masses mainly arrive from the west and north to this area. In the period covered by the present study (2014–2018), the temperature averaged 16.0 ◦C, while the precipitation was about 980 mm per year.

**Figure 1.** Location of the sampling site in the northern area of the Iberian Peninsula.

### *2.2. Alpha Activity Concentration and Radon Measurements.*

The concentration of radionuclides in soil in this area is very low based on the European Digital Atlas of Natural Radiation [12] (for <sup>40</sup>K, uranium, thorium, etc.).

− The alpha activity concentrations were measured in a monitoring station sited in Bilbao city (43.26 N, −2.94 W) at 34 m above sea level, on the roof of the Faculty of Engineering (Figure 2). −

**Figure 2.** The aerosol sampler, the meteorological station and the radon station used in the present analysis.

Atmospheric airborne aerosols were collected using a low volume sampler whose nominal flow rate is 30 L/min. It is important to remark that the aerosol sampling takes seven days. In a week the total volume circulated is about 300 m<sup>3</sup> . This sampler uses a cellulose nitrate membrane filter (diameter 47 mm and 0.8 um pore size) which was replaced weekly. The filter is taken from the sampler and placed in a dryer for six days to allow thoron and radon short-lived alpha progeny to decay, and, then, it is measured in a gas flow proportional counter. The gas flow proportional counter used is LB 770 10-Channel α-β Low-Level Counter by Berthold Technologies GmbH & Co. KG trademark. This device is equipped with 10 detectors. The efficiency of the proportional counter calibrated with <sup>241</sup>Am is close to 20%. A calibration source was prepared ad hoc with the same type of filter as used. The counting time was 6 × 10<sup>4</sup> seconds, and the detection limits obtained were between 9.71 µBq/m<sup>3</sup> and 50.68 µBq/m<sup>3</sup> and uncertainties, with a coverage factor k = 2, were between 14 ± 6.7 µBq/m<sup>3</sup> .

A high flow sampler and a radiological station are located close to this low flow sampler (Figure 2). The radiological station provides hourly radon activity concentrations in the atmosphere, which have also been used in this study.

#### *2.3. Meteorology and Backward Trajectories*

Ten-minute intervals of the following surface meteorological parameters, which were measured in the weather station placed close to the sampling station (Figure 2), were also used in the present analysis to characterize the impact of local meteorology on alpha activity concentrations: air temperature (◦C), relative humidity (%), pressure (mbar), wind direction (◦ ), wind speed (m/s), and precipitation (mm).

Backward trajectories were calculated and used in the present study. Four kinematic three-dimensional backward trajectories per day during specific sampling periods between 2014 and 2018, computed at 00:00, 06:00, 12:00, and 18:00 UTC, with a run time of 96 hours and with a final height of 100 m above ground level, were calculated by using the hybrid single particle Lagrangian integrated trajectory (HYSPLIT) model [13]. The Global Data Assimilation System National Centers for Environmental Prediction (GDAS-NCEP) model, maintained by the National Oceanic and Atmospheric Administration Air Resources Laboratory (NOAA/ARL) were used as meteorological input.

Working with a large number of trajectories requires the application of cluster analysis [14] to classify them according to their pathways, and hence, to identify the corresponding airflow patterns. The present study uses the cluster methodology implemented in the HYSPLIT model and described in Stunder [15]. This multivariate statistical technique groups trajectories with similar direction and lengths by using the Ward's minimum variance clustering method to minimize the dissimilarity in the trajectories within each cluster while maximizing the dissimilarity of different clusters.

#### **3. Results**

#### *3.1. Alpha Activity Concentration Characterization: Time Series*

The alpha activity concentration values, measured weekly from 2014 to 2018, are plotted in chronological order (Figure 3a). In total, there were 258 data, and they were available in each year, thus guaranteeing the largest statistical sample. The activity concentrations in the aerosol's samples range from 246.5 µBq/m<sup>3</sup> to 13.97 µBq/m<sup>3</sup> , with a mean of 66.49 ± 39.33 µBq/m<sup>3</sup> , where the uncertainty is given as the sample standard deviation. This Figure 3a suggests a cyclical variability every year, with low values in the first half of the year followed by an increase and low values at the end of the year. This behavior is only broken in 2015. These alpha activity concentrations are also assessed in a box-and-whisker plot (Figure 3b), which shows that there is a symmetrical distribution of the values, i.e., differences between P75 and P25 with the median (P50) are similar. The average is greater than the P50, which denotes the prevalence of low alpha values, although with large impact of occasional high ones. Using the Grubbs test, the existence of outliers among the registered values has been assessed, concluding that there are no outliers in the present set of data.

Alpha measurements in Figure 3a constitute a time series due to the existence of correlation delays of up to 12 weeks between data [16,17]. The periodogram calculated for this period is shown in Figure 4. A periodogram is a helpful tool used to identify the dominant cyclical behavior of a time series. In the Figure 4, the *X*-axis represents frequency and the Y-axis represents the spectrum. Four main peaks, i.e., periodicities of alpha concentration records, are identified in Figure 4. The first one reveals a one-year period variation (frequency 1.93 × 10−<sup>2</sup> week−<sup>1</sup> i.e., 52 weeks), while the second, third and

a)

)

100 150

250 300

−

fourth peaks indicate intra-annual variability (frequencies 3.9 × 10−<sup>2</sup> , 6.9 × 10−<sup>2</sup> and 1.6 × 10−<sup>1</sup> week−<sup>1</sup> i.e., 25, 12, and six weeks). These results highly agree with the components of the <sup>7</sup>Be time series identified and analyzed in previous studies [4,18] 0 50 01/01/2014 01/07/2014 01/01/2015 01/07/2015 01/01/2016 01/07/2016 01/01/2017 01/07/2017 01/01/2018 01/07/2018 01/01/2019 Activity concentration (µBq/m3

**Figure 3.** (**a**) Temporal variation and (**b**) box-and-whisker plot of the alpha activity concentrations from January 2014 to December 2018 at Bilbao sampling station. In **b**, the center of the box denotes the 50% (P50), and the bottom and top of the box correspond to the 25% (P25) and 75% (P75) values, respectively. The square indicates the P90 and the circle the P10 values, while the extremes of the box represent the P95 and P5 values. The cross in the box is the average value. − − − − − −

**Figure 4.** Periodogram for Bilbao time series of alpha activity concentrations from 2014 to 2018.

Figure 5 displays the box-and-whisker plots of alpha activity concentrations over the years. This Figure 5 shows that the average (square) and the variability of alpha values (P75–P25) are similar over the years, which can be associated with the minimal change in the source term of alpha activity concentrations during this five-year period. The averages are always higher than P50. The average ranges from 60.2 µBq/m<sup>3</sup> in 2018 to 72.9 µBq/m<sup>3</sup> in 2015, and the 90th percentile (P90), which is usually

taken as reference of maximum values [19], to avoid the impact of possible outliers, ranges from 108.1 µBq/m<sup>3</sup> in 2014 to 117.2 µBq/m<sup>3</sup> in 2015.

− **Figure 5.** Box-and-whisker plot of alpha activity concentrations over the years (2014–2018) at Bilbao. The meaning of the box plot is the same of Figure 3b.

The analysis of the intra-annual variability, considering winter (December, January, and February), spring (March, April, and May), summer (June, July, and August), and Autumn (September, October, and November), reports the maximum average value in autumn (Figure 6a), concretely in October (110 µBq/m<sup>3</sup> ) (Figure 6b), while the minimum average value is in winter (Figure 6a), with the monthly minimum in February (40 µBq/m<sup>3</sup> ) (Figure 6b). These seasonal and monthly variability do not agree with previous studies carried out in Spain [5,8,9], in which the activity maxima are registered in summer. This difference would point out differences in source terms and in the meteorological conditions driven alpha activity concentrations during the years among these sites, and hence, the need to analyze the link between alpha and meteorological parameters.

**Figure 6.** (**a**) Stational variability and (**b**) monthly average of alpha activity concentrations at Bilbao from 2014 to 2018.

#### *3.2. Meteorological Characterization of the Alpha Peak Activity Concentrations*

−

−

− This section presents the results of studying the impact of meteorological parameters on the highest alpha activity concentrations registered in Bilbao for the five years analyzed (2014–2018).

−

The reason for taking the highest alpha activity concentrations as reference is due to health hazards associated with the presence of radionuclides in air [20].

High alpha surface concentrations are identified as values exceeding the P90 calculated over 2014–2018 (117.1 µBq/m3 in Figure 3b). P90 has been previously used in other research as the extreme criterion [18–21]. The average of the 26 high alpha values identified during the period 2014–2018 above P90 is 152.0 ± 6.6 µBq/m<sup>3</sup> , while the highest alpha value is 246.5 µBq/m<sup>3</sup> . Seasonal distribution of these alpha peak concentrations is shown in Figure 7, as well as the corresponding seasonal mean. This Figure 7 reports the largest occurrence of these high alpha concentrations in autumn (14 out of 26 periods, 53%), and the minimum occurrence in summer (three out of 26 periods, 11.5%). However, the maximum seasonal average of these peak concentrations is registered in spring (187.9 µBq/m<sup>3</sup> ), followed by summer (166.7 µBq/m<sup>3</sup> ), and lower and similar values in winter (142.4 µBq/m<sup>3</sup> ) and autumn (144.8 µBq/m<sup>3</sup> ). This seasonal variability of alpha peak concentrations, mainly explained by differences in atmospheric factors, is different to Figure 6a but similar to the seasonal variations shown in previous articles carried out in Spain, and, at the same time, it provides insight into the large impact of occasional high alpha values during the warm (spring/summer) period of the year.

− **Figure 7.** Averaged concentration in each season considering only those sampling periods registering alpha concentrations above P90 during 2014–2018 (in brackets below the season).

The evaluation of correlation between weekly alpha peak activity concentrations and arithmetic values of local meteorological parameters (temperature, relative humidity, atmospheric pressure, and wind speed and direction) during each sampling period is performed at the significance level of 0.05. In the case of precipitation, the total amount during each sampling period is used. The correlations found are around zero, i.e., there is no relationship between the two variables used. The highest is −0.22 with precipitation, which, being negative, indicates that precipitation tends to reduce alpha activity at the site as a consequence of scavenging of airborne radionuclides by rainfall [22]. These results, with none of them significant at 0.05, indicate the weakness of the relationship between alpha peak concentrations and each individual meteorological variable at this weekly temporal scale level. Therefore, other atmospheric processes should be included in the present analysis aiming at understanding the meteorological reasons of these high alpha surface concentrations in Bilbao.

Due to the impact that airflow patterns have over the surface weather conditions and aerosols temporal distribution in a given area and over a period of time [23], backward trajectories in Bilbao are calculated and are used to determine the influence of regional meteorological patterns on alpha peak concentrations. Figure 8a shows the mean trajectories (cluster) for Bilbao at 100 m calculated considering all the 26 sampling periods above P90. Seven airflow patterns are identified by applying the cluster methodology explained in Stunder [15]. Taking into account the pathways followed by each airflow pattern as well as the frequency (number in brackets) of each one, the arrival of slow

−

**Activity concentration (µBq/m3)**

continental southern (27%, cluster 4) and maritime northern (28%, cluster 1) air masses dominates periods with high alpha concentrations against the arrival of fast northern (7%, cluster 6) and different Atlantic advections from the northwest (24%, clusters 2 and 3) and from the west (14%, clusters 5 and 7) (Figure 8a).

−

Winter Spring Summer Autumn

(5) (4) (3) (14)

− **Figure 8.** Back-trajectory cluster centers (centroids) obtained at 100 m agl for periods with alpha activity concentrations (**a**) above the 90th percentile and (**b**) below the 10th percentile during the 2014–2018 period at Bilbao. The left numbers in the centroids are an identification number of the centroid and the right numbers (in brackets) are the percentage of complete trajectories occurring in that cluster.

To establish the link of these airflow patterns with alpha peak concentrations, we have also calculated the airflow patterns associated with the lowest alpha concentrations registered in Bilbao (26 periods in which alpha concentrations were below the 10th percentile) (Figure 8b). The comparison of both (Figure 8) points out similarities between airflow patterns, but with notable differences in frequencies (number in brackets) in those similar. The lowest alpha concentrations highlights the higher occurrence of western airflows, and air masses from the south and north, with slow displacement (Cluster 1 and Cluster 4 in Figure 8a), are not identified. Therefore, both Figures mainly suggest that the arrival of continental southern and maritime northern airflows is mainly associated with measurements of high alpha concentrations in this area.

To study the relationships between airflow patterns (Figure 8a) and alpha peak concentrations, and with the purpose to not consider the same sampling period within two different airflow patterns, which would create misunderstanding in the results, we have identified persistent sampling periods, i.e., periods in which most of the trajectories (above 75% of the calculated trajectories in each period) are grouped into the same cluster. Ten out of 26 sampling periods are identified as persistent ones: six within cluster 1, one within cluster 2, and three within cluster 4. Figure 9 shows the summary of alpha peak concentration measured and the seasonal distribution of these persistent sampling periods. They occurred in autumn/winter seasons in cluster 1 and cluster 4, while the only persistent period corresponding to cluster 2 took place in spring. This result, hence, reports that sampling periods characterized by prevailing southern and northern airflows seems to drive the occurrence of alpha peak activity concentrations in this area. The fact of mainly being located in autumn/winter, which is the wettest period in Bilbao (see Section 2.1), would not be in line with previous studies [6], where alpha activity concentrations decrease with precipitation. However, and due to the geographic location of the region, winds coming from the south descending from the Cantabrian Mountains became warm and dry because the Foehn effect [24], and even within the wettest period, under the arrival of southerly winds, the chances of rain are nearly zero in the region even.

**Figure 9.** Statistical values in air for alpha peak concentrations during the persistent sampling periods at Bilbao and the corresponding seasonal distribution of the periods. The uncertainties of the means are indicated as the standard deviation of the average Sx/N1/2, where N is the number of persistent periods and the standard deviation. Clusters refer to Figure 8a.

#### *3.3. Case Studies*

Case studies are needed to evaluate on a real basis the potential influence of each airflow pattern and the associated meteorological conditions on the alpha peak concentrations in Bilbao. So, we have analyzed representative periods of the two main meteorological scenarios driving these concentrations, i.e., the arrival of continental southern and maritime northern air masses. To analyze the influence of each airflow pattern on alpha activity concentrations, Figure 10 displays surface winds (blowing from) during one persistent sampling period (02−08/10/2018) under northerly flows, and from 9 December 2015 to 14 December 2015 under southerly flows. Hourly evolution of <sup>222</sup>Rn is also shown in both periods to represent the impact of surface winds on activity concentrations.

**Figure 10.** Daily variation of 222Rn activity concentration, wind speed and direction, and total amount of precipitation at Bilbao between (**a**) 2 October 2018 and 8 October 2018 (northerly flows-cluster 1) and (**b**) 9 December 2015 and 14 December 2015 (southerly flows-cluster 4). Note the differences in scales.

The alpha activity concentration recorded in the first sampling period is 122.5 µBq/m<sup>3</sup> . According to the wind direction variability, two different periods can be identified in these seven days of sampling (Figure 10a). This fact allows considering this sampling period as suitable to analyze differences in surface winds under the same synoptic meteorological scenario, and its relationship with surface activity concentrations. In addition, at the end of this sampling period was constantly raining.

The first four days, the synoptic conditions favored the progressive development of a daily evolution cycle of surface winds, blowing from the south during the night, followed by the clockwise swing to northwesterly winds during the day. Wind speed was less than 2 m/s. This daily evolution of winds at surface level agrees with the development of sea/land breezes in this area [25], which are more common during spring and summer. Under this surface condition, and considering the marked daily cycle of <sup>222</sup>Rn [26], it is suggested an accumulation process, which progressively raises the concentration of particles close to the surface in this area, i.e., it traps the particles near the surface and limits a sufficient mixing with the air above [27]. The limitation of the dispersion processes would thus favor the increase in the alpha activity concentration in this area. This daily cycle is broken at the end of the sampling period, with a sudden change in surface winds, which started blowing constantly from the northwest-north range and with wind speeds higher than in the previous days. Under these new conditions, the dispersion processes are favored, and together with the sporadic rainfall during these two days, cause a decrease in activity concentrations in the area, as it well observed in the daily evolution of <sup>222</sup>Rn concentrations.

The alpha activity concentration recorded in the second sampling period is 149.6 µBq/m<sup>3</sup> , and there is no rainfall measured in these six days. The synoptic meteorological conditions allowed the constant arrival of surface winds from the south. In this period, changes in wind speed allow identifying two different surface dynamics. In the first five days, the southerly pattern is developed with similar and low wind speeds, which favored similar daily activity concentrations but limited the accumulation processes, not registering high <sup>222</sup>Rn activity concentrations (comparing with Figure 10a, the radon concentrations differ by about two orders of magnitude). In the last day, the increase in wind speed cause a cleaning effect by enhancing the dispersive effect in the surface atmospheric layers, as can be seen in Figure 10b following the evolution of <sup>222</sup>Rn activity concentrations.

#### **4. Conclusions**

Alpha activity concentration can be usually used as a control indicator of the radioactivity levels in air, and so a low-flow station is located in Bilbao, northern Spain. Here, we have analyzed 258 data registered from 2014 to 2018. These alpha activity concentration data can be considered as a time series after autocorrelation function confirms its behavior and periodogram shows an annual periodicity and a significant seasonal and monthly variation of alpha activity concentrations.

In order to explain the alpha peak concentrations, the influence of the airflow patterns and local meteorological parameters has been studied. Precipitation presented a high negative correlation with alpha values, indicating removal process of airborne and the decrease in concentrations. Slow air masses from the south (continental displacement) and from the north (Atlantic displacement) have been found as mainly responsible for the alpha peak concentrations in this area. Under the arrival of northern airflows, the wind directions at the surface levels displayed the development of sea-land breezes, with northern and southern winds, while under the arrival of continental airflows, the arrival of southern winds at surface levels dominates. Analyzing the hourly evolution of <sup>222</sup>Rn together with surface winds, we have identified how the first meteorological scenario favors the accumulation processes and hence drastically increases activities, while the second scenario keeps activity concentrations similar but with lower values than in the previous scenario during the sampling period.

**Author Contributions:** Conceptualization, M.Á.H.-C. and N.A.; methodology, M.Á.H.-C. and N.A.; software, M.Á.H.-C. and N.A.; formal analysis, M.Á.H.-C. and N.A.; writing—original draft preparation, M.Á.H.-C. and N.A.; writing—review and editing, M.Á.H., N.A. and F.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors gratefully acknowledge the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT transport and dispersion model and/or READY website (https://www.ready.noaa.gov) used in this publication.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

International Journal of *Environmental Research and Public Health*

## *Article* **Role of Meteorological Parameters in the Diurnal and Seasonal Variation of NO<sup>2</sup> in a Romanian Urban Environment**

## **Mirela Voiculescu , Daniel-Eduard Constantin \* , Simona Condurache-Bota , Valentina Călmuc, Adrian Ros,u and Carmelia Mariana Dragomir Bălănică**

Faculty of Sciences and Environment, European Center of Excellence for the Environment, "Dunărea de Jos" University of Galat,i, 800008 Galat,i, Romania; Mirela.Voiculescu@ugal.ro (M.V.); scondurache@ugal.ro (S.C.-B.); valentina.calmuc@ugal.ro (V.C.); adrian.rosu@ugal.ro (A.R.); carmelia.dragomir@ugal.ro (C.M.D.B.)

**\*** Correspondence: Daniel.Constantin@ugal.ro

Received: 16 July 2020; Accepted: 24 August 2020; Published: 27 August 2020

**Abstract:** The main purpose of this study was to investigate whether meteorological parameters (temperature, relative humidity, direct radiation) play an important role in modifying the NO<sup>2</sup> concentration in an urban environment. The diurnal and seasonal variation recorded at a NO<sup>2</sup> traffic station was analyzed, based on data collected in situ in a Romanian city, Braila (45.26◦ N, 27.95◦ E), during 2009–2014. The NO<sup>2</sup> atmospheric content close to the ground had, in general, a summer minimum and a late autumn/winter maximum for most years. Two diurnal peaks were observed, regardless of the season, which were more evident during cold months. Traffic is an important contributor to the NO<sup>2</sup> atmospheric pollution during daytime hours. The variability of in situ measurements of NO<sup>2</sup> concentration compared relatively well with space-based observations of the NO<sup>2</sup> vertical column by the Ozone Monitoring Instrument (OMI) satellite for most of the period under scrutiny. Data for daytime and nighttime (when the traffic is reduced) were analyzed separately, in the attempt to isolate meteorological effects. Meteorological parameters are not fully independent and we used partial correlation analysis to check whether the relationships with one parameter may be induced by another. The correlation between NO<sup>2</sup> and temperature was not coherent. Relative humidity and solar radiation seemed to play a role in shaping the NO<sup>2</sup> concentration, regardless of the time of day, and these relationships were only partially interconnected.

**Keywords:** nitrogen dioxide; in situ urban concentrations; meteorological measurements; NO<sup>2</sup> variation; partial correlation

#### **1. Introduction**

Monitoring atmospheric pollution and the possibility to predict its evolution are of high interest. One major atmospheric pollutant is nitrogen dioxide (NO2), which may cause many health problems [1]. The NO<sup>2</sup> gas has a reddish-brown color, is nonflammable, and has a detectable, pungent odor, perceptible from concentrations of approximately 190 µg/m<sup>3</sup> [2]. Nitrogen dioxide reacts with the hydroxyl radical (-OH) in the atmosphere, forming the highly-corrosive nitric acid, but it can also form toxic organic nitrates. Nitrogen oxides known as NO<sup>x</sup> (i.e., NO + NO2) are involved in the formation of tropospheric ozone and smog, mediated by light through photolysis. Due to the significant environmental impact, the NO<sup>x</sup> compounds are strictly monitored and legislative limits were set at EU and national levels [2]. Part of the NO<sup>2</sup> molecules in the atmosphere are of primary nature (i.e., directly emitted), while most NO<sup>2</sup> results from nitrogen monoxide (NO) [1,3]. The latter is produced via natural and anthropogenic processes. About one-fifth of NO<sup>x</sup> is released in the atmosphere during

thunderstorms, through lightning, in the upper troposphere [3]. Other natural NO<sup>x</sup> sources include volcanic activity emissions, bacteria decay, soil processes and stratospheric sources [1], while the rest is the result of fossil fuel combustion processes: fires, transportation, industrial—such as power generation, nitric acid manufacture, welding processes, the use of explosives and iron and steel industry and refineries [1,2]. According to [4], about 65% of NO<sup>x</sup> emissions are of anthropogenic origin. The NO<sup>2</sup> reaction with the atmospheric hydroxyl radical (-OH), and the NO reaction with hydroperoxyl (HO2), followed by the formation of nitric acid and its wet deposition are important NO sinks (almost 60%), according to [4]. The contribution of road traffic to nitrogen oxides (NOx) emissions has been described by various studies [5–11]. Some authors assume that about half of the anthropogenic NO<sup>x</sup> comes from traffic [12]. Although several measures of emission control were implemented by various European Commission (EC) regulations, emissions of nitrogen oxides from road traffic have not significantly decreased across Europe; sometimes, NO<sup>2</sup> concentrations have even increased [10,11]. Some recent papers showed that the real reduction of nitrogen oxides and other pollutants following Euro 5 and 6a/b introduction was not as important as expected [10].

The concentration of NO<sup>2</sup> (and, in general, of any atmospheric pollutants) in an urban area and its variation are influenced by the industrialization level, the number of inhabitants and traffic density, but also by the topography of the area and meteorological parameters [5,9,11]. It is expected that meteorological factors would affect both natural and anthropogenic emissions of NO<sup>2</sup> (e.g., the reduced solar radiation during wintertime will impact the rate of NO<sup>2</sup> photochemical reactions and the reduced temperatures during the cold season will increase NO<sup>2</sup> emissions associated with heating). Moreover, meteorological parameters are interlinked in various ways: the relative humidity decreases with increasing temperature, the effective and saturation air vapor pressures increase with atmospheric temperature, etc. Thus, it is difficult to assess the effect of each separate meteorological factor on NO<sup>2</sup> pollution. None of these influences were yet quantified in a reliable way and a clear relationship is far from being established. Studies of the link between meteorological factors and atmospheric pollutants, in general, or NO2, in particular, are rather scarce, inconsistent and strongly depend on local factors [5,13–19]. According to [13], the NO<sup>2</sup> concentration is slightly higher at a lower relative humidity, whereas other authors found that the NO<sup>2</sup> concentration correlates positively with the relative humidity in all seasons, especially during winter [20]. The dependence between NO<sup>2</sup> and temperature was found to be weak, in general, except for two considerable positive correlations, which were obtained for July and December, and for which no particular explanation was found in [13]. Other studies found strong positive correlations between the NO<sup>2</sup> concentration and temperature for all seasons [20]. A negative correlation between wind speed and NO<sup>2</sup> concentration in all seasons except for summer was found by [19].

This paper presents a study of the diurnal, monthly, seasonal and interannual evolution of nitrogen dioxide concentration at an urban site in the southeast of Romania (Braila, 45.3◦ N, 27.9◦ E), which is influenced by road traffic mainly during the daytime. The role played by anthropogenic sources and the most relevant meteorological parameters that may affect the level of NO<sup>2</sup> pollution at the local scale is investigated. A comparison between in situ measurements of the NO<sup>2</sup> concentration near the ground and satellite measurements of NO<sup>2</sup> vertical columns is also performed. The timeframe for both in situ and satellite data is 2009–2014.

#### **2. Data and Methods**/**Experimental**

#### *2.1. Data*

The experimental data used in this study are from the automated monitoring air quality station called "BRAILA 1", denoted as "BR1", which is a traffic type station, one of the four available stations in Braila town, which are part of the National Network for Monitoring Air Quality (RNMCA, The Agency of Environmental Protection Braila: http://apmbr.anpm.ro/). The county of Braila, situated on the right shore of the Danube River, is considered to be the 11th Romanian city by the number of inhabitants. The number of vehicles in Braila has increased during the last few years, which led to a rather increased concentration of atmospheric pollutants. BR1 is a traffic station; it is located on Calea Gala¸ti No. 53 (The Agency of Environmental Protection Braila: http://apmbr.anpm.ro/) and monitors on a continuous basis the pollution levels generated mainly by traffic emissions, with medium and high flows, from the neighboring streets. Calea Galati Street is one of the busiest traffic streets in Braila. The width of the street is 16 m and it is covered with an asphalt carpet. The area around the air quality monitoring station has apartment blocks with four floors and some green spaces. The air quality monitoring station is located approximately 20 m from Calea Galati Street and 15 m from the apartment buildings. The surface is flat, characteristic of a plain area. The sources of pollution in the area consist mostly of road traffic and domestic heating.

The hourly atmospheric concentrations of NO<sup>2</sup> are determined in situ, by using the chemiluminescence technique [8,21]. This method implies the reduction of NO<sup>2</sup> to NO in the presence of a Molybdenum surface, at a temperature of roughly 310 ◦C. The resulted NO reacts with ozone (O3), leading to the formation of fluorescent NO2, whose emission is detected by a sensor [4,21]. Data between 2009 and 2014 were used in this study. Besides, concentrations of NO2, meteorological parameters were also recorded.

Satellite measurements were provided by the Ozone Monitoring Instrument (OMI) onboard the Earth Observing System Aura satellite [22,23]. The OMI measures several important pollutants, such as O3, NO2, SO2, and aerosols, with a spatial resolution of 13 km × 24 km, i.e., at a near urban scale resolution, https://aura.gsfc.nasa.gov/omi.html [24]. OMI is a remote sensing space instrument onboard the Aura satellite that provides a raw NO<sup>2</sup> product called DSCD (Differential Slant Column Density), which is retrieved using Differential Optical Absorption Spectroscopy (DOAS) in the 405–465 nm range. Satellite measurements provide valuable data about atmospheric pollutants, including NO<sup>2</sup> [25–32] at a large scale. The most refined product of OMI is the tropospheric Vertical Column Density (VCD), which is the result of a near real-time retrieval algorithm that gives a 0.7 × 10<sup>15</sup> molec./cm<sup>2</sup> uncertainty for each individual pixel [26]. In our study, we have used level 3 data, which are a conversion of OMI 13 × 24 km<sup>2</sup> data to a 0.25◦ × 0.25◦ resolution [33]. Such data collected between 2005 and 2015 by OMI confirmed a reduction in NO<sup>2</sup> concentration over Northern China, Eastern Europe, and USA, while an increase in NO<sup>2</sup> concentration was detected over the Persian Gulf and India [30].

#### *2.2. Methods*

Meteorological parameters recorded at an hourly rate were used, in order to verify to what extent they are relevant for modulating the NO<sup>x</sup> variability by means of correlation analysis. Correlations/anticorrelations (i.e., positive/negative correlation coefficients) between two variables are not the decisive proof for direct cause–effect relationships; however, they suggest that a link may exist if a physical mechanism supports it. When two variables are not correlated, the first thought is that there is no link between the two variables, which, however, may not always be true. Moreover, meteorological parameters (temperature, humidity, solar radiation) are not independent (e.g., solar radiation at the ground is a measure of the cloud cover, which is linked to relative humidity but also to temperature). Consequently, the correlation analysis was refined and a partial correlation analysis was used [34,35], complementary to the bivariate correlation. Such an analysis is useful when the relationship between two variables, e.g., *X* and *Z*, measured by the bivariate correlation coefficient, *CXZ*, may be induced or suppressed by a third (intervening) variable, e.g., *Y*, which affects both variables *X* and *Z*. The partial correlation *PX(Y)Z* corresponds to the link between *X* and *Z* when *Y* is constant. To be more specific, in this particular case the main variables, *X* and *Z*, were the NO<sup>2</sup> concentration and one meteorological factor (e.g., temperature) and the intervening variable, *Y*, was another meteorological parameter (e.g., radiation). The radiation (*Y*) may affect both the NO<sup>2</sup> concentration (*X*) and the temperature (*Z*). The difference *DX(Y)Z* = *PX(Y)Z* − *CXZ* between the partial correlation coefficient, *PX(Y)Z,* and the direct correlation coefficient, *CXZ,* measures the degree of intervention for *Y*. If *PX(Y)Z* is smaller than *CXZ*, i.e., the difference and the direct coefficient have opposite signs and *CXZ* ×*DX(Y)Z* < 0, then *Y* is responsible for part of the correlation, or it is said that Y "intervenes". If the two coefficients, *C* and *D*, have the same sign, i.e., *CXZ* × *DX(Y)Z* > 0, then the correlation is real and *Y* even suppresses, partially, the correlation. If the partial and direct coefficients are equal, i.e., *D* = 0, then the *Y* variable does not intervene. Table 1 summarizes possible results and interpretations for various combinations of direct and partial correlation coefficients.


**Table 1.** Type of correlation for various results of the partial correlation analysis.

The correlation analysis has been applied to standardized data (by their standard deviation) in order to identify possible links between NO<sup>2</sup> concentration in the atmosphere and the variations of several meteorological parameters.

#### **3. Results and Discussion**

#### *3.1. Diurnal and Seasonal Variation*

Figure 1 shows the diurnal and seasonal variation of NO<sup>2</sup> concentration (top) and temperature (bottom) during 2009–2014. The hourly averaged NO<sup>2</sup> concentration varied between 14 µg/m<sup>3</sup> (July 2009) and 41 µg/m<sup>3</sup> (January 2011). These values are smaller than those measured by traffic stations in Bucharest [28].

The NO<sup>2</sup> content varied with year; it was low in 2009 and in the first part of 2010, reached a maximum in the winter between 2010 and 2011, and then slightly decreased towards the end of the interval. The period is too short for assessing whether a trend does exist.

The seasonal variation of NO<sup>2</sup> was evident; NO<sup>2</sup> concentrations were clearly higher during winter than during summer. The seasonal variation may be explained by: (1) the strong variation of the anthropogenic sources contributing to NO<sup>2</sup> emissions, i.e., traffic and combustion, but also by (2) a longer lifetime of the NO<sup>2</sup> during winter, taking into account that the NO<sup>2</sup> lifetime and its concentration in the atmosphere are affected by the seasonal variation of the photochemical activity [7,8,26,31], which is reduced during winter. Additionally, during winter, the soil is cold, and thus the nearby air is heavy so that the emitted pollutants may stay close to the ground longer [2]. Another explanation for the seasonal NO<sup>x</sup> variation relates to the variation of the boundary layer height with the season, which influences the wind pattern that, in turn, has a very important role in pollutant dispersion [7].

#### *Int. J. Environ. Res. Public Health* **2020**, *17*, 6228

Two diurnal peaks of the NO<sup>2</sup> concentration could be observed, centered around 09:00 Local Time (LT) and 20:00 LT, which were more evident during cold months. The diurnal peaks are in agreement with previous findings [8] and are associated with peaks in traffic [7]. Morning peaks, observed around 9:00, were smaller than the evening peaks. The NO<sup>2</sup> concentration increased during the afternoon-evening period, between 17:00 and 24:00, with maxima progressing towards earlier time during winter months. This was seen in all seasons except summer. In December, these diurnal maxima were smaller, for all years, which can be linked to the reduced traffic and industrial activity due to holidays. *Int. J. Environ. Res. Public Health* **2020**, *17*, x 5 of 15 earlier time during winter months. This was seen in all seasons except summer. In December, these diurnal maxima were smaller, for all years, which can be linked to the reduced traffic and industrial activity due to holidays.

**Figure 1.** Diurnal and seasonal variation of NO2 concentration (**bottom**) and of temperature (**top**) during 2009–2014 (monthly means). **Figure 1.** Diurnal and seasonal variation of NO<sup>2</sup> concentration (**bottom**) and of temperature (**top**) during 2009–2014 (monthly means).

Figure 2 shows the variation of the tropospheric NO2 VCD (Vertical Column Density) over Braila, derived from OMI, together with in situ measurements at 11:00 LT, which is the approximate hour when the satellite overpasses the city. The tropospheric NO2 VCD varied between 0.8 × 10<sup>15</sup> molec./cm<sup>2</sup> for February 2011 and 4.7 × 1015 molec./cm<sup>2</sup> for December 2012. Note that the comparison in Figure 2 is strictly qualitative, since OMI measures the number of NO2 molecules in a column over a grid of 0.25° [33], while in situ instruments measure the volumetric concentration near the ground; thus a direct quantitative comparison makes no sense. The annual NO2 VCD ranges between 1.73 and 2.23 × 1015 molec./cm<sup>2</sup> . Previous studies have shown that OMI measurements give reliable results about NO2 emissions of anthropogenic sources at a large scale in cities [26,32,36], above large industrial platforms (located away from cities) [25,37] and given by vehicle emissions and also by residential emissions [38]. However, other studies have shown that satellite observations do not correctly evaluate the NO2 pollution caused by traffic or by other point sources [27,28]. Figure 2 shows the variation of the tropospheric NO<sup>2</sup> VCD (Vertical Column Density) over Braila, derived from OMI, together with in situ measurements at 11:00 LT, which is the approximate hour when the satellite overpasses the city. The tropospheric NO<sup>2</sup> VCD varied between 0.8 × 10<sup>15</sup> molec./cm<sup>2</sup> for February 2011 and 4.7 × 10<sup>15</sup> molec./cm<sup>2</sup> for December 2012. Note that the comparison in Figure 2 is strictly qualitative, since OMI measures the number of NO<sup>2</sup> molecules in a column over a grid of 0.25◦ [33], while in situ instruments measure the volumetric concentration near the ground; thus a direct quantitative comparison makes no sense. The annual NO<sup>2</sup> VCD ranges between 1.73 and 2.23 × 10<sup>15</sup> molec./cm<sup>2</sup> . Previous studies have shown that OMI measurements give reliable results about NO<sup>2</sup> emissions of anthropogenic sources at a large scale in cities [26,32,36], above large industrial platforms (located away from cities) [25,37] and given by vehicle emissions and also by residential emissions [38]. However, other studies have shown that satellite observations do not correctly evaluate the NO<sup>2</sup> pollution caused by traffic or by other point sources [27,28].

Figure 2 shows that the two time series agreed relatively well except for a period between the winter of 2009–2010 and the summer of 2011. The correlation coefficient between the two series was close to 0.40 (after values in December 2010 were removed) and did not change when NO<sup>2</sup> concentrations measured earlier, between 08:00 and 11:00 LT, were considered for the mean. Figure 2 shows that the two time series agreed relatively well except for a period between the winter of 2009–2010 and the summer of 2011. The correlation coefficient between the two series was close to 0.40 (after values in December 2010 were removed) and did not change when NO<sup>2</sup> concentrations measured earlier, between 08:00 and 11:00 LT, were considered for the mean.

The seasonal variation of the tropospheric VCD was more regular: it was small during warm months and high during cold months. In 2010, the two time series did not coincide, which is mainly due to a peculiar behavior of the in situ measurements that showed an unexpected wide maximum during summer and autumn, and a minimum in December. The OMI instrument "sees" a large surface area (312 km<sup>2</sup> ), and thus integrates the emissions of many sources. Consequently, it cannot accurately measure the local variability, captured by the in situ measurements. This opposing The seasonal variation of the tropospheric VCD was more regular: it was small during warm months and high during cold months. In 2010, the two time series did not coincide, which is mainly due to a peculiar behavior of the in situ measurements that showed an unexpected wide maximum during summer and autumn, and a minimum in December. The OMI instrument "sees" a large surface area (312 km<sup>2</sup> ), and thus integrates the emissions of many sources. Consequently, it cannot accurately measure the local variability, captured by the in situ measurements. This opposing variation

variation suggests that a punctual, temporary, source of NO2 was added in the second part of 2010

close to the measuring station.

suggests that a punctual, temporary, source of NO<sup>2</sup> was added in the second part of 2010 close to the measuring station. *Int. J. Environ. Res. Public Health* **2020**, *17*, x 6 of 15

**Figure 2.** Variation of standardized monthly averages of tropospheric NO2 (VCD measured by OMI, blue, dash) and of NO2 ground concentration measured in situ at 11:00 LT (black, solid) starting with January 2009 until December 2014. **Figure 2.** Variation of standardized monthly averages of tropospheric NO<sup>2</sup> (VCD measured by OMI, blue, dash) and of NO<sup>2</sup> ground concentration measured in situ at 11:00 LT (black, solid) starting with January 2009 until December 2014.

Figure 3 shows the average diurnal variation of the NO2 concentration during each season for the entire period. Seasons are defined here in a slightly different manner than usual: spring and fall/autumn consist of two months each (March–April and September–October), while solstice seasons have four months: the summer covers the period between May and August, and the winter starts in November and ends in February. This definition of seasons is justified by the weather profile during the last 20 years in Braila County and in South-Eastern Romania, which shows a fast transition from winter to summer and then back to winter. Figure 3 shows the average diurnal variation of the NO<sup>2</sup> concentration during each season for the entire period. Seasons are defined here in a slightly different manner than usual: spring and fall/autumn consist of two months each (March–April and September–October), while solstice seasons have four months: the summer covers the period between May and August, and the winter starts in November and ends in February. This definition of seasons is justified by the weather profile during the last 20 years in Braila County and in South-Eastern Romania, which shows a fast transition from winter to summer and then back to winter.

In general, the diurnal variation was relatively regular for fall, when a relatively small data spread was seen. A higher variability was observed in winter, spring and summer, when the data spread was larger. There were no outliers in the middle part of the day, when the traffic is slightly reduced compared to the morning and afternoon. The outliers during morning and afternoon may be associated with peaks or drop-offs of the traffic. The least regular season was summer, when the number of outliers was high regardless of the hour and their large majority lied in the upper part of the plot. This suggests that the NO2 loading was significantly higher than the average for a short period of time. This is confirmed by the next analysis (Figure 4). The explanation might relate to road construction works that resulted in deviation of the traffic. Starting with 2010, landslides occurred on streets near the measuring station. In 2011, the asphalt was removed and some excavations were done, in order to check the underground water and sewerage pipes. Subsequently, utility pipes were installed, the foundation was compacted and the traffic flow returned to normal values. In general, the diurnal variation was relatively regular for fall, when a relatively small data spread was seen. A higher variability was observed in winter, spring and summer, when the data spread was larger. There were no outliers in the middle part of the day, when the traffic is slightly reduced compared to the morning and afternoon. The outliers during morning and afternoon may be associated with peaks or drop-offs of the traffic. The least regular season was summer, when the number of outliers was high regardless of the hour and their large majority lied in the upper part of the plot. This suggests that the NO<sup>2</sup> loading was significantly higher than the average for a short period of time. This is confirmed by the next analysis (Figure 4). The explanation might relate to road construction works that resulted in deviation of the traffic. Starting with 2010, landslides occurred on streets near the measuring station. In 2011, the asphalt was removed and some excavations were done, in order to check the underground water and sewerage pipes. Subsequently, utility pipes were installed, the foundation was compacted and the traffic flow returned to normal values.

Apparently, the NO2 concentration and the temperature were anticorrelated, especially during daytime: high NO2 concentrations during morning corresponded to the lowest temperatures, while the afternoon reduction in NO2 coincided with the highest temperatures. However, this is just coincidental, since the daytime variation of NO2 is governed by traffic [8], whose peaks happen to occur at times when the temperature is lowest. Moreover, during nighttime, the NO2 concentration and the temperature seemed to be correlated. Apparently, the NO<sup>2</sup> concentration and the temperature were anticorrelated, especially during daytime: high NO<sup>2</sup> concentrations during morning corresponded to the lowest temperatures, while the afternoon reduction in NO<sup>2</sup> coincided with the highest temperatures. However, this is just coincidental, since the daytime variation of NO<sup>2</sup> is governed by traffic [8], whose peaks happen to occur at times when the temperature is lowest. Moreover, during nighttime, the NO<sup>2</sup> concentration and the temperature seemed to be correlated.

*3.2. Monthly Deviation of NO2 Concentration from the Seasonal Mean* 

Obviously, the diurnal variation of NO2 was not the same for each year. Figure 4 shows the difference between the hourly NO2 concentrations at ground level and the hourly seasonal mean during each month of each year. Hourly seasonal means are averages of hourly NO2 values measured during 2009 and 2014 over those months that define a season, i.e., November–February for

that particular time.

**Figure 3.** Diurnal variation of NO2 concentration (boxplots, left axis) for each season (top–down: spring, summer, fall and winter) for the entire interval, 2009–2014. See text for definition of seasons. The average temperature is superimposed (orange line, right axis). **Figure 3.** Diurnal variation of NO<sup>2</sup> concentration (boxplots, left axis) for each season (top–down: spring, summer, fall and winter) for the entire interval, 2009–2014. See text for definition of seasons. The average temperature is superimposed (orange line, right axis).

The analysis of the temporal variation of the aforementioned differences may be useful for identifying months or periods of the day when the NO2 loading departed from the expected variability, which in turn may help in finding the cause for the outliers seen in Figure 3. Values lying in the positive/negative part mean that the NO2 concentration were higher/lower than the average at

A clear pattern could be observed for equinox seasons: the afternoon peak was higher during colder months (March and October) than during warmer months (April and September) for each year. This was not true for the morning peak or for solstice seasons. The month of May, which is part of the summer season, was the least regular; the NO2 content in 2010 was much lower, while in 2011 it was much higher than the average. The unusual increase during the second part of 2010, shown in seasonal mean were positive for all years except 2011.

*Int. J. Environ. Res. Public Health* **2020**, *17*, x 8 of 15

Figure 2, is confirmed by the higher values seen in Figure 4 for July 2010 and, partially, for June 2010. The explanation might relate to the construction works previously described, which resulted in significant alterations of the traffic flow during 2010–2011. In general, the major contributor to the

**Figure 4.** Departures of hourly NO2 concentration from the hourly seasonal means for each month, January to December. Each year is shown with different colors/line styles (see legend). **Figure 4.** Departures of hourly NO<sup>2</sup> concentration from the hourly seasonal means for each month, January to December. Each year is shown with different colors/line styles (see legend).

#### *3.2. Monthly Deviation of NO<sup>2</sup> Concentration from the Seasonal Mean*

Obviously, the diurnal variation of NO<sup>2</sup> was not the same for each year. Figure 4 shows the difference between the hourly NO<sup>2</sup> concentrations at ground level and the hourly seasonal mean during each month of each year. Hourly seasonal means are averages of hourly NO<sup>2</sup> values measured during 2009 and 2014 over those months that define a season, i.e., November–February for winter, March–April for spring, May–August for summer and September–October for fall.

The analysis of the temporal variation of the aforementioned differences may be useful for identifying months or periods of the day when the NO<sup>2</sup> loading departed from the expected variability, which in turn may help in finding the cause for the outliers seen in Figure 3. Values lying in the positive/negative part mean that the NO<sup>2</sup> concentration were higher/lower than the average at that particular time.

A clear pattern could be observed for equinox seasons: the afternoon peak was higher during colder months (March and October) than during warmer months (April and September) for each year. This was not true for the morning peak or for solstice seasons. The month of May, which is part of the summer season, was the least regular; the NO<sup>2</sup> content in 2010 was much lower, while in 2011 it was much higher than the average. The unusual increase during the second part of 2010, shown in Figure 2, is confirmed by the higher values seen in Figure 4 for July 2010 and, partially, for June 2010. The explanation might relate to the construction works previously described, which resulted in significant alterations of the traffic flow during 2010–2011. In general, the major contributor to the summer mean for all years came from the NO<sup>2</sup> content in August, since the departures from the seasonal mean were positive for all years except 2011.

The NO<sup>2</sup> content in winter also depended on month and year. November seemed to be the month with the highest contribution to the winter seasonal average of NO<sup>2</sup> concentration for the first part of the interval (2009–2011). In December, the NO<sup>2</sup> diurnal variability changed with the year: the NO<sup>2</sup> content was lower during the first three years and higher afterwards. We assume that the explanation lies in a combination of meteorological conditions and important variations of traffic and industrial emissions during these particular years. During the analyzed period, the distribution of thermal energy and hot water in Braila municipality was mainly controlled by the heating operator SC "CET" SA. Starting with 2012, the lack of investments in modernizing the heating system affected the control of emissions: severe losses and the evolution of the methane gas tariff led to an increase of classical housing heating systems, whose impact is important especially during the cold season [39].

#### *3.3. Correlation between NO<sup>2</sup> Concentration and Meteorological Parameters*

The effect of meteorological parameters on the variability of NO<sup>2</sup> concentration for urban sites is, still, a conundrum. Table 2 shows some examples of correlations between NO<sup>2</sup> and meteorological factors for different locations [13–20]. The correlation between temperature and humidity, on the one side, and NOx, on the other side, was positive, negative or insignificant. The correlation with the wind speed was more consistent, i.e., is negative for all sites, which is rather normal, since a stronger wind will disperse the NO<sup>2</sup> at a local urban site. This is the reason for investigating whether meteorological parameters may be linked to the variability of the NO<sup>2</sup> in a relatively small city, with an average level of pollution [31] and whether this relationship depends on seasons or on time of day.


**Table 2.** Correlations between NO<sup>2</sup> and meteorological parameters for various locations.

Figures 2–4 show that during about 07:00 and 21:00 LT, the NO<sup>2</sup> variability was strongly influenced by traffic, which is also confirmed by [8,37,40–42]. The correlation coefficients were computed separately for the full 24 h time (black bars), for daytime (8–21 LT, red bars) and for the nighttime (22–7 LT, blue bars). The separation between daytime and nighttime was done because the influence of the road traffic should be lower during nighttime, and thus meteorological parameters may play a different role in modulating the NO<sup>2</sup> concentration during the night. Obviously, this is not valid for radiation, which is absent during nighttime. However, one should keep in mind that the NO<sup>2</sup> content is largely controlled by traffic, especially during the day, and this will definitely affect the correlation with meteorological parameters (or lack thereof). However, we assumed that the traffic pressure does not change for the analysis period.

Figure 5 shows the variation of the direct bivariate correlation between NO<sup>2</sup> and temperature (left), and humidity (right) with time. Correlation coefficients calculated for the full day are shown by black bars, while for day (night) these are shown by red (blue) bars. Only coefficients that were significant at 90% are shown. *Int. J. Environ. Res. Public Health* **2020**, *17*, x 10 of 15

**Figure 5.** Full bivariate correlation of NO2 with the temperature (**left**) and humidity (**right**). Coefficients (>90% confidence level) are shown for 24 h by black bars, while red/blue bars correspond to day/night, respectively—see text for details. **Figure 5.** Full bivariate correlation of NO<sup>2</sup> with the temperature (**left**) and humidity (**right**). Coefficients (>90% confidence level) are shown for 24 h by black bars, while red/blue bars correspond to day/night, respectively—see text for details.

There was no clear NO2 dependence on temperature, since correlation coefficients were positive for March, May, August and November, while for February and September these were negative. During daytime most correlations were negative; however, this was already discussed as being artificially induced by rush-hour traffic. Correlations did not change from day to night, except for May, when the correlation changed from negative for the full-day to positive for the full day and night. The full 24 h correlation was rather small or insignificant for most months and changed from positive to negative during months that had similar meteorological characteristics, e.g., March (positive), April (negative) and May (positive). Unsurprisingly, NO2 and temperature were anticorrelated during the day, but this was already discussed as being mostly artificial. Significant correlation during the night was positive in March, May, August and November, and negative in January, February and September. All in all, temperature seemed to play no clear role in the There was no clear NO<sup>2</sup> dependence on temperature, since correlation coefficients were positive for March, May, August and November, while for February and September these were negative. During daytime most correlations were negative; however, this was already discussed as being artificially induced by rush-hour traffic. Correlations did not change from day to night, except for May, when the correlation changed from negative for the full-day to positive for the full day and night. The full 24 h correlation was rather small or insignificant for most months and changed from positive to negative during months that had similar meteorological characteristics, e.g., March (positive), April (negative) and May (positive). Unsurprisingly, NO<sup>2</sup> and temperature were anticorrelated during the day, but this was already discussed as being mostly artificial. Significant correlation during the night was positive in March, May, August and November, and negative in January, February and September. All in all, temperature seemed to play no clear role in the variability of NO2.

variability of NO2. Negative coefficients were found for the NO2–temperature dependence by [14,16–18], without a clear association with seasons. The NO2 lifetime is higher during winter; thus the "night" analysis may be less relevant during winter because of the lower concentrations of the N2O and N2O5 species. These species are key factors in the removal of NO2 during nighttime and in its transformation into nitric acid [43]. The amount of NOy species is higher during warmer seasons when bacteria and agriculture activities are intensified [44]. However, [15] and [19] found that an increase in Negative coefficients were found for the NO2–temperature dependence by [14,16–18], without a clear association with seasons. The NO<sup>2</sup> lifetime is higher during winter; thus the "night" analysis may be less relevant during winter because of the lower concentrations of the N2O and N2O<sup>5</sup> species. These species are key factors in the removal of NO<sup>2</sup> during nighttime and in its transformation into nitric acid [43]. The amount of NO<sup>y</sup> species is higher during warmer seasons when bacteria and agriculture activities are intensified [44]. However, [15,19] found that an increase in temperature would be followed by an increase in NOx.

temperature would be followed by an increase in NOx. The correlation of the NO2 concentration with the humidity, when significant, was positive for The correlation of the NO<sup>2</sup> concentration with the humidity, when significant, was positive for most months. The only exception was May, during which the NO<sup>2</sup> variability departed significantly

In the following the intervening effect on direct correlations is analyzed, to see whether the existing correlations were spuriously induced, especially for the link to humidity. The partial correlation was used to identify whether the effects of meteorological parameters (temperature, relative humidity and solar radiation) on NO2 concentration were interconnected. The main and intervening variables (described in Section 2) are shown in Table 3. The results are shown only for four out of the six possible combinations, because these are the most relevant for identifying possible intervening effects. The correlation between NO2 and temperature was inconsistent and when assessing the intervening effect of radiation on correlation with humidity one can infer that the

Figure 6 shows the results of the partial correlation analysis for the combinations in Table 3. Full bars describe the direct correlation, while empty bars stand for the differences between partial and direct coefficients. According to Table 1, when these two are opposite, the correlation is partially mediated by the intervening variable. When these both have the same sign, the correlation is real.

September.

intervening effect of humidity was similar.

from the expected behavior (Figure 4). The correlation did not change from day to night except in September.

In the following the intervening effect on direct correlations is analyzed, to see whether the existing correlations were spuriously induced, especially for the link to humidity. The partial correlation was used to identify whether the effects of meteorological parameters (temperature, relative humidity and solar radiation) on NO<sup>2</sup> concentration were interconnected. The main and intervening variables (described in Section 2) are shown in Table 3. The results are shown only for four out of the six possible combinations, because these are the most relevant for identifying possible intervening effects. The correlation between NO<sup>2</sup> and temperature was inconsistent and when assessing the intervening effect of radiation on correlation with humidity one can infer that the intervening effect of humidity was similar.


**Table 3.** Meteorological parameters as main/intervening variables, for the partial correlation analysis.

**Figure 6.** Partial **c**orrelation between NO2 concentration and meteorological parameters. Full bars describe the bivariate correlation and empty bars describe the difference (*D*) between the partial and bivariate correlation. Correlation coefficients for the whole day are shown in black, while red/blue correspond to day/night, respectively—see text for details. Coefficients are shown for the >90% confidence level. Partial correlation between **NO2 and humidity** is shown in (**a**) (*temperature* constant) and (**b**) (*radiation* constant). Partial correlation between **NO2 and radiation** is shown in (**c**) (temperature constant) and the effect of *radiation* on the correlation with **temperature** is in (**d**). **Figure 6.** Partial correlation between NO<sup>2</sup> concentration and meteorological parameters. Full bars describe the bivariate correlation and empty bars describe the difference (*D*) between the partial and bivariate correlation. Correlation coefficients for the whole day are shown in black, while red/blue correspond to day/night, respectively—see text for details. Coefficients are shown for the >90% confidence level. Partial correlation between **NO2 and humidity** is shown in (**a**) (*temperature* constant) and (**b**) (*radiation* constant). Partial correlation between **NO2 and radiation** is shown in (**c**) (temperature constant) and the effect of *radiation* on the correlation with **temperature** is in (**d**).

NO2 relative humidity temperature a NO2 relative humidity radiation b NO2 radiation temperature c NO2 temperature radiation d

No consistent direct correlation between NO2 and humidity was found for the whole 24 h period (Figure 6a, full black bars are completely absent). Notably, *DNO(t)H* was rather high and positive for the whole 24 h day during May, June and December. This means that the temperature suppressed a potentially positive NOx–humidity correlation. During daytime, NO2 and humidity were positively correlated (red full bars) for a large part of the year, except for February, August and November. The differences between partial and direct correlation, *DNO(t)H*, albeit small, were of opposing signs for most cases. This suggests that the temperature was partially responsible for the

The important role played by the humidity in the variation of NO2 is supported by Figure 6b, where the effect of radiation on the same correlation (NO2–humidity) is shown only for daytime. The solar radiation partially induced a positive correlation during some months, but the effect was not important, since all *DNO(R)H* were pretty small. Our results agree with [19] or [20], but contradict [13], who showed that NO2 was negatively correlated with relative humidity. They argued that NO<sup>2</sup>

analysis.

NO2–humidity correlation.

**Table 3.** Meteorological parameters as main/intervening variables, for the partial correlation

Figure 6 shows the results of the partial correlation analysis for the combinations in Table 3. Full bars describe the direct correlation, while empty bars stand for the differences between partial and direct coefficients. According to Table 1, when these two are opposite, the correlation is partially mediated by the intervening variable. When these both have the same sign, the correlation is real.

No consistent direct correlation between NO<sup>2</sup> and humidity was found for the whole 24 h period (Figure 6a, full black bars are completely absent). Notably, *DNO(t)H* was rather high and positive for the whole 24 h day during May, June and December. This means that the temperature suppressed a potentially positive NOx–humidity correlation. During daytime, NO<sup>2</sup> and humidity were positively correlated (red full bars) for a large part of the year, except for February, August and November. The differences between partial and direct correlation, *DNO(t)H*, albeit small, were of opposing signs for most cases. This suggests that the temperature was partially responsible for the NO2–humidity correlation.

The important role played by the humidity in the variation of NO<sup>2</sup> is supported by Figure 6b, where the effect of radiation on the same correlation (NO2–humidity) is shown only for daytime. The solar radiation partially induced a positive correlation during some months, but the effect was not important, since all *DNO(R)H* were pretty small. Our results agree with [19] or [20], but contradict [13], who showed that NO<sup>2</sup> was negatively correlated with relative humidity. They argued that NO<sup>2</sup> concentrations are slightly higher at a lower relative humidity because the reactions between NO<sup>2</sup> and OH are less frequent, and thus the NO<sup>2</sup> persists more in the atmosphere. However, this was not confirmed by our results.

The NO<sup>2</sup> concentration was negatively correlated with solar radiation during the entire year. A higher direct radiation implies, usually, a higher air temperature; thus the intervening effect of temperature and radiation was analyzed. Figure 6c shows that the temperature did not artificially induce the anticorrelation with solar radiation (small or absent empty bars), except for May and September, when the effect was, however, small. This holds for the intervening effect of solar radiation on the anticorrelation with temperature, (Figure 6d), which was also small.

Based on observations in [45], which showed that, for a site in India, O<sup>3</sup> correlated negatively with both NO<sup>x</sup> and the humidity during all seasons, we suggest that the positive correlation between NO<sup>x</sup> and humidity may be an indirect result of the photochemical effect of solar radiation and humidity on ozone. Unfortunately, ozone measurements were not available to confirm this hypothesis. This also may partly explain the observed anticorrelation between NO<sup>2</sup> and solar radiation. Increased solar radiation favors the production of O3, which, in turn, reduces the NO<sup>2</sup> loading in the atmosphere [45].

In general, comparisons with other studies are not straightforward, since we are not aware of similar investigations of the intervening effect of meteorological parameters. Additionally, most studies did not consider monthly changes. Moreover, correlations between the NO<sup>2</sup> concentration and the meteorological parameters should be different for different cities, since the anthropogenic landscape and microclimate change significantly from one urban location to another [10,31].

#### **4. Conclusions**

This paper describes the diurnal, monthly, seasonal and annual evolution of the NO<sup>2</sup> concentration for 2009–2014, measured in situ by an urban traffic station in southeast Romania. The role that meteorological parameters might play in modulating the NO<sup>2</sup> variability was investigated in an attempt to separate the anthropogenic effects (which are well-known) from the effect of the local microclimate.

As expected, the NO<sup>2</sup> was higher during the cold season, except for one year, 2010, when summer NO<sup>2</sup> levels were highest. This suggests that natural factors, such as the effect of temperature on the NO<sup>2</sup> lifetime, are less important than the anthropogenic ones at urban sites. Some annual variation also existed, with low values at the beginning of the interval (2009–2010), most likely caused by a severe reduction of industrial activity.The summer minima and winter maxima have both anthropogenic and natural causes and the departure from these may relate to temporary changes of the local traffic and/or construction activities. The NO<sup>2</sup> diurnal variability was clearly shaped by the local transport: two diurnal peaks were observed, one around 8–9 LT and another one around 20–21 LT, and both were associated with increased road traffic, confirming previous observations at other urban sites. The afternoon peak was higher during the colder months (March and October) than during the warmer months (April and September) for each year. An irregular diurnal variation of the NO<sup>2</sup> concentration was seen in May and December. The most consistent season was autumn, with a relatively similar diurnal variation in all years.

Additionally, we found that over Braila, space observations of OMI followed the in situ observations during most of the selected interval (R = 0.4).

The analysis of the correlations between the NO<sup>2</sup> concentration and temperature, relative humidity and radiation has shown that the association with temperature is the least relevant. The correlation changed from positive to negative throughout the year without a clear pattern. Obviously, the contribution of traffic cannot be disregarded and may mask or suppress the impact of temperature variations on the NO<sup>2</sup> concentration. Our assumption that during the night, the situation may change due to traffic disappearance, was not confirmed. The correlations with the humidity and radiation, on the other hand, were notably consistent: the NO<sup>2</sup> concentration correlated with the relative humidity and was anticorrelated with radiation for almost the entire year. Moreover, most of these relationships were real and the intervening effect of the other meteorological parameters was small.

Our results showed that finding a link between meteorological parameters and NO<sup>2</sup> variability for an urban site is a difficult task. Attempts to predict the NO<sup>2</sup> behavior based on meteorological data, even combined with traffic flux data, cannot be very successful at the local or regional scale or on a short-term basis, since landscape, infrastructure, traffic, local activities, and population clearly affect the NO<sup>2</sup> concentration. However, this may be useful for assessing trends or long-tern variability at a large scale, since these are averaging over a large area, thus reducing local and short-term contributions.

**Author Contributions:** Conceptualization, M.V.; methodology, M.V. and D.-E.C.; formal analysis, M.V.; investigation, M.V. and V.C.; resources, D.-E.C., C.M.D.B., S.C.-B., V.C. and A.R.; data curation, C.M.D.B.; writing—original draft preparation, M.V.; writing—review and editing, M.V., D.-E.C., S.C.-B.; funding acquisition, M.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the EXPERT project, contract no. 14PFE/17.10.2018 with the Romanian Ministry of Research and Innovation.

**Acknowledgments:** We acknowledge the free use of in situ measurements of the NO<sup>2</sup> concentrations of the National Network for Monitoring Air Quality (RNMCA). Space-based VDC of NO<sup>2</sup> (OMI) were provided by the Giovanni online data system, developed and maintained by the NASA GES DIS (https://disc.gsfc.nasa.gov/ datasets/OMNO2d\_V003/summary).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **Key Points in Air Pollution Meteorology**

## **Isidro A. Pérez \* , Mª Ángeles García , Mª Luisa Sánchez, Nuria Pardo and Beatriz Fernández-Duque**

Department of Applied Physics, Faculty of Sciences, University of Valladolid, Paseo de Belén, 7, 47011 Valladolid, Spain; magperez@fa1.uva.es (M.Á.G.); mluisa.sanchez@uva.es (M.L.S.); npardo@fa1.uva.es (N.P.); beatriz.fernandez.duque@uva.es (B.F.-D.)

**\*** Correspondence: iaperez@fa1.uva.es

Received: 22 September 2020; Accepted: 9 November 2020; Published: 11 November 2020 -

**Abstract:** Although emissions have a direct impact on air pollution, meteorological processes may influence inmission concentration, with the only way to control air pollution being through the rates emitted. This paper presents the close relationship between air pollution and meteorology following the scales of atmospheric motion. In macroscale, this review focuses on the synoptic pattern, since certain weather types are related to pollution episodes, with the determination of these weather types being the key point of these studies. The contrasting contribution of cold fronts is also presented, whilst mathematical models are seen to increase the analysis possibilities of pollution transport. In mesoscale, land–sea and mountain–valley breezes may reinforce certain pollution episodes, and recirculation processes are sometimes favoured by orographic features. The urban heat island is also considered, since the formation of mesovortices determines the entry of pollutants into the city. At the microscale, the influence of the boundary layer height and its evolution are evaluated; in particular, the contribution of the low-level jet to pollutant transport and dispersion. Local meteorological variables have a major influence on calculations with the Gaussian plume model, whilst some eddies are features exclusive to urban environments. Finally, the impact of air pollution on meteorology is briefly commented on.

**Keywords:** particulate matter; atmospheric boundary layer; weather types; Gaussian plume model; low-level jet; recirculation; microscale; macroscale; mesoscale

#### **1. Introduction**

Air pollution is the subject of active research due to its marked impact on human life in a number of different environments, with regard to materials as well as living beings. Its impact on materials may have not only aesthetic but also economic implications. For example, the marginal costs of cleaning and repairing building façades due to air pollution have mainly been attributed to PM<sup>10</sup> and SO<sup>2</sup> [1], while in the plant environment, a key role is played by biomonitors. *Cupressus macrocarpa* leaves were recently used to evaluate the pollution load, and high values of anthropogenic elements, such as Cu and As, were observed near an industrial area [2]. Finally, the effects on human health have also been extensively investigated due to the global scope of the air pollution problem [3]. For instance, a positive and significant association between air pollutant concentration and the incidence of breast cancer has been described by Hwang et al. [4]. Moreover, the relationship between atmospheric levels of pollutants such as particulate matter, ozone, nitrogen dioxide, and sulphur dioxide and mortality causes has been reported in Portugal [5], with the increase in health expenditure due to air pollution being the economic implication of this issue [6].

Substances released into the atmosphere become pollutants when their concentrations exceed certain levels. These values may be reached not only by the increase in emissions but also under given meteorological conditions. The main difference between the roles of emissions and meteorological conditions is that emissions may be controlled when inmission levels are exceeded [7,8], whereas meteorological conditions cannot be controlled. Consequently, air pollution control strategies should take into account meteorological variables, since these will shape any such strategies to a certain degree. Kanawade et al. [9] analysed severe air pollution episodes in New Delhi, India, which were linked to stagnant atmospheric conditions. During the study period, low-level winds were below 8 ms−<sup>1</sup> , there was no precipitation, and 500 hPa wind speeds were below 13 ms−<sup>1</sup> . However, meteorology may also influence air pollution transport. Dong et al. [10] analysed the seasonal impact of meteorological variations on the air quality of the Beijing–Tianjin–Hebei region, China, which is a populated metropolitan region comparable to the Pearl River Delta or the Yangtze River Delta regions, where pollution transport plays a key role. One result of this study is the successful pollution control undertaken by joint regional controls. Moreover, the noticeable impact of meteorological conditions, up to 40% on regional transport over two seasons, i.e., spring and winter, suggested that regional air pollution needed to be controlled in both seasons.

The objective of this paper is to highlight the feedback between air pollution and meteorology. Consequently, the first goal of this paper is to simplify this topic, in which a range of disciplines, such as chemistry, physics, numerical modelling, biology, medicine, or even architecture, are involved. Moreover, this paper is not only aimed at pointing out the research lines that are currently active but also at highlighting possible future studies in order to attract potential researchers to a field that stands out thanks to its social impact and constant expansion.

Although the contribution of meteorological processes to air pollution events is usually assumed, it is frequently relegated to the background. Physical processes in the atmosphere are varied and include clear-sky surface solar radiation reduction by aerosols [11] or the study of their scattering angle [12], temperature inversions associated with air pollution and health effects [13], or analysis of electricity in the convective boundary layer [14]. Moreover, meteorological processes may have contrasting effects. Meteorological conditions sometimes determine air pollution situations, since they may aggravate pollution episodes, such as the recirculation of air parcel trajectories studied in Los Angeles [15]. However, certain meteorological conditions may also improve air quality, with the inverse relationship between particulate matter concentrations and the ventilation factor (the product of wind speed at 10 m and the mixed-layer height estimate) observed in Santiago City, Chile, being one noticeable example [16].

This study differs from previous analyses in the approach used, since it takes dynamic meteorology as a basis, due to the significant weight this subject carries in atmospheric analyses. Since energy in the atmosphere flows from large eddies to molecular processes, the structure of this paper follows the scale of atmospheric motions: macroscale, mesoscale, and microscale. Macroscale considers weather phenomena with a diameter above 1000 km and with a life of several days or even weeks, such as synoptic and tropical cyclones or fronts. Mesoscale corresponds to middle-sized atmospheric structures with a diameter between 2 and 1000 km and which last from hours to a couple of days, such as mountain waves or sea breezes. Finally, microscale is linked to processes of between less than a metre and a couple of kilometres and lasting from a few seconds to several minutes. Energy exchanges at the low atmosphere together with surface conditions determine atmospheric processes at this scale, such as tornadoes, dust devils, or plumes. This classification must be considered as the result of a simplification method to separate spatial and temporal ranges. However, in agreement with atmospheric evolution, processes may sometimes appear mixed, and their influence on air pollution needs to be evaluated.

Urban environments merit special consideration since large numbers of people live, work, meet, or spend time at these sites. Urban activities are sources of pollutants, and their emissions depend on city populations and facilities. Pollutant dispersion is conditioned by urban features, and the spread of cities may vary following cultural and geographical characteristics. Consequently, the impact of cities ranges from the micro to the mesoscale, and this review considers both possibilities.

A second key point of this study is its review of the influence of air pollution on meteorology, which has scarcely been analysed to date. This impact may be noticeable at a local scale, such as the effects observed by Saha et al. [17] during the eve of the Deepawali festival at Kolkata, India, when a massive firework display triggered a marked increase in particulate matter, the consequences of which were increased surface temperature and changes in the temperature vertical profile and diurnal pattern of relative humidity. However, the effects of air pollution may also be observed at a global scale. Wang et al. [18] simulated the long-range transport of anthropogenic aerosols from Asia, which moved on the north Pacific and led to the increase in cloud top height against preindustrial conditions, resulting in an invigoration of mid-latitude cyclones.

#### **2. Macroscale**

Two groups of papers investigate the influence of meteorology at this scale on the concentrations recorded. One considers the effect of the synoptic pattern, while the second focuses on pollutant transport, which is analysed using air parcel trajectory models.

#### *2.1. Synoptic Pattern*

Atmospheric flow may be classified by discrete circulation patterns in agreement with certain modes of atmospheric circulation. Both subjective and objective procedures may be used to achieve this goal. Dayan et al. [19] highlighted the importance of these classification methods in environmental phenomena such as air pollution, desert dust intrusions, and flash floods.

Dai et al. [20] used a classification method for regional particulate pollution days in Eastern China based on a subjective approach that employs only three patterns: equalised pressure (EQP), the advancing edge of a cold front (ACF), and the inverted trough of low pressure (INT). The seasonal features of these patterns were also considered to form eight types, three EQPs (in summer, autumn, and winter), two ACFs (in autumn and winter), and three INTs (in spring, autumn, and winter). Most of the pollution days appeared in winter, with the lowest number appearing in spring. Moreover, most of the pollution days were linked to the EQP, whereas the INT was the least frequent. The main weakness of subjective procedures is that they rely on the skill and experience of the classifier, since certain atmospheric flows may display features that are compatible with different patterns.

In contrast to this relatively simple and subjective classification, objective procedures are based on mathematical expressions that could be applied not only by users who have experience in atmospheric evolution but also, in some cases, by users who are beginning their career. One example is the Lamb–Jenkinson weather-type approach, which was considered by Li et al. [21] in their study into PM2.5 formation in North China. In this classification, the location of pressure centres was determined using gridded pressure data of a 16-point moveable region. Additionally, a small number of empirical rules was considered to classify the daily synoptic pattern as one of 26 weather types. As a result, three weather types were associated with high PM2.5 pollution levels. Another example of an objective classification is the obliquely rotated principal component analysis in T-mode, used by Mao et al. [22] to investigate air pollution over Wuhan in wintertime from 2014 to 2018. Three of five dominant synoptic patterns were linked to heavy pollution. A third example of objective weather synoptic classification is based on the self-organising map approach. Shu et al. [23] employed this procedure to analyse summertime ozone pollution over the Yangtze River Delta in China and found six predominant synoptic weather patterns: four ozone-polluted types and two ozone-clean types.

The synoptic pattern has a noticeable impact at the microscale level, with high-pressure systems being an example. They are linked to low wind speed, quite limited expansion of the boundary layer, and a reinforced inversion layer, which thereby increases pollutant concentration [24].

#### *2.2. Cold Fronts*

Yang et al. [25] analysed the weather systems associated with nine pollution episodes in Hong Kong. Four were linked with monsoons, two with tropical cyclones, one with a saddle-type pressure field, one with a warm high pressure, and one with a cold-front passage. Cold fronts are evolving structures that have a varied impact on air pollution. Jia et al. [26] analysed the "sawtooth cycles" in Beijing, although these features have also been observed in other mid-latitude areas, such as eastern North America, that are regularly crossed by the polar front. Sawtooths are linked with double fronts, where a weak cold front precedes a second and stronger cold front. Their period is irregular, frequently around five days, although some sawtooths extend up to a week and have been confused with regular, weekly cycles. The impact of this irregular synoptic cycle on particulate matter concentration is characterised by a slow increase over a few days followed by a sharp drop in just a few hours.

Rainfall associated with cold fronts determines the decrease of particulate concentrations [27]. Cold frontal passages usually favour atmospheric pollutant removal. Hu et al. [28] observed air pollution purification during the cold frontal passage associated with a strong cold-air outbreak, since surface wind speed and atmospheric boundary layer height increased. However, there is also evidence of quite different results. The formation and evolution of air pollution episodes is influenced by wind speeds and wind shears caused by nocturnal low-level jets and cold fronts [29]. Moreover, the chemical production of sulphate, nitrate, and ammonium may increase due to the uplift of air pollutant together with humidity, which facilitates the heterogeneous and aqueous-phase oxidation of precursors [30]. Another example is presented by Kang et al. [31], who observed that particulate matter may rise to the free troposphere from the North China Plain due to a cold frontal passage. These particles were transported by strong north-westerly frontal airflow and returned to the boundary layer at the Yangtze River Delta by synoptic subsidence under the high-pressure system that dominates once the frontal zone moves downstream. They concluded that cold fronts may be seen as potential carriers of atmospheric pollutants. Similarly, Zhou et al. [32] considered the influence of large-scale cold fronts and warm conveyor belts on the transport of continental aerosols located near the surface mixed with aerosols from the upper planetary boundary layer and middle troposphere. Finally, Song et al. [33] analysed a dust storm event in Northern China in early May 2017 and concluded that the front's activity was associated with increased wind speed and dust emission.

#### *2.3. Pollution Transport*

Air pollution transport has an impact not only on people's health but also on political relationships between countries who share common borders, with particulate matter transport from Poland to the Czech Republic being one notable example when the wind direction may be established, which occurs in around 50% of cases. This is true despite the prevailing flow being from the Czech Republic to Poland, which is typical of northeast Moravia due to the orographic influence of the Moravian Gate, and though emissions from both countries contribute equally with low wind speeds, since wind direction frequently changes under such conditions. Most of the pollution was attributed to days with unclear airflow or with a significant wind change [34]. Another example is the pollution from China that affects South Korea's air quality with a seasonal evolution [35]. In this analysis, agricultural burning and coal-fired heating were noticeable in summer and winter, respectively, while dust storms from the deserts impacted in spring.

Different models have been developed to explore air parcel trajectories. The HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) model has been widely used. It was developed by the NOAA (National Oceanic and Atmospheric Administration) and, since it was first described by Draxler and Hess [36], it has been regularly updated [37]. However, other models have also gained certain relevance. The FLEXPART (FLEXible PARTicle dispersion model) may simulate long-range transport, turbulent diffusion, deposition, or linear chemistry. Since it first appeared in the mid-1990s [38], the model has also improved its performance, physical–chemical parameterisations, input and output formats, or available software [39]. Another example is the METEX (METeorological data EXplorer) model [40], whose main feature is its easy application, since only the coordinates of the site of interest are required and air parcel trajectories may be systematically calculated.

Determining air parcel trajectories is essential in terms of knowing which locations are affected by hazardous substances that have been accidentally released, such as radionuclides [41]. Air parcel trajectory models may be a necessary tool to determine the atmospheric corridors that control air flow. However, a large number of air parcel trajectories need be considered if any firm conclusions are to be reached, since handling isolated trajectories may not prove to be easy. Procedures for grouping trajectories, such as clustering analysis, are required and possible applications vary. Salmabadi et al. [42] employed the total spatial variance method to identify four main dust corridors affecting Ahvaz, Iran. Three of these, which had westerly or north-westerly orientations, were associated with the Shamal winds and were particularly active in summer. Another corridor, with a southerly orientation, was linked to a prefrontal dust-storm mechanism and was active in spring and winter. Rana et al. [43] considered the "angle distance matrix procedure" to form six trajectory groups and obtained high particulate matter concentration when the air parcels followed the route from the Middle East through the Himalayan valley to Dhaka, Bangladesh, during the dry season. Moreover, theoretical trajectory calculations must be combined with measurements and satellite imagery in order to obtain a precise knowledge of potential sources and transport times [44]. Trajectory calculations are the basis of certain useful calculations, such as the potential source contribution function and concentration-weighted trajectories [45]. In both procedures, the study region is divided into cells and the number of trajectories calculated must be high enough to obtain reliable results. The potential source contribution function considers the relationship between trajectory endpoints above a given concentration level and the total endpoints in each cell, whereas concentration-weighted trajectories calculate the average weighted concentration in each cell.

#### **3. Mesoscale**

Typical mesoscale processes are breezes induced by land–sea interactions. Orographic features condition the transport and dispersion of air pollution, and this interaction is the subject of ongoing research [46]. Moreover, air parcel recirculation and urban heat island effects also lie within this scale.

#### *3.1. Breezes*

Land–sea breezes are formed at the interface of surfaces with contrasting features. At these sites, transport of particulate matter and its composition reveal its origin. Batchvarova et al. [47] measured aerosol concentration at Ahtopol (Bulgaria) and found a lower number of coarse particles for marine air masses. Another example is the Pearl River Delta, where interactions of synoptic winds and mesoscale breezes determined weak winds and air pollution accumulation, favouring intense photochemical reactions and the formation of an ozone "pool" [48]. Similarly, in Manila, Philippines, strong sea breezes disperse fine particulate matter inland during the day in spring, although breezes are the outflow path of these aerosols during the night. However, stagnant air determines high accumulated particulate concentration in autumn when the sea breeze is weaker [49]. Dasari et al. [50] indicated that differential heating between land and sea caused strong breezes that led to diurnal dispersion and distribution of pollutants in Neom, Saudi Arabia.

Several examples illustrate the role of orographic features in pollutant transport and dispersion. Bei et al. [51] analysed the influence of mountain–valley breeze circulations in the Guanzhong basin, China, and found that this breeze determines a convergence zone in the basin where pollutants return from the mountain and concentrations increase in the evening. Alternatively, different indicators may be used to characterise the relief effect on atmospheric flow. Lai and Lin [52] analysed wind direction, the Froude number, and the pitching angle of the surface level among other variables to investigate air pollution events in Taiwan. They considered ten categories of air pollution events following the airflow and the location of the air pollution event. Moreover, the Froude number was below 0.25 in most of these events, revealing the high stability of the airflow, which tends to split around the mountains. Another noticeable effect induced by the relief was described by Quan et al. [53], who observed that a

vertical vortex elevated ground pollutants at an altitude between 1.4 and 1.7 km where an elevated pollutant layer is formed and transported to Beijing by the southerly wind.

#### *3.2. Recirculation Processes*

Allwine and Whiteman [54] proposed integral quantities to describe stagnation, recirculation, and ventilation to characterise airflow at specific sites. Stagnation conditions may favour high pollutant concentrations. However, recirculation is quantified by the recirculation factor where 0 corresponds to straight-line transport without recirculation, and 1 is reached when the air parcel returns to its origin without net transport. In this latter case, concentrations may change when pollutants are affected by chemical reactions. One example is recirculation formed by sea breezes combined with upslope winds in the western Mediterranean basin during summer, and which may be considered as natural-photochemical reactors, since nitrogen oxides and precursors are transformed into oxidants, such as ozone, under strong insolation with residence times in the order of days [55]. Recirculation has also been observed in the Hong Kong area, another coastal site, where ozone pollution may be aggravated by this process [56]. Dimitriou and Kassomenos [57] analysed the impact of Saharan dust intrusions on continental Greece and concluded that certain infrequent events observed in August–September were due to dust recirculation by the Etesian winds linked to the extension of the Azores anticyclone in continental Europe. Since recirculation is a kind of eddy, orographic features, such as valleys, favour these processes. Quimbayo-Duarte et al. [58] analysed the air flow in a valley with two sections and observed recirculations when the air from the slopes converges at the valley floor and flows to the narrow section in the up-valley direction.

Another example of recirculation is the Fresno eddy [59], which is a cyclonic vortex formed during the night in the San Joaquin Valley of California. During the daytime, winds mainly enter through the San Francisco Bay area and transport precursors to the valley where ozone is formed. The air rises into the mountains at the southern end and leaves the valley. However, during the night-time, the air cannot leave the valley due to the cool wind drainage from the mountains and so returns north in a circular flow pattern, known as the Fresno eddy. Pollutants are forced to return to urban areas where more precursors are added [60]. Observations indicate that this eddy dominates ozone production in the central part of this valley and triggers many ozone exceedances if it is strong on apparently ventilated days [61].

#### *3.3. Urban Heat Island*

The features as well as the extension of certain urban structures can modify both thermal and mechanical turbulences. Absorption and emission of radiation by urban surfaces in urban environments differs from the same processes when they occur in rural surfaces. Moreover, the roughness determined by buildings is also greater. The urban heat island may enhance the existing vorticity, which produces direct thermal circulation, causing surface mesovortices. Mainhart et al. [62] observed that ozone formed outside the Saint Louis metropolitan area returned to the city due to such mesovortices.

The beneficial effect of urban heat island mitigation by urban cooling strategies on air quality is assumed. Green roofs and living walls, for example, are useful for removing air pollutants [63]. However, competing feedbacks should be considered since green and cool roofs reduce wind speed and vertical mixing during the day, which might lead to air stagnation and air pollution situations [64]. Similarly, Henao et al. [65] indicated that air pollution transport mechanisms in urban valleys could be altered by urban heat island mitigation and that air quality could worsen. In fact, Li et al. [66] presented the "urban aerosol pollution island", which is negatively correlated with the urban heat island in summer. Moreover, mitigation of the urban heat island by an increasing albedo weakened turbulent mixing, made the planetary boundary layer height decrease, and led to a stronger urban aerosol pollution island. Consequently, a balance needs to be found between urban heat island mitigation and its impact on air quality.

Droste et al. [67] introduced the "urban wind island effect" since, under certain atmospheric conditions, wind speed in the city may be higher than that corresponding to rural environments. This effect would be observed in the afternoon and would be due to the contrast between the city and the countryside in the atmospheric boundary layer growth, surface roughness, and ageostrophic wind. However, this effect has only been suggested and needs to be verified prior to gauging its influence on air pollution.

#### **4. Microscale**

At the shortest spatial scale, pollutant concentration is affected by the atmospheric boundary layer height. However, determining its depth may not prove to be easy, and the issue has been avoided in various studies where concentrations are related with meteorological variables measured at a local scale.

#### *4.1. Atmospheric Boundary Layer*

Su et al. [68] found strong interactions between the planetary boundary layer height and particulate matter concentration in winter, when the planetary boundary layer height is shallow and particulate matter concentration is high. Correlations between the planetary boundary layer height and particulate concentrations tend to be negative in general. However, the magnitude, significance, or sign of these correlations depend on the location, season, and meteorological conditions. Results indicate that this relationship is weak over clean regions but nonlinear over polluted regions.

The development of a thermal internal boundary layer at coastal sites is linked with pollution episodes. Wei et al. [69] obtained the boundary layer height with a ceilometer at Qinhuangdao, North China. Determining the thermal internal boundary layer was based on wind direction, when the wind blew from the sea, and boundary layer height. The boundary layer was thinner with a thermal internal boundary layer. If this layer was observed, particulate concentrations increased when the boundary layer height rose, since the dilution effect of the breeze is weak, and the boundary layer height was small. However, if the thermal internal boundary layer was not formed, particulate concentrations fell when the boundary layer height increased.

The atmospheric boundary layer, which develops during the day, evolves to a stable nocturnal boundary layer due to the radiative cooling of the ground and is topped by a residual layer, which is a reservoir of pollutants from daytime emissions and photochemical production. Observations indicate a wind speed maximum, known as low-level jet, near the top of the nocturnal boundary layer. This nocturnal jet may favour mixing between the two layers despite stable stratification.

The influence of low-level jets on pollutant concentrations varies. These jets may transport pollutants trapped in the residual layer and trigger episodes when mixed at the surface by convective turbulence after sunrise. One example of this process is presented by Sullivan et al. [70], who reported an ozone episode in the Baltimore–Washington DC urban corridor in June 2015, with concentrations between 70 and 100 ppb at around 1 km above the surface in the nocturnal residual layer. These concentrations were transported by the nocturnal low-level jet, and both the simulations and the measurements indicated high "next-day" ozone concentrations at sites along the southern New England region. A similar case is presented by Li et al. [71], where pollutants from the North China Plain were transported by the nocturnal low-level jet to Shenyang. Caputi et al. [72] suggested that the ozone of the residual layer is actively mixed in the nocturnal boundary layer when the low-level jet is strong. Under these conditions, ozone dry deposition causes low concentrations the following day. However, ozone is retained by the residual layer when the low-level jet is weak since the residual layer is more decoupled. In this situation, fumigation processes the following morning would determine high ozone concentrations near the ground. Wei et al. [73] also suggested an abrupt reduction in particulate matter concentration and an improvement in air quality due to enhanced vertical mixing linked to the low-level jets formed at the end of the pollution periods.

Wind speed and vertical turbulence control air pollution during the day. However, vertical turbulence is weak at night and vertical mixing is controlled by large-scale turbulence eddies whose increased number, due to low wind speed, might lead to a fall in concentration near the ground [74].

#### *4.2. Local Meteorological Variables*

Specialised devices are not usually available to most researchers. In these situations, some analyses relate pollutant concentrations with variables provided by meteorological stations.

Chen et al. [75] observed that temperature and air pressure had a marked influence on ozone concentrations in Beijing, although meteorological influence was weak in summer, since emissions of volatile organic compounds and nitrogen oxides increased in this season and led to ozone pollution episodes. The close relationship between particulate matter and meteorological variables has been the subject of several analyses. For instance, Meng et al. [76] concluded that the planetary boundary layer height and temperature difference of the inversion layer and its depth are variables that impact particulate matter concentrations. This analysis presented the ranges of these meteorological variables, where the correlation with air pollution was observed. In this line, certain thresholds of the meteorological variables should be established, since the relationship with air pollution may be the opposite. Wang et al. [77] considered 27 cities in China and studied surface solar irradiance and wind speed. Both variables displayed a similar trend between 1961 and 2011. The authors found that winds below 2.5 ms−<sup>1</sup> disperse air pollutants and enhance solar surface irradiance. However, winds above 3.5 ms−<sup>1</sup> enhance aerosol concentration, perhaps from dust storms, and surface solar irradiance was seen to be attenuated. Between the two thresholds, no relationship between wind speed and surface solar radiation could be established.

Local meteorology may have a marked influence when a hazardous substance is accidentally released or when an industrial facility is close to a population centre, and various mathematical models have been developed to investigate air pollutant dispersion. The Gaussian model is a classical formulation of plumes that presents a simple structure and allows multiple simulation possibilities. Moreover, it may easily be extended to more complex configurations of sources. ¸Sahin and Ali [78] determined the emergency planning zones around two nuclear power plants using this model, and Rohmah and Nurokhim [79] simulated I<sup>131</sup> dispersion around a nuclear facility. Similarly, Maués et al. [80] applied this model to evaluate the air quality of Volta Redonda, Brazil, where a large steel plant is located, and reported concentration curves around the plant showing the plume impact. Dispersion coefficients play a key role in this model. Although their calculation relies on empirical schemes, the parameterisations used may give unequal results and should be compared. Mao et al. [81] considered four dispersion schemes applied to the Prairie Grass experiments in 1956. They evaluated the performance of these schemes following the stability conditions and selected the Chinese National Standard from the results obtained for the different stability classes.

#### *4.3. Urban Micrometeorology*

Eddy formation is favoured by the rough elements and varied covers present in cities. Han et al. [82] compared the air quality between urban and rural areas in 35 major Chinese metropolitan regions. Although urban pollution is usually more marked than rural pollution in summer, both pollutions are noticeable in winter, with rural pollution being even stronger than urban pollution in certain cities. The authors found that the key factors for the longer spatial inhomogeneity in summer versus winter were emissions, urban structures, meteorology, topography, and humidity. At a smaller, daily timescale, Choi et al. [83] observed in Los Angeles that limited dispersion capacity and the relatively stable surface layer in the morning led to higher ultrafine particle concentrations than in the afternoon. In this period, wind speeds were higher and turbulence was stronger, with vertical turbulence intensity being the most effective factor controlling ultrafine particle concentrations. The effect of tall buildings on pollutant dispersion was studied by Aristodemou et al. [84], who found a concentration increase at high floors, whereas dispersion was favoured by tall buildings at low floors.

#### **5. Influence of Air Pollution on Meteorology**

Although the spread and intensity of atmospheric processes make them insensitive to external influences, the strength of air pollution might condition an area's weather or even its climate.

Wang et al. [85] studied the effect of aerosols on radiation. Surface temperature is seen to be reduced by solar radiation reduction and the temperature increase of the upper layer by solar radiation absorption or backscattering. As a result, atmospheric stratification becomes more stable. In fact, simulations with increasing aerosol optical depth have revealed a transition from unstable to stable stratification at rural sites, to weakly stable at suburban sites and to less unstable or neutral stratification at urban sites. Consequently, energy transfer from the surface is reduced, as is the boundary layer height. This impact of aerosol pollution on meteorology has been reported in numerous analyses. Nguyen et al. [86] also observed a decrease in shortwave radiation, temperature, planetary boundary layer height, and wind speed. Visibility degradation in Eastern Thailand has been studied by Aman et al. [87], who observed a marked seasonality as well as a relationship with wind direction. Moreover, this visibility degradation was linked to air pollution and meteorological conditions, such as small wet scavenging, low mixing height, and a high recirculation factor.

Evidence concerning the relationship between aerosols and precipitations is conflicting. Rosenfeld and Woodley [88] indicated that the number of cloud condensation nuclei in polluted air is over ten times higher than this number for the same volume of clean air. Consequently, under air pollution conditions, small droplets would be formed, which could not easily grow to reach precipitation size. Jirak and Cotton [89] observed this effect at elevated sites west of Denver and Colorado Springs, Colorado. However, Rosenfeld et al. [90] have explained both the decrease and the increase in rainfall by air pollution from a combination of radiative effects on surface heating and microphysical effects on invigorating the convection. Recently, Yang et al. [91] concluded that tropical cyclone precipitation was invigorated by aerosols over the Chinese mainland from 1980 to 2014. Moreover, light rain presented a decreasing trend against time, whereas the trend for heavy rains was seen to be increasing. In contrast, post-monsoon precipitation reduction in South Asia between 2004 and 2014 was explained by two factors. The first was aerosol loading, which triggered surface cooling, and the second was the temperature decrease by evaporation around irrigation hotspots [92].

#### **6. Conclusions**

Meteorological conditions have a direct influence on emission control in air pollution events. Although atmospheric processes are varied and frequently mixed, one or a few sometimes prevail and a direct relationship with high pollutant concentrations may be established.

At the synoptic scale, objective procedures provide for a systematic classification of weather types that can be linked to pollution events. However, pressure centres, which are associated with certain air mass features, move slowly and their fronts sweep the low atmosphere, with contrasting effects on air pollution. Air parcel trajectory models use meteorological observations and are a useful tool to investigate the sources, paths, and destinations of air pollution.

Air flow may be disturbed by coasts, mountains, or big cities. These features may determine horizontal or vertical vortices where pollutants may be trapped. Recirculation processes have been the subject of some studies, and their analysis continues to offer a promising line of research. Moreover, the contrast between large cities and their surrounding areas is evolving rapidly, and the corresponding changes in meteorological aspects or consequences in pollutant transport and dispersion are fields that remain open to study.

At the lowest scale, the boundary layer height is inversely correlated with pollutant concentrations. Boundary layer evolution is closely related with air pollution since a polluted residual layer may be formed at night and the low-level jet favours both pollutant transport and dispersion. The relationship between concentrations and the usual meteorological variables is common in air pollution studies. Cities are particular sites that favour turbulence and where the low atmosphere is more unstable than at rural sites. The Gaussian plume model enables pollutant concentration to be calculated and

is especially useful when dealing with the release of hazardous substances. Moreover, this model may be easily adapted to greater distances and varied conditions such as trapped plumes, fumigation, or orographic features. At this scale, thresholds of meteorological variables should be investigated following air pollution levels.

Finally, the impact of air pollution on meteorology should not be excluded, with radiation absorption by aerosols being one example. The relationship between aerosols and precipitation depends on the balance between microphysical and radiative effects, and the conditions linked to this balance should be analysed.

Consequently, the examples presented in this study demonstrate that the relationship between meteorology and air pollution merits in-depth enquiry, and air pollution control strategies can no doubt benefit from this research. Although these studies may appear to be local analyses, ascertaining common features remains open to further research in order to gain insights into the feedback between meteorology and air pollution.

**Author Contributions:** Conceptualization, I.A.P. and M.Á.G.; methodology, I.A.P.; writing—original draft preparation, I.A.P.; writing—review and editing, N.P. and B.F.-D.; funding acquisition, M.L.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This review was funded by the Regional Government of Castile and Leon, project number VA027G19.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*International Journal of Environmental Research and Public Health* Editorial Office E-mail: ijerph@mdpi.com www.mdpi.com/journal/ijerph

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-2552-5