Estimation of Particulate Matter Contributions from Desert Outbreaks in Mediterranean Countries (2015–2018) Using the Time Series Clustering Method

atmosphere Article Estimation of Particulate Matter Contributions from Desert Outbreaks in Mediterranean Countries (2015–2018) Using the Time Series Clustering Method Álvaro Gómez-Losada 1,* and JoséC. M. Pires 2   Citation: Gómez-Losada, Á.; Pires, J.C.M. Estimation of Particulate Matter Contributions from Desert Outbreaks in Mediterranean Countries (2015–2018) Using the Time Series Clustering Method. Atmosphere 2021,12, 5. https://dx.doi.org /10.3390/atmos12010005 Received: 17 November 2020 Accepted: 21 December 2020 Published: 23 December 2020 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2020 by the authors. LicenseeMDPI, Basel, Switzerland. This articleisanopenaccessarticledistributed under the terms and conditions of the CreativeCommonsAttribution(CCBY) license(https://creativecommons.org/ licenses/by/4.0/). 1Departamento de Estadística e Investigación Operativa, Facultad de Matemáticas, Universidad de Sevilla, 41012 Sevilla, Spain 2LEPABE—Laboratory for Process Engineering, Environment, Biotechnology and Energy, Faculty of Engineering, University of Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal; [email protected] *Correspondence: [email protected] or [email protected] Abstract: North African dust intrusions can contribute to exceedances of the European PM 10 and PM 2.5 limit values and World Health Organisation standards, diminishing air quality, and increased mortality and morbidity at higher concentrations. In this study, the contribution of North African dust in Mediterranean countries was estimated using the time series clustering method. This method combines the non-parametric approach of Hidden Markov Models for studying time series, and the definition of different air pollution profiles (regimes of concentration). Using this approach, PM 10 and PM 2.5 time series obtained at background monitoring stations from seven countries were analysed from 2015 to 2018. The average characteristic contributions to PM 10 were estimated as 11.6 ±10.3 µg·m−3 (Bosnia and Herzegovina), 8.8 ± 7.5 µ g · m −3 (Spain), 7.0 ±6.2 µg·m−3 (France), 8.1 ± 5.9 µ g · m −3 (Croatia), 7.5 ± 5.5 µ g · m −3 (Italy), 8.1 ± 7.0 µ g · m −3 (Portugal), and 17.0 ±9.8 µg·m−3 (Turkey). For PM 2.5 , estimated contributions were 4.1 ± 3.5 µ g · m −3 (Spain), 6.0 ±4.8 µg·m−3 (France), 9.1 ±6.4 µg·m−3 (Croatia), 5.2 ± 3.8 µ g · m −3 (Italy), 6.0 ± 4.4 µ g · m −3 (Portugal), and 9.0 ± 5.6 µ g · m −3 (Turkey). The observed PM 2.5 /PM 10 ratios were between 0.36 and 0.69, and their seasonal variation was characterised, presenting higher values in colder months. Principal component analysis enabled the association of background sites based on their estimated PM10 and PM2.5 pollution profiles. Keywords: African dust; air pollution; hidden Markov models; particulate matter; PM 2.5 /PM 10 ratio; principal component analysis 1. Introduction Air pollution is a mixture of particles and gases, which can reach unsafe concentrations for human health, the environment, and vegetation. It has become one of the main sustainability issues and a concerning topic in atmospheric science. According to the World Health Organisation (WHO), 90% of the world population lives in highly polluted environments, and about 7 million deaths are caused every year by outdoor and indoor air pollution [1]. Particulate matter (PM) is one of the most concerning air pollutants in Europe [ 2 , 3 ]. It comprises suspended particles in the air with variable composition and size, which results from several anthropogenic activities. The main sources are industrial facilities, power plants, traffic emissions, wood stove combustion, dust, and forest fires. PM can be divided into two main components [ 4 , 5 ]: (i) primary PM emitted directly by pollution sources; and (ii) secondary PM that is produced in the atmosphere by chemical reactions with primary gaseous pollutants (e.g., organic compounds, sulphur dioxide, and nitrogen oxides). PM is usually classified according to its size, with fine particles being designated as PM 2.5 (PM with an aerodynamic diameter <2.5 µ m). The group of fine (with an aerodynamic diameter Atmosphere 2021,12, 5. https://dx.doi.org/10.3390/atmos12010005 https://www.mdpi.com/journal/atmosphere Atmosphere 2021,12, 5 2 of 18 less than 2.5 µ m) and coarse (with an aerodynamic diameter between 2.5 and 10 µ m) particles is designated as PM 10 . The monitoring of PM concentrations according to their sizes is relevant due to the corresponding impact, in particular on human health [ 6 ]. Concerning the respiratory system, the upper respiratory tract is affected by PM 10 , while fine particles have the capacity to deeply penetrate in the human respiratory system until lung alveoli (therefore, more hazardous to human health than the coarser particles) [ 7 ]. Besides human health, PM is associated with environmental concerns (e.g., visibility impairment—haze) and the damage of ecosystems [ 8 , 9 ]. The impact depends on size, surface, number, and composition of particles, of which characteristics are mainly affected by anthropogenic and natural sources, meteorology (that influences their dispersion, transport, and deposition), and urbanisation level. The monitoring and the understanding of the temporal and spatial behaviours of PM concentrations are essential for the implementation of air quality policies and the definition of effective measures to mitigate PM presence and its effects. Besides the study of PM 2.5 and PM 10 , the ratio between their concentrations (PM 2.5 /PM 10 ) may give relevant information about PM sources and effects [ 10 – 13 ]. High values indicate an increased contribution of anthropogenic sources, while low values are associated with natural sources due to a high fraction of coarse particles. In Europe, the Directive 2008/50/EC defines the limits for the concentrations of several air pollutants, including PM (daily limit of 50 µ g · m −3 not to be exceeded more than 35 times per year for PM 10 , and the annual limit value of 20 µ g · m −3 for PM 2.5 ). However, the limit values for PM concentrations in Mediterranean countries are not easily fulfilled. The main reason is the strong influence of dust transport from the Sahara Desert (North Africa). The contribution of this natural source should be known for each country in order for the European legislation to be fairly applied. The directive allows the subtraction of daily limit value exceedances of PM 10 concentrations (50 µ g/m 3 ) when these episodes are attributed to natural events (articles 2.15 and 20.1). There are several statistical techniques that were applied to quantify the contribution of this natural source. Gómez-Losada et al. [ 14 ] applied the Hidden Markov Model (HMM) to estimate PM 10 contributions using the data of 33 background sites in Spain and Portugal. The model results were then compared with the method widely accepted to estimate the daily African PM 10 load: the 40th percentile (P40) method introduced by Escudero et al. [ 15 ]. The HMM model presented a better estimation of the Sahara dust contribution in the monitored values of PM 10 concentrations. Stafoggia et al. [ 16 ] distinguished PM 10 from desert and non-desert sources in 13 Southern European cities for 2001–2010. They combined modelling tools, back-trajectories, and satellite data to identify desert dust advection days. In 15% of the days, PM 10 concentrations were influenced by desert dust, with most episodes occurred in spring–summer. Gama et al. [ 17 ] estimated the contribution of desert dust to the PM 10 and PM 2.5 concentrations in background sites in Portugal during 2016. Two different approaches were applied: the 40th percentile (P40) method and the WRF-CHIMERE model (multi-scale chemistry-transport model). The P40 method identified higher African dust contributions than the modelling approach. The PM 10 contributions in most of the regions were between 0 and 90 µ g · m −3 , while for PM 2.5 this natural source was responsible for 0 to 30 µg·m−3. This study aims to estimate the particulate matter contributions from desert outbreaks. Specific objectives are: (i) to characterise PM pollution behaviour in different Mediterranean counties, (ii) to estimate the different PM 10 and PM 2.5 pollution profiles from Southern Europe remote monitoring stations, and once they are obtained, (iii) estimate these air pollutant contributions from desert dust outbreaks. 2. Data and Methods 2.1. Data PM 10 and PM 2.5 validated hourly time series (TS) from background stations of southern regions of Europe were obtained from the European Environment Agency data repository (DiscoMap, https://discomap.eea.europa.eu) for the 2015 to 2018 period. Only TS Atmosphere 2021,12, 5 3 of 18 with time coverage higher than 80% (7008 hourly observations per annual TS) were considered in this study and independently analysed for each year. The number of analysed TS from background sites fulfilling this data quality criterion is indicated in Table 1, together with their distribution by countries and years. A relation of the studied PM 10 and PM 2.5 background monitoring sites is provided as Supplementary Material (SM) in two comma-separated value files. In these files, the typology, main pollution source, and the geographical coordinates of monitoring stations are given in EPSG 4979 projection. The length of the analysed TS from each station by year is also provided. Table 1. Country distribution of PM 10 and PM 2.5 time series from background stations analysed in this study. Country Abbreviation 2015 2016 2017 2018 PM10 PM2.5 PM10 PM2.5 PM10 PM2.5 PM10 PM2.5 Bosnia and Herzegovina BA - - 1 - 1 - - - Spain ES 15 8 20 5 22 1 25 9 France FR 22 16 21 15 22 12 25 19 Croatia HR 3 3 3 2 4 3 3 3 Italy IT 7 2 8 3 7 - 11 4 Portugal PT 9 6 10 4 10 4 10 3 Turkey TR 5 - 7 4 7 3 12 6 For the years 2015, 2016, 2017, and 2018, a total of 61, 70, 73, and 86 PM 10 TS were studied, respectively; with respect to PM 2.5 , 35, 33, 23, and 44 TS, respectively, were studied for those years. All of these PM 10 and PM 2.5 TS were obtained from background monitoring stations. Tables 2and 3show summary statistics for dust desert days based on the exceedance of PM 10 and PM 2.5 daily average concentrations of 150 and 43 µ g · m −3 , respectively. These values were selected as threshold values due to the maximum levels of these pollutant concentrations in Southern Europe during dust outbreaks [18,19]. They correspond to a lower bound of ranges typically found in this region. Table 2. Average concentration by country for PM 10 time series during desert and non-desert days in background monitoring sites (in µ g · m −3 ). n, number of days with exceedance of 150 µ g · m −3 PM 10 daily average. BA- Bosnia and Herzegovina; ES—Spain; FR—France; HR—Croatia; IT—Italy; PT—Portugal; TR—Turkey. Country Year Dust Days Non-Dust Days nMin. p50 Max. Min. p50 Max. ES 2015 2 255.0 291.2 327.4 2.8 13.1 64.6 FR 0 - - - 2.5 12.9 68.4 HR 0 - - - 2.5 15.6 74.3 IT 0 - - - 2.4 17.7 75.7 PT 0 - - - 1.9 13.6 71.9 TR 21 167.0 261.6 444.9 7.0 31.3 125.3 BA 2016 3 158.3 193.0 196.8 2.1 16.1 107.1 ES 10 229.0 241.5 257.1 3.8 14.2 93.0 FR 0 - - - 2.5 10.9 56.2 HR 0 - - - 2.3 13.0 66.4 IT 0 - - - 2.8 14.7 98.2 PT 2 200.9 200.9 200.9 2.7 12.2 82.0 TR 19 263.7 267.2 338.4 4.5 21.4 105.9 Atmosphere 2021,12, 5 4 of 18 Table 2. Cont. Country Year Dust Days Non-Dust Days nMin. p50 Max. Min. p50 Max. BA 2017 0 - - - 2.3 13.9 68.9 ES 10 212.9 223.5 226.3 2.8 12.9 78,3 FR 0 - - - 2.3 10.8 55.6 HR 0 - - - 1.0 11.4 70.5 IT 0 - - - 1.9 15.2 72.1 PT 2 451.7 451.7 451.7 3.0 15.4 85.2 TR 24 152.6 232.0 334.3 7.5 27.4 113.1 ES 2018 5 179.7 192.8 202.3 2.7 11.9 63.4 FR 0 - - - 2.1 11.4 49.3 HR 0 - - - 3.1 16.0 77.0 IT 0 - - - 2.9 20.6 74.9 PT 0 - - - 1.9 11.4 70.5 TR 14 180.9 192.3 208.7 5.4 21.1 110.0 Table 3. Average concentration by country for PM 2.5 time series as in Table 1.n, number of days with exceedance of 43 µ g · m −3 PM 2.5 daily average. ES—Spain; FR—France; HR—Croatia; IT—Italy; PT—Portugal. Country Year Dust Days Non-Dust Days nMin. p50 Max. Min. p50 Max. ES 2015 0 - - - 1.4 7.7 30.8 FR 56 47.8 53.9 66.2 1.2 7.4 38.8 HR 24 43.7 49.2 64.0 0.7 10.7 39.6 IT 10 43.4 46.7 57.1 1.0 9.0 34.6 PT 5 75.7 78.1 78.1 1.4 7.7 34.8 ES 2016 2 78.5 78.5 78.5 2.2 5.9 30.1 FR 30 48.5 50.9 54.9 1.6 6.4 35.6 HR 21 51.5 55.8 66.0 1.3 10.1 39.8 IT 11 47.8 51.3 55.9 1.7 8.5 31.7 PT 4 50.7 51.3 51.9 1.4 6.7 31.2 TR 48 44.7 53.6 78.5 4.3 14.6 40.3 ES 2017 0 - - - 1.0 4.8 24.3 FR 42 45.1 51.6 61.0 1.9 6.4 39.7 HR 36 43.6 51.7 72.9 0.3 7.8 42.2 PT 7 48.9 56.6 64.4 2.2 9.4 39.9 TR 71 44.1 50.3 65.1 4.8 19.1 42.6 ES 2018 1 49.8 49.8 49.8 2.0 6.4 29.9 FR 12 47.3 47.8 49.3 1.2 6.4 34.8 HR 32 43.5 51.4 79.3 0.7 10.6 40.1 IT 29 43.6 47.6 61.9 1.7 12.0 36.6 PT 1 43.0 43.0 43.0 0.8 5.1 32.1 TR 60 43.7 47.4 63.8 3.1 14.4 42.5 2.2. Time Series Clustering In this study, HMM were used to detect clusters in PM 10 and PM 2.5 TS data. It represents a model-based strategy for clustering by assuming that each cluster of data is described by a different probability distribution (components of the mixture). This methodology is explained in Gómez-Losada et al. [ 14 , 20 ] and illustrated graphically in Figure 1. A complete definition and elements of HMM are provided in SM1. HMM represents a flexible method of modelling TS that exhibits dependence over time as those obtained in air quality monitoring networks. The model estimation lies in finding the parameter of the HMM that is most likely to have generated a given TS. This is referred Atmosphere 2021,12, 5 5 of 18 to as the maximum likelihood estimation problem, which in this study was solved using the Expectation-Maximisation algorithm [ 21 , 22 ]. The computational implementation of HMM in this study was accomplished using the depmixS4 package [ 23 ] from open-source software R[24]. The basic script of this implementation is provided in SM2. Atmosphere 2021, 12, x FOR PEER REVIEW 5 of 19 the HMM that is most likely to have generated a given TS. This is referred to as the maximum likelihood estimation problem, which in this study was solved using the Expectation-Maximisation algorithm [21,22]. The computational implementation of HMM in this study was accomplished using the depmixS4 package [23] from open-source software R [24]. The basic script of this implementation is provided in SM2. Figure 1. Time series (TS) clustering method showing the four pollution profiles (regimes) detected in the PM10 TS from the ES1517 Spanish station (2017). First to fourth pollution regimes are represented in blue, green, orange, and red, respectively. (A) TS observations (grey). Below TS, each observation is associated with a cluster. Separately, this grouping is represented above the TS. Horizontal lines represent the average values of each cluster (µ1, µ2, µ3, µ4). (B) The estimated density of the mixture of distributions (black line) is superimposed onto the histogram of TS data (grey line). Components of the mixture are represented for each cluster (coloured shadows). (C) Distribution of values of each TS pollution profile is represented as box-whisker plots. For many practical applications, there is often some physical significance attached to the hidden states of the HMM. For instance, in economics, states of the HMM can be related to expansion and recession periods and the interest is in studying the dynamics between them [25]; in developmental psychology, the states of the HMM are used to quantify knowledge that subjects express in an implicit learning task [26]. In this work, clusters of TS data are assigned to represented different PM concentration ranges (regimes of concentration or pollution profiles) in each modelled annual TS. HMM provides, inter alia, the mean and standard deviation values for each regime. These values represent the different parameterisation of the Gaussian distributions assigned to the regimes, providing Figure 1. Time series (TS) clustering method showing the four pollution profiles (regimes) detected in the PM 10 TS from the ES1517 Spanish station (2017). First to fourth pollution regimes are represented in blue, green, orange, and red, respectively. ( A ) TS observations (grey). Below TS, each observation is associated with a cluster. Separately, this grouping is represented above the TS. Horizontal lines represent the average values of each cluster ( µ1 , µ2 , µ3 , µ4 ). ( B ) The estimated density of the mixture of distributions (black line) is superimposed onto the histogram of TS data (grey line). Components of the mixture are represented for each cluster (coloured shadows). ( C ) Distribution of values of each TS pollution profile is represented as box-whisker plots. For many practical applications, there is often some physical significance attached to the hidden states of the HMM. For instance, in economics, states of the HMM can be related to expansion and recession periods and the interest is in studying the dynamics between them [ 25 ]; in developmental psychology, the states of the HMM are used to quantify knowledge that subjects express in an implicit learning task [ 26 ]. In this work, clusters of TS data are assigned to represented different PM concentration ranges (regimes of concentration or pollution profiles) in each modelled annual TS. HMM provides, inter alia, the mean and standard deviation values for each regime. These values represent the different parameterisation of the Gaussian distributions assigned to the regimes, providing very useful information as they fully characterise every analysed TS for a given time period Atmosphere 2021,12, 5 6 of 18 and monitoring site. Each cluster’s average values were represented as µi (i= 1, . . . ,n) with nthe number of detected clusters in the TS data (Figure 1A). 2.3. Time Series Clustering Once the pollution profiles in the TS data are determined, a meaning for each profile must be then assigned to estimate the PM 10 and PM 2.5 contributions from desert outbreaks. The following definitions are provided for the average values of each profile [ 14 ], applied to the remote stations studied in this work and considering as an example, the PM 10 TS from Figure 1: •µ1 (16 µ g · m −3 ) represents the underlying or threshold (background) concentration over which great changes in value are not expected over the years if atmospheric and pollution conditions remain relatively constant, being a characteristic of the studied area. Daily average PM 10 concentrations assigned to this regime are supposedly not caused by any direct influence of natural or anthropogenic sources, or if they are, they are negligible. •µ2 (22 µ g · m −3 ) is the average PM 10 concentration on the days affected by moderate contributions of anthropogenic sources due to activities that take place in the region. The value of µ2 is subject to slightly more variation than µ1 between years. The referenced days may be affected by contributions from natural sources attributable to African dust transport episodes that have a minor impact on the observed PM 10 concentrations. •µ3 (30 µ g · m −3 ) is the average PM 10 concentration on days affected by characteristic, usual contributions from African outbreaks. These contributions are highly variable in concentration and are the main factor responsible for the exceedances of the 50 µg·m−3 limit value established by the Directive 2008/50/EC on ambient air quality and cleaner air for Europe (Directive, 2008) for this pollutant. •µ4 (40 µ g · m −3 ) is the average PM 10 concentration on days with unusual but severe episodes of natural contributions from North African episodes. Some useful quantities can be then roughly estimated using these definitions, namely: •µ2 – µ1 (6 µ g · m −3 ): average concentration due to anthropogenic contributions from the region. •µ3 – µ2 (8 µ g · m −3 ): average concentration associated with characteristic contributions from dry regions from North Africa when they occur. •µ4 – µ2 (18 µ g · m −3 ): average concentration from severe contributions from dry regions from North Africa when they occur. It must be noted that the above definitions coincide from a conceptual point of view with those of horizontal profiles given by Lenschow et al. [ 27 ] when studying the PM 10 regimes of concentrations in the Berlin area. These definitions have to be adapted then to the geographical location under study and adapted to each monitoring network. Given the transnational scope of this study and it being unfeasible to characterise every PM regime for each station and country, these definitions were adopted and used generically. Furthermore, it must also be considered that contributions from other sources of PM such as those from wildfires or other anthropogenic sources of pollution could also be represented in the aforementioned defined pollution profiles. In this study, the above definitions are extended to the PM2.5 case. As with other available modelling techniques deriving information about pollution sources and the amount they contribute to ambient air pollution levels (e.g., chemical mass balance and positive matrix factorisation models), the time series clustering method is also subject to intrinsic limitations. Among them are that it provides estimates together with its uncertainties only based on the observed concentrations of pollutants (in this study, PM 10 and PM 2.5 ). To overcome this setback, definitions given to each of the pollution regimes were incorporated into this method when created (Gómez-Losada et al., 2015) to provide a reasonable physical explanation to these regimes. However, even though these Atmosphere 2021,12, 5 7 of 18 definitions work well in general, it can ignore some local or regional peculiarities from areas where emission sources originate. In this study, this limitation was partially overcome by selecting only background monitoring stations from the studied countries when analysing the contributions from North African arid regions. Another understandable disadvantage of the time series clustering method is that it does not take into consideration the origin or location of the emission sources, as air mass trajectory or wind-based models can do. For that reason, only countries bordering the Mediterranean Sea were selected for this study (Portugal was included due to being strongly influenced by Saharan dust) due to their proximity to the North African regions. Having exposed all of the above, it is also convenient to remember their advantages. As a modelling based on a statistical clustering technique, it does not require from expensive and time-consuming chemical assays, modelling outputs are comprehensive and easily interpreted, its computational implementation requires limited competence and low resources, and can be deployed using open-source software (R or Python programming languages). Moreover, scripts for the computational implementation have been provided by the authors on several occasions in other studies, and again, are available here upon request. Furthermore, the time clustering method can be used to complement other source apportionment methods, validate their modelling outputs, or simply used as an exploratory technique for deriving information from pollution sources. 2.4. Principal Component Analysis To study the clustering of remote monitoring sites according to the annual average values of their pollution profiles ( µ1 , µ2 , µ3 , µ4 ) for PM 10 and PM 2.5 , a biplot principal component analysis (PCA) was performed for each year. To facilitate the analysis, 95% confidence interval ellipses were calculated by studied country. The computational implementation of PCA was performed using the ggplot2 function from the ggbiplot library from R [ 28 ]. Datasets containing the estimated pollution profiles by year are available upon request. 3. Results and Discussion 3.1. Behaviour of PM2.5/PM10 Ratio PM 2.5 and PM 10 have several sources, thus presenting different physical and chemical properties. The study of PM 2.5 /PM 10 ratio is a useful tool to identify PM sources and effects on human health [ 10 , 11 ]. Evolution of the PM 2.5 /PM 10 ratios from 2015 to 2018 is illustrated in Figure 2. Remarkably, the pattern of the ratio in countries is similar to alternate years (2015 and 2017, 2016 and 2018). Even this analysis is absent for Turkey and Italy in 2015 and 2017, respectively, due to the lack of data fulfilling the quality criteria applied in this study; the highest ratios are obtained regularly for Croatia (HR—50% of the ratios were higher than 0.69 in 2016), followed by Italy (IT—with a median of 0.61 in 2016 and 2018), and the lowest ones for Portugal (PT—50% of the ratios were lower than 0.36 in 2016 and 2018) and Spain (ES—with a median of 0.5 in 2015 and 2017). The achieved values are in agreement with others already reported in the literature. For instance, Adães and Pires [ 9 ] analysed the PM 2.5 concentrations measured at 23 stations distributed in 12 different European cities. Median ratio values were between 0.31 and 0.73 with high values attributed to sites with high traffic influences. Xu et al. [ 10 ] analysed the temporal and spatial variability of PM 2.5 /PM 10 ratio. Mean values of 0.62 and 0.68 were determined for urban and background sites, respectively. Munir [ 11 ] analysed PM 2.5 /PM 10 ratio at 46 monitoring sites in the UK, where medium values were between 0.4 and 0.8. Atmosphere 2021,12, 5 8 of 18 Atmosphere 2021, 12, x FOR PEER REVIEW 8 of 19 Figure 2. Distribution of the general evolution of PM 2.5 /PM 10 ratios in studied countries as box-whisker plots (2015–2018). As a reference, medians between countries are joined, and a dotted horizontal line is added at 0.5 ratio. Figure 3 shows, as a more detailed analysis, a monthly analysis of PM 2.5 /PM 10 ratio. Higher values (with median higher than 0.65) are indicated with a red median as a reference in the box-whisker plots. As a general rule, higher ratios are observed in the coldest months (October to February), in particular for countries like Croatia (HR) and Italy (IT). Other countries like France (FR) may show the highest ratios in February. This seasonal pattern is justified by the increase in fuel consumption in this period of the year for domestic and industrial heating [29]. Additionally, the stable atmospheric conditions promote the dry deposition of coarse particles, maintaining the fine particles in the air [10,30]. It can be slightly appreciated that PM 2.5 /PM 10 ratio values higher than the P75 and lower than the P25 are prone to extending to the maximum and minimum ratios, respectively, in the warmest month. Figure 2. Distribution of the general evolution of PM 2.5 /PM 10 ratios in studied countries as box-whisker plots (2015–2018). As a reference, medians between countries are joined, and a dotted horizontal line is added at 0.5 ratio. Figure 3shows, as a more detailed analysis, a monthly analysis of PM 2.5 /PM 10 ratio. Higher values (with median higher than 0.65) are indicated with a red median as a reference in the box-whisker plots. As a general rule, higher ratios are observed in the coldest months (October to February), in particular for countries like Croatia (HR) and Italy (IT). Other countries like France (FR) may show the highest ratios in February. This seasonal pattern is justified by the increase in fuel consumption in this period of the year for domestic and industrial heating [ 29 ]. Additionally, the stable atmospheric conditions promote the dry deposition of coarse particles, maintaining the fine particles in the air [ 10 , 30 ]. It can be slightly appreciated that PM 2.5 /PM 10 ratio values higher than the P75 and lower than the P25 are prone to extending to the maximum and minimum ratios, respectively, in the warmest month. 3.2. Study of Pollution Profiles over Time Every PM 10 and PM 2.5 TS from 2015 to 2018 obtained from the studied countries (Table 1) was analysed independently with HMM to estimate their regimes of concentration. Regimes of monitoring stations by countries are illustrated as box-whisker plots in Figure 4 , for PM 10 (A) and PM 2.5 (B). The first to fourth pollution profiles are shown in blue, green, orange, and red, respectively. To detect the most representative particulate matter concentrations (PM 10 and PM 2.5 ), outliers (1.5 times below and higher the 25th and 75th percentile in distributions) were not considered in this study since the aim is to estimate the general trend of PM contributions quantitatively and not isolated outbreaks characterised by the highest concentrations. In Figure 4A, the highest pollution profile (red box-whisker plots) is estimated in Turkey (TR) and Spain (ES), while remaining with similar concentrations for the rest of countries (approximately around 50 µ g · m −3 ), except for Italy (IT) in 2016 and Portugal (PT) in 2016 and 2017. Bosnia and Herzegovina (BA) lacks representativeness regarding its pollution profiles estimation since only one TS from a single monitoring station was available for this country for 2016 and 2017. The rest of the pollution profiles (one to three) remains similar to little variations for all countries. Atmosphere 2021,12, 5 9 of 18 Atmosphere 2021, 12, x FOR PEER REVIEW 9 of 19 Figure 3. Monthly distribution of PM 2.5 /PM 10 ratios as in Figure 2. As a reference, countries exceeding the 0.65 ratio in 50% of the data are indicated by a red median, otherwise in blue. ES—Spain; FR—France; HR—Croatia; IT—Italy; PT—Portugal; TR—Turkey. 3.2. Study of Pollution Profiles over Time Every PM 10 and PM 2.5 TS from 2015 to 2018 obtained from the studied countries (Table 1) was analysed independently with HMM to estimate their regimes of concentration. Regimes of monitoring stations by countries are illustrated as box-whisker plots in Figure 4, for PM 10 (A) and PM 2.5 (B). The first to fourth pollution profiles are shown in blue, green, orange, and red, respectively. To detect the most representative particulate matter concentrations (PM 10 and PM 2.5 ), outliers (1.5 times below and higher the 25th and 75th percentile in distributions) were not considered in this study since the aim is to estimate the general trend of PM contributions quantitatively and not isolated outbreaks characterised Figure 3. Monthly distribution of PM 2.5 /PM 10 ratios as in Figure 2. As a reference, countries exceeding the 0.65 ratio in 50% of the data are indicated by a red median, otherwise in blue. ES—Spain; FR—France; HR—Croatia; IT—Italy; PT—Portugal; TR— Turkey. Atmosphere 2021,12, 5 16 of 18 Atmosphere 2021, 12, x FOR PEER REVIEW 17 of 19 Figure 9. PCA biplot analysis for PM 2.5 (2015–2018), showing the relative position of monitoring stations by countries concerning the principal components. BA—Bosnia and Herzegovina; ES— Spain; FR—France; HR—Croatia; IT—Italy; PT—Portugal; TR—Turkey. Table 4. Proportion of variance explained (in %) by pollutant and principal components (PC1 and PC2). Year PM 10 PM 2.5 PC1 PC2 PC1 PC2 2015 90.7 7.8 92.2 4.7 2016 80.4 16.1 90.3 7.1 2017 85.5 11.4 87.9 9.2 2018 88.4 10.0 90.9 8.0 4. Conclusions This study aims to estimate the PM 10 and PM 2.5 pollution profiles at background locations in Mediterranean countries and to quantify the contributions to PM 10 and PM 2.5 caused by dust outbreaks from North African desert regions. The time series clustering method based on Hidden Markov Models was used to estimate such pollution profiles, and contributions were obtained according to their given definitions. A discussion on the advantages and limitations of the time series clustering methods was provided in the methodological section. Countries like Turkey and Spain show the PM 10 profiles with the highest concentration, and Turkey, Croatia, Portugal, Italy, and France with respect to Figure 9. PCA biplot analysis for PM 2.5 (2015–2018), showing the relative position of monitoring stations by countries concerning the principal components. BA—Bosnia and Herzegovina; ES—Spain; FR—France; HR—Croatia; IT—Italy; PT—Portugal; TR—Turkey. Table 4. Proportion of variance explained (in %) by pollutant and principal components (PC1 and PC2). Year PM10 PM2.5 PC1 PC2 PC1 PC2 2015 90.7 7.8 92.2 4.7 2016 80.4 16.1 90.3 7.1 2017 85.5 11.4 87.9 9.2 2018 88.4 10.0 90.9 8.0 4. Conclusions This study aims to estimate the PM 10 and PM 2.5 pollution profiles at background locations in Mediterranean countries and to quantify the contributions to PM 10 and PM 2.5 caused by dust outbreaks from North African desert regions. The time series clustering method based on Hidden Markov Models was used to estimate such pollution profiles, and contributions were obtained according to their given definitions. A discussion on the advantages and limitations of the time series clustering methods was provided in the methodological section. Countries like Turkey and Spain show the PM 10 profiles with the highest concentration, and Turkey, Croatia, Portugal, Italy, and France with respect to PM 2.5 . Atmosphere 2021,12, 5 17 of 18 From 2015 to 2018, average characteristic contributions to PM 10 ( µ3 – µ2 ) were estimated as 11.6 ± 10.3 µ g · m −3 (Bosnia and Herzegovina), 8.8 ± 7.5 µ g · m −3 (Spain), 7.0 ± 6.2 µ g · m −3 (France), 8.1 ± 5.9 µ g · m −3 (Croatia), 7.5 ± 5.5 µ g · m −3 (Italy), 8.1 ±7.0 µg·m−3 (Portugal), and 17.0 ± 9.8 µ g · m −3 (Turkey). For the same period and PM 2.5 , the estimated average characteristic contributions were 4.1 ± 3.5 µ g · m −3 (Spain), 6.0 ± 4.8 µ g · m −3 (France), 9.1 ±6.4 µg·m−3 (Croatia), 5.2 ± 3.8 µ g · m −3 (Italy), 6.0 ± 4.4 µ g · m −3 (Portugal), and 9.0 ±5.6 µg·m−3 (Turkey). The analysed evolution of the PM 2.5 /PM 10 ratios in background stations from 2015 to 2018 was in agreement with other studies in the literature. Later, principal component analysis grouped the monitoring stations according to the different estimated PM10 and PM2.5 pollution profiles. Supplementary Materials: The description of the applied methods and the list of monitoring stations are available online at https://www.mdpi.com/2073-4433/12/1/5/s1. Author Contributions: Conceptualisation, Á.G.-L.; methodology, Á.G.-L. and J.C.M.P.; software, Á.G.-L.; validation, Á.G.-L. and J.C.M.P.; formal analysis, Á.G.-L.; investigation, Á.G.-L. and J.C.M.P.; resources, Á.G.-L.; data curation, Á.G.-L.; writing—original draft preparation, Á.G.-L. and J.C.M.P.; writing—review and editing, Á.G.-L. and J.C.M.P.; visualization, Á.G.-L.; supervision, Á.G.-L. and J.C.M.P.; project administration, Á.G.-L. and J.C.M.P.; funding acquisition, Á.G.-L. and J.C.M.P. All authors have read and agreed to the published version of the manuscript. Funding: This work was financially supported by Base Funding—UIDB/00511/2020 of the Laboratory for Process Engineering, Environment, Biotechnology and Energy—LEPABE—funded by national funds through the FCT/MCTES (PIDDAC). J.C.M. Pires acknowledges the FCT Investigator 2015 Programme (IF/01341/2015). Acknowledgments: The authors thank the European Environment Agency for providing the air quality data used in this study. Conflicts of Interest: The authors declare no conflict of interest. References 1. WHO. 9Out of10PeopleWorldwide Breathe PollutedAir, But More Countries are Taking Action. Available online: https://www.who. int/news-room/detail/02-05-2018-9-out-of-10-people-worldwide-breathe-polluted-air-but-more-countries-are-taking-action (accessed on 19 November 2020). 2. Guerreiro, C.B.B.; Foltescu, V.; de Leeuw, F. Air quality status and trends in europe. Atmos. Environ. 2014 ,98, 376–384. [CrossRef] 3. Juda-Rezler, K.; Reizer, M.; Maciejewska, K.; Blaszczak, B.; Klejnowski, K. Characterization of atmospheric pm(2.5 )sources at a central european urban background site. Sci. Total Environ. 2020,713, 136729. [CrossRef] [PubMed] 4. Yang, W.; Chen, H.; Wu, J.; Wang, W.; Zheng, J.; Chen, D.; Li, J.; Tang, X.; Wang, Z.; Zhu, L.; et al. Characteristics of the source apportionment of primary and secondary inorganic pm2.5 in the pearl river delta region during 2015 by numerical modeling. Environ. Pollut. 2020,267, 115418. [CrossRef] [PubMed] 5. Duan, L.; Yan, L.; Xiu, G. Online measurement of pm2.5 at an air monitoring supersite in yangtze river delta: Temporal variation and source identification. Atmosphere 2020,11, 789. [CrossRef] 6. Eeftens, M.; Tsai, M.Y.; Ampe, C.; Anwander, B.; Beelen, R.; Bellander, T.; Cesaroni, G.; Cirach, M.; Cyrys, J.; de Hoogh, K.; et al. Spatial variation of pm2.5, pm10, pm2.5 absorbance and pmcoarse concentrations between and within 20 european study areas and the relationship with no2—results of the escape project. Atmos. Environ. 2012,62, 303–317. [CrossRef] 7. Neophytou, A.M.; Yiallouros, P.; Coull, B.A.; Kleanthous, S.; Pavlou, P.; Pashiardis, S.; Dockery, D.W.; Koutrakis, P.; Laden, F. Particulate matter concentrations during desert dust outbreaks and daily mortality in nicosia, cyprus. J. Expo. Sci. Environ. Epid. 2013,23, 275–280. [CrossRef] 8. Green, M.; Xu, J. Causes of haze in the columbia river gorge. J. Air Waste Manag. 2007,57, 947–958. [CrossRef] 9. Adaes, J.; Pires, J.C.M. Analysis and modelling of pm2.5 temporal and spatial behaviors in european cities. Sustainability 2019 , 11, 6019. [CrossRef] 10. Xu, G.; Jiao, L.M.; Zhang, B.E.; Zhao, S.L.; Yuan, M.; Gu, Y.Y.; Liu, J.F.; Tang, X. Spatial and temporal variability of the pm2.5/pm10 ratio in wuhan, central china. Aerosol. Air Qual. Res. 2017,17, 741–751. [CrossRef] 11. Munir, S. Analysing temporal trends in the ratios of pm2.5/pm10 in the uk. Aerosol. Air Qual. Res. 2017,17, 34–48. [CrossRef] 12. Sugimoto, N.; Shimizu, A.; Matsui, I.; Nishikawa, M. A method for estimating the fraction of mineral dust in particulate matter using pm2.5-to-pm10 ratios. Particuology 2016,28, 114–120. [CrossRef] 13. Speranza, A.; Caggiano, R.; Margiotta, S.; Trippetta, S. A novel approach to comparing simultaneous size-segregated particulate matter (pm) concentration ratios by means of a dedicated triangular diagram using the agri valley pm measurements as an example. Nat. Hazard. Earth Syst. 2014,14, 2727–2733. [CrossRef] Atmosphere 2021,12, 5 18 of 18 14. Gomez-Losada, A.; Pires, J.C.M.; Pino-Mejias, R. Time series clustering for estimating particulate matter contributions and its use in quantifying impacts from deserts. Atmos. Environ. 2015,117, 271–281. [CrossRef] 15. Escudero, M.; Querol, X.; Pey, J.; Alastuey, A.; Perez, N.; Ferreira, F.; Alonso, S.; Rodriguez, S.; Cuevas, E. A methodology for the quantification of the net african dust load in air quality monitoring networks. Atmos. Environ. 2007,41, 5516–5524. [CrossRef] 16. Stafoggia, M.; Zauli-Sajani, S.; Pey, J.; Samoli, E.; Alessandrini, E.; Basagana, X.; Cernigliaro, A.; Chiusolo, M.; Demaria, M.; Diaz, J.; et al. Desert dust outbreaks in southern europe: Contribution to daily pm10 concentrations and short-term associations with mortality and hospital admissions. Environ. Health Persp. 2016,124, 413–419. [CrossRef] 17. Gama, C.; Pio, C.; Monteiro, A.; Russo, M.; Fernandes, A.P.; Borrego, C.; Baldasano, J.M.; Tchepel, O. Comparison of methodologies for assessing desert dust contribution to regional pm10 and pm2.5 levels: A one-year study over portugal. Atmosphere 2020 , 11, 134. [CrossRef] 18. Goudie, A.S. Desert dust and human health disorders. Environ. Int. 2014,63, 101–113. [CrossRef] [PubMed] 19. Querol, X.; Tobias, A.; Perez, N.; Karanasiou, A.; Amato, F.; Stafoggia, M.; Garcia-Pando, C.P.; Ginoux, P.; Forastiere, F.; Gumy, S.; et al. Monitoring the impact of desert dust outbreaks for air quality for health studies. Environ. Int. 2019 ,130, 104867. [CrossRef] 20. Gomez-Losada, A.; Pires, J.C.M.; Pino-Mejias, R. Modelling background air pollution exposure in urban environments: Implications for epidemiological research. Environ. Modell Softw. 2018,106, 13–21. [CrossRef] 21. Dempster, A.P.; Laird, N.M.; Rubin, D.B. Maximum likelihood from incomplete data via the em algorithm. J. R. Stat. Soc. Ser. B 1977,39, 1–38. 22. Wu, X.D.; Kumar, V.; Quinlan, J.R.; Ghosh, J.; Yang, Q.; Motoda, H.; McLachlan, G.J.; Ng, A.; Liu, B.; Yu, P.S.; et al. Top 10 algorithms in data mining. Knowl. Inf. Syst. 2008,14, 1–37. [CrossRef] 23. Visser, I.; Speekenbrink, M. Depmixs4: An r package for hidden markov models. J. Stat. Softw. 2010,36, 1–21. [CrossRef] 24. R: A Language and Environment for Statistical Computing. Available online: http://www.R-project.org/ (accessed on 5 October 2020). 25. Dias, J.G.; Vermunt, J.K.; Ramos, S. Mixture Hidden Markov Models in Finance Research, Advances in Data Analysis, Data Handling and Business Intelligence; Fink, A., Lausen, B., Seidel, W., Ultsch, A., Eds.; Springer: Berlin/Heidelberg, Germany, 2010; pp. 451–459. 26. Visser, I.; Raijmakers, M.E.J.; Molenaar, P.C.M. Fitting hidden markov models to psychological data. Sci. Program. 2002 , 10, 185–199. [CrossRef] 27. Lenschow, P.; Abraham, H.J.; Kutzner, K.; Lutz, M.; Preuss, J.D.; Reichenbacher, W. Some ideas about the sources of pm10. Atmos. Environ. 2001,35, S23–S33. [CrossRef] 28. Vu, V.Q. A ggplot2 Based Biplot. Available online: http://github.com/vqv/ggbiplot (accessed on 5 October 2020). 29. Liu, J.; Mauzerall, D.L.; Chen, Q.; Zhang, Q.; Song, Y.; Peng, W.; Klimont, Z.; Qiu, X.H.; Zhang, S.Q.; Hu, M.; et al. Air pollutant emissions from chinese households: A major and underappreciated ambient pollution source. Proc. Natl. Acad. Sci. USA 2016 , 113, 7756–7761. [CrossRef] [PubMed] 30. Huang, W.; Long, E.S.; Wang, J.; Huang, R.Y.; Ma, L. Characterising spatial distribution and temporal variation of pm10 and pm2.5 mass concentrations in an urban area of southwest china. Atmos. Pollut. Res. 2015,6, 842–848. [CrossRef] 31. Querol, X.; Perez, N.; Reche, C.; Ealo, M.; Ripoll, A.; Tur, J.; Pandolfi, M.; Pey, J.; Salvador, P.; Moreno, T.; et al. African dust and air quality over spain: Is it only dust that matters? Sci. Total Environ. 2019,686, 737–752. [CrossRef] 32. Perez, L.; Tobias, A.; Querol, X.; Kunzli, N.; Pey, J.; Alastuey, A.; Viana, M.; Valero, N.; Gonzalez-Cabre, M.; Sunyer, J. Coarse particles from saharan dust and daily mortality. Epidemiology 2008,19, 800–807. [CrossRef] 33. Salvador, P.; Alonso-Perez, S.; Pey, J.; Artinano, B.; de Bustos, J.J.; Alastuey, A.; Querol, X. African dust outbreaks over the western mediterranean basin: 11-year characterisation of atmospheric circulation patterns and dust source areas. Atmos. Chem. Phys. 2014 , 14, 6759–6775. [CrossRef] 34. Pérez, N.; Querol, X.; Alastuey, A.; Olivares, I.; Campos, A.; Hervás, M.; Cornide, M.J.; Javato, R.; Salvador, P.; Artínano, B.; et al. Episodios Naturales de Partículas 2018; Consejo Superior de Investigaciones Científicas (CSIC), Ministerio para la Transición Ecológica: Madrid, Spain, 2019. 35. Pérez, N.; Querol, X.; Alastuey, A.; Orío, A.; Olivares, I.; Reina, F.; Hervás, M.; Cornide, M.J.; Javato, R.; Salvador, P.; et al. Episodios Naturales de Partículas 2016; CSIC, CIEMAT, Ministerio de Agricultura y Pesca, Alimentación y Medio Ambiente: Madrid, Spain, 2017. 36. Pérez, N.; Querol, X.; Alastuey, A.; Orío, A.; Olivares, I.; Reina, F.; Hervás, M.; Cornide, M.J.; Javato, R.; Salvador, P.; et al. Episodios Naturales de Partículas 2017; CSIC, CIEMAT, Ministerio de Agricultura y Pesca, Alimentación y Medio Ambiente: Madrid, Spain, 2018. 37. Pérez, N.; Querol, X.; Alastuey, A.; Orío, A.; Reina, F.; Pallarés, M.; Salvador, P.; Artínano, B.; de la Rosa, J. Episodios Naturales de Partículas 2015; CSIC, CIEMAT, Ministerio de Agricultura, Alimentación y Medio Ambiente: Madrid, Spain, 2016. 38. Wang, Q.; Gu, J.; Wang, X. The impact of Sahara dust on air quality and public health in European countries. Atmos. Environ. 2020,241, 117771. [CrossRef] 39. Sun, P.; Nie, W.; Chi, X.; Xie, Y.; Huang, X.; Xu, Z.; Qi, X.; Xu, Z.; Wang, L.; Wang, T.; et al. Two years of online measurement of fine particulate nitrate in the western Yangtze River Delta: Influences of thermodynamics and N2O5 hydrolysis. Atmos. Chem. Phys. 2018,18, 17177–17190. [CrossRef]