Meteorology-driven variability of air pollution (PM1) revealed with explainable machine learning

Air pollution, in particular high concentrations of particulate matter smaller than 1μm in diameter (PM1), continues to be a major health problem, and meteorology is known to substantially influence atmospheric PM concentrations. However, the scientific understanding of the ways by which complex interactions of meteorological factors lead to high pollution episodes is inconclusive. In this study, a novel, data-driven approach based on empirical relationships is used to characterise, and better understand the meteorology-driven component of PM1 variability. A tree-based machine learning model is set up 5 to reproduce concentrations of speciated PM1 at a suburban site southwest of Paris, France, using meteorological variables as input features. The model is able to capture the majority of occurring variance of mean afternoon total PM1 concentrations (coefficient of determination (R) of 0.58), with model performance depending on the individual PM1 species predicted. Based on the models, an isolation and quantification of individual, season-specific meteorological influences for process understanding at the measurement site is achieved using SHapley Additive exPlanation (SHAP) regression values. Model results suggest 10 that winter pollution episodes are often driven by a combination of shallow mixed layer heights (MLH), low temperatures, low wind speeds or inflow from northeastern wind directions. Contributions of MLHs to the winter pollution episodes are quantified to be on average ∼5μg/m for MLHs below <500 m agl. Temperatures below freezing initiate formation processes and increase local emissions related to residential heating, amounting to a contribution to predicted PM1 concentrations of as much as ∼9μg/m. Northeasterly winds are found to contribute ∼5μg/m to predicted PM1 concentrations (combined effects of u15 and v-wind components), by advecting particles from source regions, e.g. central Europe or the Paris region. Meteorological drivers of unusually high PM1 concentrations in summer are temperatures above ∼25 ◦C (contributions of up to ∼2.5μg/m), dry spells of several days (maximum contributions of ∼1.5μg/m) and wind speeds below ∼2 m/s (maximum contributions of ∼3μg/m), which cause a lack of dispersion. High-resolution case studies are conducted showing a large variability of processes that can lead to high pollution episodes.The identification of these meteorological conditions that increase air pollution 20 could help policy makers to adapt policy measures, issue warnings to the public, or to assess the effectiveness of air pollution measures.


Introduction
Air pollution has serious implications on human well-being, including deleterious effects on the cardiovascular system and the 25 lungs (Hennig et al., 2018;Lelieveld et al., 2019), and an increased number of asthma seizures (Hughes et al., 2018). This includes particles smaller than 1 µm in diameter (PM 1 ), which are associated with fits of coughing (Yang et al., 2018) and an increase in emergency hospital visits (Chen et al., 2017b). The adverse health effect lead to an increase in mortality of people exposed to high particle concentrations (Samoli et al., 2008(Samoli et al., , 2013Lelieveld et al., 2015). People living in urban areas are particularly affected by poor air quality and with increasing urbanization, their number is projected to grow (Baklanov et al.,30 2016; Li et al., 2019). These developments have motivated several countermeasures to improve air quality. Proposed efforts to reduce anthropogenic particle emissions include partial traffic bans (Su et al., 2015;Dey et al., 2018) and the reduction of solid fuel use for domestic heating (Chafe et al., 2014). Although emissions play an important role for PM concentrations in the atmosphere, meteorological conditions related to large-scale circulation patterns as well as local-scale boundary layer processes and interactions with the land surface are major drivers of PM variability as well (Cermak and Knutti, 2009;Bressi 35 et al., 2013;Megaritis et al., 2014;Dupont et al., 2016;Petäjä et al., 2016;Yang et al., 2016;Li et al., 2017). Wind speed and direction generally have a strong influence on air quality as they determine the advection of pollutants (Petetin et al., 2014;Petit et al., 2015;Srivastava et al., 2018). Limiting the vertical exchange of air masses, the mixed layer height (MLH) governs the volume of air in which particles are typically dispersed. Although some authors indicate that mixed layer height cannot be related directly to concentrations of pollutants and that other meteorological parameters and local sources need to be considered 40 (Geiß et al., 2017), a lower MLH can increase PM concentrations as particles are not mixed into higher atmospheric levels and accumulate near the ground (Gupta and Christopher, 2009;Schäfer et al., 2012;Stirnberg et al., 2020).
Higher MLHs in combination with high wind speeds increase atmospheric ventilation processes, thus decreasing near-surface particle concentrations (Sujatha et al., 2016;Wang et al., 2018). Air temperature can influence PM concentrations in multiple ways, e.g. by modifying the emission of secondary PM precursors such as volatile organic compounds (VOCs) during summer 45 (Fowler et al., 2009;Megaritis et al., 2013;Churkina et al., 2017), and by condensating high saturation vapour pressure compounds such as nitric acid and sulfuric acid (Hueglin et al., 2005;Pay et al., 2012;Bressi et al., 2013;Megaritis et al., 2014). The wet removal of particles by precipitation is known to be an efficient atmospheric aerosol sink (Radke et al., 1980;Bressi et al., 2013), while moisture in the atmosphere can stimulate secondary particle formation processes (Ervens et al., 2011). Although all these atmospheric conditions and processes have been identified as drivers of local air quality, it is usually 50 a complex combination of meteorological and chemical processes that lead to the formation of high-pollution events (Petit et al., 2015;Dupont et al., 2016;Stirnberg et al., 2020).
The metropolitan area of Paris is one of the most densely populated and industrialised areas in Europe. Thus, air quality is a recurring issue and has been at the focus of many studies in the past years (Bressi et al., 2014;Petetin et al., 2014;Petit et al., 2015;Dupont et al., 2016;Petit et al., 2017;Srivastava et al., 2018). Results indicate that the Paris metropolitan region is 55 often affected by mid-to long-range transport of pollutants, as due to the city's flat orography, an efficient horizontal exchange of air masses is frequent (Bressi et al., 2013;Petit et al., 2015). High-pollution events commonly occur in late autumn, winter, and early spring. Often, these episodes are characterised by stagnant atmospheric conditions and a combination of local contributions, e.g. traffic emissions, residential emissions, or regionally transported particles, e.g. ammonium nitrates from manure spreading, or sulfates from point sources (Petetin et al., 2014;Petit et al., 2014Petit et al., , 2015Srivastava et al., 2018). High-60 pressure conditions with air masses originating from continental Europe (Belgium, Netherlands, West Germany) are generally associated with an increase in particle concentrations, especially of secondary inorganic aerosols (SIA, Bressi et al. (2013); Srivastava et al. (2018). The regional contribution has been found to be in the range of 70 % for background concentrations in Paris of particles with a diameter smaller 2.5 µm (Petetin et al., 2014). Hence the variability between high-pollution episodes in terms of timing, sources and meteorological boundary conditions is considerable (Petit et al., 2017). Previous approaches to 65 determine meteorological drivers of air pollution included, for example, the use of chemical transport models (CTMs), which, however, require comprehensive knowledge on emission sources and secondary particle formation pathways and are associated with considerable uncertainties (Sciare et al., 2010;Petetin et al., 2014;Kiesewetter et al., 2015). Further methods rely on data exploration, e.g. the statistical analysis of time-series (Dupont et al., 2016), which can be coupled with positive matrix factorization (PMF, Paatero and Tapper, 1994) to derive PM sources (Petit et al., 2014;Srivastava et al., 2018). To take into account 70 the interconnected nature of PM drivers, multivariate statistical approaches such as principal component analysis (PCA) have been applied (Chen et al., 2014;Leung et al., 2017). In recent years, machine learning techniques have been increasingly used to expand the analysis of PM concentrations with respect to meteorology, allowing to retrace general patterns (Hu et al., 2017;Grange et al., 2018).
Here, the multivariate and highly interconnected nature of meteorology-dependent atmospheric processes influencing local 75 PM 1 concentrations at a suburban site southwest of Paris is analysed in a data-driven way. Therefore, a state-of-the-art explainable machine learning model is set up to reproduce the variability of PM 1 concentrations, thereby capturing empirical relationships between PM 1 concentrations and meteorological parameters. The goal is to separate and quantify influences of the meteorological variables on PM 1 concentrations to advance the process understanding of the complex mechanisms that govern pollution concentrations at the measurement site. Localised (i.e. situation-based) and individualised attributions of fea-80 ture contributions are performed using SHapley Additive exPlanation regression (SHAP) values (Lundberg and Lee, 2017;Lundberg et al., 2018aLundberg et al., , 2020, allowing to infer meteorology-dependent processes driving PM concentrations at high temporal resolution. Typical situations that lead to high PM 1 concentrations are identified, serving as a decision support to policymakers to issue preventative warnings to the public if these situations are to be expected. In addition, by directly accounting for meteorological effects on PM1 concentrations, such a machine learning-based framework could help in assessing the effective-85 ness of measures towards better air quality. Furthermore, the proposed ML framework can be viewed as a first step towards a data-driven, prognostic tool in operational air quality forecasting, complementary to CTM approaches.

Data sets
Seven years (2012-2018) of meteorological and air quality data from the Site Instrumental de Recherche par Télédétection Atmosphérique (SIRTA, Haeffelin et al., 2005) supersite are the basis of this study. The SIRTA Atmospheric Observatory is 90 located about 25km southwest of Paris (48.713 • N and 2.208 • E, Fig. 1). This study focuses on day-to-day variations of total and speciated PM 1 , a highly health relevant fraction of PM including small particles that can penetrate deep into the lungs (Yang et al., 2018;Chen et al., 2017a). To separate diurnal effects e.g. the development of the boundary layer during morning hours (Petit et al., 2014;Dupont et al., 2016;Kotthaus and Grimmond, 2018a) from day-to-day variations of PM 1 , mean concentrations of total and speciated PM 1 for the afternoon period 12-15 UTC are considered, when the boundary layer is fully 95 developed. In sections 2.1 and 2.2, the PM 1 and meteorological data and preprocessing steps before setting up the machine learning model are described. The applied machine learning model and data analysis techniques are presented in sections 3.1 and 3.2.

Submicron particle measurements
Aerosol chemical speciation monitor (ACSM, Ng et al., 2011) measurements are conducted at SIRTA in the framework of the 100 ACTRIS project. The ACSM provides continuous and near real-time measurements of the major chemical composition of nonrefractory submicron aerosols, i.e., organics (Org), ammonium (NH + 4 ), sulfate (SO 2− 4 ), nitrate (NO − 3 ) and chloride (Cl − ). A detailed description of its functionality can be found in Ng et al. (2011). Data processing and validation protocol can be found in Petit et al. (2015) and Zhang et al. (2019). In addition, black carbon (BC) has been monitored by a seven-wavelength Magee Scientific Aethalometer AE31 from 2011 to mid-2013, and a dual-spot AE33 (Drinovec et al., 2015) from mid-2013 onwards.
Consistency of both instruments have been checked in Petit et al. (2014). Using the multispectral information, a differentiation into fossil fuel-based BC (BCff) and BC from wood burning (BCwb) is achieved (Sciare et al., 2010;Healy et al., 2012;Petit et al., 2014;Zhang et al., 2019). Here, the sum of all measured species is assumed to represent the total PM 1 content (see Petit et al., 2014Petit et al., , 2015. The consistency of ACSM and Aethalometer measurements is checked by comparing the sum of all monitored species with measurements of a nearby Tapered Element Oscillating Microbalance equipped with a Filter Dynamic 110 Measurement System (TEOM-FDMS). PM 1 measurements are representative of suburban background pollution levels of the region of Paris (Petit et al., 2015). As an additional input to the machine learning model, the average fraction of NO − 3 of the previous day is added (NO3_frac). Pollution events dominated by NO − 3 are often linked to regional-scale events, which depend on anthropogenically-influenced processes in the source regions of NO − 3 precursors (Petit et al., 2017). This is approximated by the inclusion of the average fraction of NO − 3 of the previous day, assuming that a high fraction of NO − 3 indicates the occurrence 115 of such an anthropogenically-influenced regime.

Meteorological data
Following the objective of this study, a set of meteorological variables is chosen as inputs for the ML model that either influence PM concentrations directly via dilution (MLH, wind speed (ws), and wet scavenging of particles (precipitation)) and particle transport (wind direction as u, v components, air pressure (AirPres)), as a proxy for emissions (e.g. from residential heating: 120 temperature at a height of 2 m (T)), and as a proxy for transformation processes (total incoming solar radiation (TISR), relative humidity (RH), T). Data are taken from the quality-controlled and 1h averaged re-analysed observation (ReObs) dataset.
Further information on the instrumentation used for the acquisition of these variables is provided in Chiriaco et al. (2018). counted and used as input to capture the particle accumulation effect between precipitation events (Rost et al., 2009;Petit et al., 2017).

Machine learning model: technique and application
Gradient Boosted Regression Trees (GBRT, used here in a python 3.6.4 environment with the scikit-learn module, Friedman,140 2002; Pedregosa et al., 2012) are applied to predict daily total and speciated PM 1 concentrations. As a tree-based method, GBRTs use a tree regressor, which sets up decision trees based on a training data set. The trees split the training data along decision nodes, creating homogeneous subsamples of the data by minimizing the variance of each subsample. For each subsample, regression trees fit the mean response of the model to the observations (Elith et al., 2008). To increase confidence in the model outputs, decision trees are combined to form an ensemble prediction. Trees are sequentially added to the ensemble (Elith 145 et al., 2008;Rybarczyk and Zalakeviciute, 2018) and each new tree is fitted to the predecessor's previous residual error using gradient descent (Friedman, 2002). This is an advantage of GBRT over standard ensemble tree methods (e.g. Random Forests (RF), Just et al., 2018) as trees are built systematically and fewer iterations are required (Elith et al., 2008). Characteristics of the meteorological training data set with respect to observed total and speciated PM 1 concentrations are conveyed to the statistical model. The learned relationships are then used for model interpretation and to produce estimates of PM 1 based on 150 unseen meteorological data to test the model. The architecture of the statistical model is determined by the hyperparameters, e.g. the number of trees, the maximum depth of each tree (i.e., the number of split nodes on each tree) and the learning rate (i.e., the magnitude of the contribution of each tree to the model outcome, which is basically the step size of the gradient descent).
The hyperparameters are tuned by executing a grid search, systematically validating testing previously defined hyperparameter combinations and determining the best combination via a three-fold cross validation. Note that PM 1 data is not uniformly 155 distributed, i.e. there is more data available for mid-range PM 1 concentrations. To avoid that the model primarily optimizes its predictions on these values, a least-squares loss function was chosen. This loss function is more sensitive to higher PM 1 values (i.e. outliers of the PM 1 data distribution), as it strongly penalises high absolute differences between predictions and observations. Accordingly, the model is adjusted to reproduce higher concentrations as well.
For each PM species, a specific GBRT model is set up and used for the analysis of meteorological influences on individual 160 PM 1 species (see section 4.2). Additionally, a quasi-total PM 1 model is used to reproduce the sum of all species at once, which is used for an analysis of meteorological drivers of high-pollution events (see sections 4.3 and 4.4). Train and test data sets to evaluate each model are created by randomly splitting the full data set. These splits, however, are the same for the species models and the full PM 1 model to ensure comparability between the models. Three quarters of the data are used for training and hyperparameter tuning with cross-validation (n=1086), and one quarter for testing (n=363). In addition, the robustness 165 of the model results is tested by repeating this process ten times, resulting in ten models with different train-/ test-splits and different hyperparameters.

Explaining model decisions to infer processes: SHapley Additive exPlanation (SHAP) values
While being powerful predictive models, tree-based machine learning methods also have a high interpretability (Lundberg et al., 2020). In order to understand physical mechanisms on the basis of model decisions, the contributions of the meteorological 170 input features to the model outcome are analysed. Feature contributions are attributed using SHAP values, which allow for an individualised, unique feature attribution for every prediction (Shapley, 1953;Lundberg and Lee, 2017;Lundberg et al., 2018aLundberg et al., , 2020. SHAP values provide a deeper understanding of model decisions than the relatively widely used partial dependence plots (Friedman, 2001;Goldstein et al., 2015;Fuchs et al., 2018;Lundberg et al., 2018a;McGovern et al., 2019;Stirnberg et al., 2020). Partial dependence plots show the global mean effect of an input feature to the model outcome, while SHAP 175 values quantify the feature contribution to each single model output, accounting for multicollinearity. Feature contributions are calculated from the difference in model outputs with that feature present, versus outputs for a retrained model, without the feature. Since the effect of withholding a feature depends on other features in the model due to interactive effects between the features, differences are computed for all possible feature subset combinations of each data instance (Lundberg and Lee, 2017). 4 Results and discussion

Model performance
The performance of the species and total PM 1 models, each with ten model iterations (of which each has different hyperpa-195 rameters) is assessed by comparing the coefficient of determination (R 2 ) and normalised root mean square error (NRSME) for the independent test data that was withheld during the training process (Fig. 3). While the models for BCwb, BCff and total PM 1 show small spread, Cl − and NO − 3 exhibit larger variations between model runs (indicated by horizontal and vertical lines in Fig. 3). This suggests that while drivers of variations in BCff concentration are well covered by the model, this is less so in the case of Cl − and NO  Depending on whether positive or negative SHAP values dominate, the prediction is higher or lower than the base value (Lundberg et al., 2018b). Adapted from https://github.com/slundberg/shap.
The mean input feature importance, ordered from high to low, of the total PM 1 model run by means of the SHAP feature attribution values is shown in Fig. 4, The NO − 3 fraction of the previous day has the highest impact on the model, followed by 205 temperature, wind direction information, and MLH. To some extent, NO − 3 fraction can be related to PM 1 mass concentrations (Petit et al., 2015;Beekmann et al., 2015). This means that the higher the PM 1 levels one day, the greater the chances of having higher PM 1 levels the next day See Fig. B1). Lower wind speeds generally lead to higher particle concentrations (see Fig. B2) due to a lack of dispersion (Sujatha2016). Temperature, MLH and wind direction require an in-depth analysis, as changes of these variables cause nonlinear responses in PM1 predictions, which vary also between species.

Influence of temperature
The impact of ambient air temperature on modelled species concentrations is highly non-linear (Fig. 5). All species show increased contributions to model outcomes at temperatures below ∼4 • C while the contribution of high temperatures on model 220 outcomes differs substantially between species. The statistical model is able to reproduce well-known characteristics of species concentration variations related to temperature. For example, sulfate formation is enhanced with increasing temperatures (Fig.   5d) due to an increased oxidation rate of SO 2 (see Dawson et al., 2007;Li et al., 2017) and strong solar irradiation due to photochemical oxidation (Gen et al., 2019). Dawson et al. (2007) reported an increase of 34 ng/m 3 K for PM 2.5 concentrations using a CTM. The increase in sulfate at low ambient temperatures as suggested by Fig. 5d is not reported in this study. It 225 is likely linked to increased aqueous phase particle formation in cold and foggy situations (Rengarajan et al., 2011;Petetin et al., 2014;Cheng et al., 2016). Considerable local formation of nitrate at low temperatures (Fig. 5b) is consistent with results from previous studies in western Europe and enhanced formation of ammonium nitrate at lower temperatures (Fig. 5c) by the shifting gas-particle equilibrium is a well-known pattern (e.g., Clegg et al., 1998;Pay et al., 2012;Bressi et al., 2013;Petetin et al., 2014;Petit et al., 2015). The increase in organic matter and BCwb concentrations at low temperatures (Fig. 5g) 230 is likely related to the emission intensity, as biomass burning is often used for domestic heating in the study area (Favez et al., 2009; Sciare et al., 2010;Healy et al., 2012;Jiang et al., 2019). In addition, organic matter concentrations are linked to the condensation of semi-volatile organic species at low temperatures (Putaud et al., 2004;Bressi et al., 2013). The sharp increase in modelled concentrations of organics above 25 • C (Fig 5a) could be due to enhanced biogenic activity leading to a rise in biogenic emissions and secondary aerosol formation (Guenther et al., 1993;Churkina et al., 2017;Jiang et al., 2019).

235
The contribution of temperature on modelled total PM 1 concentrations (Fig. 6h) is consistent with the response patterns to changes in temperatures described for the individual species in panels 6a-6g, with positive contributions at both low (<4 • C) and high air temperatures (>25 • C). For temperatures below freezing, the model allocates maximum contributions to modelled total PM 1 concentrations of up to 12 µg/m 3 . The spread of SHAP values between model iterations is generally higher for low temperatures (vertical grey bars in Figs 5-7), where SHAP values are of greater magnitude, but in all cases the signal contained 240 in the feature contributions far exceeds the spread between model runs.

Influence of the mixed layer height (MLH)
Variations in MLH can have substantial impact on near-surface particle concentrations, as the mixed layer is the atmospheric volume in which the particles are dispersed (see Klingner and Sähn, 2008;Dupont et al., 2016;Wagner and Schäfer, 2017).
The effect of MLH variations on modelled particle concentrations is highly nonlinear and varies in magnitude for all species ( Fig. 6). Similar to the patterns observed for temperature SHAP values, the inter-model variation of predictions is highest for low MLHs where predicted particle concentrations have the highest variation. For predicted total PM 1 concentrations, the maximum positive contribution of the MLH is as high as 5.5 µg/m 3 while negative contributions can amount to -2 µg/m 3 .
While the maximum influence of MLH is lower than the maximum influence determined for air temperature, the frequency of shallow MLH is far greater than that of the minimum temperatures that have the largest effect (Figs 5d & 6d). Contributions of 250 MLH to predicted particle concentrations are highest for very shallow mixed layers due to the accumulation of particles close to the ground (Dupont et al., 2016;Wagner and Schäfer, 2017). In addition to causing particles to accumulate near the surface, low MLH can also provide effective pathways for local new particle formation. Secondary pollutants, such as ammonium nitrate, are increased at low MLHs when conditions favorable to their formation usually coincide with reduced vertical mixing (i.e., low temperatures, often in combination with high RH, Pay et al., 2012;Petetin et al., 2014;Dupont et al., 2016;Wang 255 et al., 2016). BC concentrations, on the other hand, are dominated by primary emissions, as is a substantial fraction of organic matter (Petit et al., 2015). Hence, the accumulation of these particles during low buoyancy conditions can explain the strong influence of MLH on BCwb and BCff. A relatively distinct transition from positive contributions during shallow boundary layer conditions (∼0-800 m) towards negative contributions at high MLHs is evident for all species except SO 2− 4 . Modelled SO 2− 4 concentrations show a less distinct response to changes in MLH as they are largely driven by gaseous precursor sources and 260 particle advection, both rather independent of MLH (Pay et al., 2012;Petit et al., 2014Petit et al., , 2015, so that the accumulation effect is less important. The increase of SO 2− 4 concentrations with higher MLHs (>∼ 1500 m agl) could be linked to the effective transport of SO 2− 4 and its precursor SO 2 . In agreement with results from previous studies focusing on PM 10 ( Grange et al., 2018;Stirnberg et al., 2020) or PM 2.5 , SHAP values do not change much for MLH above ∼800-900 m, i.e. boundary layer height variations above 265 this level do not influence submicron particle concentrations. Positive contributions of MLHs above ∼800-900 m on predicted PM 1 concentrations, as visible in Fig. 6 for some species, have been previously reported by Grange et al. (2018), who relate this pattern to enhanced secondary aerosol formation in a very deep and dry boundary layer. The positive influence of high MLHs on species that are partly secondarily formed, e.g. SO 2− 4 and Org, could be explained following this argumentation. The increase in SHAP values observed for BCff at high MLHs could be also related to secondary aerosol formation processes, 270 causing an "encapsulation" of BC within a thick coating of secondary aerosols (Zhang et al., 2018).

Influence of wind direction
To analyse the contribution of wind direction to predicted particle concentrations, SHAP values of normalised 3-day mean u and v wind components were added up and transformed to units of degrees (Fig. 7). Generally, wind direction has a positive contribution to the model outcome when winds from the northern to northeastern sectors prevail, while negative contributions 275 are evident for southwesterly directions. Given the location of the measurement site, this pattern undoubtedly reflects the advection of particles emitted from continental Europe and/or the Paris metropolitan area under high pressure system conditions versus cleaner marine air masses during southwesterly flow (Bressi et al., 2013;Petetin et al., 2014;Petit et al., 2015;Srivastava et al., 2018). Increased concentrations of organic matter are predicted for northerly, northeasterly and easterly winds. These patterns suggest a significant contribution of advected organic particles from a specific wind sector. This is in agreement with 280 the findings of Petetin et al. (2014) who estimated that 69 % of the PM25 organic matter fraction is advected by northeasterly winds, which is related to advected particles from wood burning sources in the Paris region and SOA formation along the transport trajectories. While Petit et al. (2015) did not find a wind direction dependence of organic matter measured at SIRTA using wind regression, they reported the regional background of organic matter to be of importance. Comparing upwind rural stations to urban sites, Bressi et al. (2013) concluded organic matter is largely driven by mid-to long-range transport. Influences Germany), which are advected towards SIRTA from northeastern directions (Petetin et al., 2014;Petit et al., 2014). It is to be 290 expected that the influence of mid-to long-range transport on the particle observations at SIRTA is rather substantial, with most high pollution days affected by particle advection from continental Europe (Bressi et al., 2013). Concerning BCff and BCwb, model results suggest a dependence on wind direction during northwestern to northeastern inflow. Although BC concentrations are expected to be largely determined by local emissions (Bressi et al., 2013), e.g. from local residential areas, a substantial contribution of imported particles from wood burning and traffic emissions from the Paris region (Laborde et al., 2013;Petetin 295 et al., 2014) and continental sources is likely (Petetin et al., 2014).

Influence of feature interactions
Pairwise interaction effects, where the effect of a specific predictor on the total PM 1 prediction is dependent on the state of a second predictor, are analysed in the model. Strong pairwise interactive effects are found between MLH vs. time since last precipitation and MLH vs. maximum wind speed and shown in Figs 8a and 8b. SHAP interaction effects between MLH 300 and time since last precipitation are most pronounced for MLHs below ∼ 500 m agl (Fig. 8a). Interaction values are negative for low MLHs paired with time since last precipitation close to zero hours. With increasing time since last precipitation, interaction effects become positive, thus increasing the contribution of Tprec and MLH to the model outcome. An explanation of this pattern concerning underlying processes could be that due to the lack of precipitation, a higher number of particles is available in the atmosphere for accumulation, hence increasing the accumulation effect of a shallow MLH. In case of recent 305 precipitation, the accumulation effect of a shallow MLH is weakened. For higher MLHs, interactive effects with time since the last precipitation event are marginal. Interactive effects between MLH and wind speed are shown in Fig. 8b. Positive SHAP values for maximum wind speeds below ∼2 m/s reflect stable situations, favoring the accumulation of particles, whereas high wind speeds enhance the ventilation of particles (Sujatha et al., 2016). This can also be deduced from

Meteorological conditions of high-pollution events
To further identify conditions that favor high pollution episodes, the data set is split into situations with exceptionally high total PM 1 concentrations (>95th percentile) and situations with typical concentrations of total PM 1 (interquartile range, IQR). This   largest contributor to high pollution situations in winter is air temperature (Fig. 9). SHAP values for temperature are substantially increased during high pollution situations, when temperatures are systematically lower. Further contributing factors to high pollution situations are the lows MLHs, low wind speeds, a high average NO − 3 fraction of the previous day and negative u (i.e., winds from the east) and v (i.e., winds from the north) wind components. In winter, the PM 1 composition shows a relatively large fraction of nitrates, which is increased during high pollution situations (Fig. 9, lower panel). High concentra-330 tions of nitrate in winter can be linked to advection or to enhanced formation due to the temperature-dependent low volatility of ammonium nitrate (Petetin et al., 2014). The organic matter fraction is slightly decreased during high pollution situations. frequent in winter (Dupont et al., 2016). Positive influence of wind direction for inflow from the northern and eastern sectors are dominant during high pollution situations while inflow from the southern and western sectors prevails during average pol-335 lution situations (see Fig. 7, Bressi et al., 2013;Petetin et al., 2014;Srivastava et al., 2018). Note that the time since the last precipitation is increased during high pollution situations, but the effects on the model outcome is weak. This suggests that lacking precipitation is not a direct driver of modelled total PM 1 concentrations, but increases the contribution of other input features (see Fig. 8a) or is a meaningful factor in only some situations.
Summer total PM 1 composition (Fig. 10) is characterised by a larger fraction of organics compared to the winter season (Fig.   340 9). As a considerable fraction of organic matter is formed locally (Petetin et al., 2014), the increased proportion of organics could be due to more frequent stagnant synoptic situations that may limit the advection of transported SIA particles. In addition, the positive SHAP values of solar irradiation and temperature highlight that the solar irradiation stimulates transformation processes and increases the number of biogenic SOA particles (Guenther et al., 1993;Petetin et al., 2014). As mean temperatures are highest in summer, positive temperature SHAP values are associated with increased organic matter concentrations ( Fig.   345 5). The higher importance (i.e. higher SHAP values) of time since the last precipitation event during high pollution situations points to an accumulation of particles in the atmosphere. Dry situations can also enhance the emission of dust over dry soils (Hoffmann and Funk, 2015). The negative influences of MLH during both typical and high pollution situations reflects seasonality, as afternoon MLHs in summer are usually too high to have a substantial positive impact on total PM 1 concentrations (see Fig. 6). MLH is thus not expected to be a driver of day-to-day variations of summer total PM 1 concentrations. Note that the 350 average MLH is higher during high pollution situations, which likely points to increased formation of SO 2− 4 (see Fig. 6).

Day-to-day variability of selected pollution events
Analysing the combination of SHAP values of the various input features on a daily basis allows for direct attribution of the respective implications for modelled total PM 1 concentrations (Lundberg et al., 2020). Here, four particular pollution episodes are selected to analyse the model outcome with respect to physical processes (Figs 11-14). The examples highlight

January 2016
Prior to the onset of the high-pollution episode in January 2016 (Fig. 11), the situation is characterised by MLHs in the range 360 of 1000m, temperatures above freezing (∼5-10 • C), frequent precipitation and winds from the southwest. The organic matter fraction dominates the particle speciation. The episode itself is reproduced well by the model. According to the model results, the event is largely temperature-driven, i.e. SHAP values of temperature explain a large fraction of the total PM 1 concentration variation (note the adjusted y-axis of the temperature SHAP values). On 18th January, temperatures drop below freezing, coupled with a decrease in MLH. As a consequence, both modelled and observed PM 1 concentrations start to rise. A further 365 increase in total PM 1 concentrations is driven by a sharp transition from stronger southwestern to weaker northeastern winds  (strong negative u component, weak negative v component) on January 19th. The combined effects of these changes lead to a marked increase in total modelled PM 1 concentrations, peaking at ∼37 µg/m 3 on 20th January. On the following days, temperatures increase steadily, thus the contribution of temperature decreases. At the same time, although values of MLH remain almost constant, the contribution of MLH drops substantially from ∼5 µg/m 3 to ∼2 µg/m 3 . This is due to interactive 370 effects between MLH and the features wind speed, time since last precipitation and normalised v-wind-component. All of these features increase the contribution of MLH on 20th January, but decrease its contribution on 21st-23rd January. The physical explanation behind this pattern would be that lacking wet deposition and low wind speeds increase particle numbers in the atmosphere, while inflow from northeasterly directions increase particle numbers in the atmosphere. Given that there is now a large number of particles present, the accumulation effect of a low MLH is more efficient. The high pollution episode 375 ceases after a shift to southeastern winds and the increasing temperatures. The pollution episode is characterised by a relatively large fraction of NO − 3 and NH + 4 , which explains the strong feature contribution of temperature to the modeled total PM 1 concentration, as the abundance of these species is temperature dependent (see Fig. 5) and points to a large contribution of locally formed inorganic particles. Still, the contribution of wind direction and speed also suggests that advected secondary particles and their build-up in the boundary layer are relevant factors during the development of the high pollution episode 380 (Petetin et al., 2014;Petit et al., 2014;Srivastava et al., 2018).

December 2016
A high-pollution episode with several peaks of total PM 1 is observed in November and December 2016. The first peak on 26th Movember is followed by an abrupt minimum in total PM 1 concentrations on 28th November, and a build-up of pollution in a shallow boundary layer towards the second peak on 2nd December with total PM 1 concentrations exceeding 40 µg/m 3 . In the 385 following days, total PM 1 concentrations continuously decrease, eventually reaching a second minimum on 11th December.
A gradual increase in total PM 1 concentrations follows, resulting in a third (double-)peak total PM 1 concentration on 17th December. Total PM 1 concentrations drop to lower levels afterwards. Throughout the 3.5 week-long episode, high pollution is largely driven by shallow MLH (<∼500m), and weak north-northeasterly winds, i.e. a regime of low ventilation associated with high pressure conditions favorable for emission accumulation and possibly some advection of polluted air from the 390 Paris region. During the brief periods with lower total PM 1 concentrations, these conditions are disrupted by a higher MLH (∼28th November), or a change in prevailing winds (∼11th December). In contrast to the pollution episode in January 2016, this December 2016 episode is not driven by temperature changes. Temperatures range between ∼5-12 • C and have a minor contribution to predicted total PM 1 concentrations (see also Fig. 5), emphasizing the different processes causing air pollution in the Paris region. Note that the model is not able to fully reproduce the pollution peak on December 2nd, which may be 395 indicative of missing input features in the model. Judging from the PM 1 species composition during this time (relatively high fraction of NO − 3 and BC), it seems likely that missing information on particle emissions may be the reason for the difference between modeled and observed total PM 1 concentration.

June 2017
A period of above average total PM 1 concentrations occurred in June 2017. The episode is very well reproduced by the model, 400 suggesting a strong dependence of the observed total PM 1 concentration to meteorological drivers. Although absolute total PM 1 concentrations are substantially lower than during the previously described winter pollution episodes, the event is still above average for summer pollution levels. Organic matter particles dominate the PM 1 fraction throughout the episode, with a relatively high SO 2− 4 fraction. Conditions during this episode are characterised by strong solar irradiation (positive SHAP values) and high MLHs (mostly negative SHAP values), which show low day-to-day variability and reflect characteristic 405 summer conditions. A lack of precipitation (no rain for a period of more than 2 weeks) and high temperatures also contribute to the total PM 1 concentrations during this episode. While solar irradiation and time since last precipitation are associated  to enhanced secondary particle formation (Megaritis et al., 2014;Jiang et al., 2019). As suggested by response patterns of species to changes in MLH shown in Fig. 7, this effect is linked to an increase in SO 2− 4 concentrations. The main fraction of

March 2015
High particle concentrations are measured in early March 2015 with high day-to-day variability. This modelled course of the pollution episode is chosen to compare results to previous studies focusing on the evolution of this episode (Petit et al., 2017;420 Srivastava et al., 2018). The episode is characterised by high fractions of SIA particles, in particular SO 2− 4 , NH + 4 and NO − 3 (Fig.  14, upper panel) and similar concentrations observed at multiple measurement sites in France (Petit et al., 2017). Contributions of local sources are low and much of the episode is characterised by winds blowing in from the northwest, advecting aged SIA particles (Petit et al., 2017;Srivastava et al., 2018) and organic particles of secondary origin (Srivastava et al., 2019) towards SIRTA. A widespread scarcity of rain probably enhanced the large-scale formation of secondary pollution across western 425 Europe (in particular western Germany, The Netherlands, Luxemburg, Petit et al., 2017), which were then transported towards SIRTA. This is reflected by the SHAP values of the u and v wind components, which are positive throughout the episode (see Fig. 14g & 14h). Concentration peaks of total PM 1 are measured on 18th and 20th March. Both peaks are characterised by a rapid development of total PM 1 concentrations. As described in Petit et al. (2017), these strong daily variations of total PM 1 , which are mainly driven by the SIA fraction, could be due to varying synoptic cycles, especially the passage of cold fronts.

430
The influence of MLH and temperature is relatively small, which is consistent with the high influence of advection on total PM 1 concentrations during the episode. The exceptional character of the episode (see Petit et al., 2017) partly explains the bad performance of the model in capturing total PM 1 variability during the event. Unusual rain shortage is observed in large areas of Western Europe prior to the episode (Petit et al., 2017). While time since precipitation at the SIRTA-site is a large positive contributor to the model outcome (see Fig. 14d), it is not driving the day-to-day variations. The unusual nature of this event and 435 lacking information on emission in the source regions and formation processes along air mass trajectories in the model likely explain why the model has difficulties in reproducing this pollution episode. While this has implications for the application of explainable machine learning models for rare events, this is not expected to be an issue for the other cases and seasonal results presented here.

440
In this study, dominant patterns of meteorological drivers of PM 1 species and total PM 1 concentrations are identified and analysed using a novel, data-driven approach. A machine learning model is set up to analyse measured speciated and total PM 1 concentrations based on meteorological measurements from the SIRTA supersite, southwest of Paris. The machine learning model is able to reproduce daily variability of particle concentrations well, and is used to analyse and quantify the atmospheric processes causing high-pollution episodes during different seasons using a SHAP-value framework. As interactions 445 between the meteorological variables are accounted for, the model enables the separation, quantification and comparison of their respective impacts the individual events. It is shown that ambient meteorology can substantially exacerbate air pollution.
Results of this study point to the distinguished role of shallow MLHs, low temperatures and low wind speeds during peak PM 1 episodes in winter. These conditions are often amplified by northeastern wind inflow under high-pressure synoptic circuation. A detailed analysis reveals how the meteorological drivers of winter high-pollution episodes interact. For an episode in January 2016, model results show a strong influence of temperature to the elevated PM 1 concentrations during this episode (up to 11 µg/m 3 are attributed to temperature), suggesting enhanced local, temperature-dependent particle formation. During a different, prolonged pollution episode in December 2016, temperature levels were relatively stable and had no influence. Here, MLH (<500 m asl) was quantified to be the main driver of modelled PM 1 peak concentrations with contributions up to 6 µg/m 3 , along with wind direction contributions of up to ∼6 µg/m 3 . Total PM 1 concentrations in spring can be as high as 50 µg/m 3 . processes along the air mass trajectories, in particular of nitrate. Summer PM 1 concentrations are lower than in other seasons.
Model results suggest that summer peak concentrations are largely driven by high temperatures, particle advection from Paris and continental Europe with low wind speeds and prolonged periods without precipitation. For an example episode in June 2017, temperatures above 30 • C contribute ∼3 µg/m 3 to the total PM 1 concentration. On site scarcity of rain increases air pol-460 lution, but does not appear to be a main driver of strong day-to-day variations in particle concentrations. Presumably, this is because droughts are synoptic and are spread over several days or even weeks. Thus, they present very low inter-daily variability on the local scale. Nonetheless, Petit et al. (2017) have highlighted the link between extreme PM concentrations (especially during spring) and extreme precipitation deficit (compared to average conditions). Main drivers of day-to-day variability of predicted PM 1 concentrations are changes in wind direction, air temperature and MLH. These changes often superimpose the 465 influence of time without precipitation. Individual PM 1 species are shown to respond differently to changes in temperature.
While SO 2− 4 and organic matter concentrations are increased during both high and low temperature situations, NH + 4 and NO − 3 are substantially increased only at low temperatures. Model results indicate that SIA particle formation is enhanced during shallow MLH conditions.
The presented findings refer to the SIRTA supersite but the results are nevertheless transferable to other regions as well. For 470 example, the importance of temperature-induced particle formation processes have been shown for the U.S.A. (Dawson et al., 2007), Europe (Megaritis et al., 2014), and China (Wang et al., 2016). Hence, it is likely that the detailed, species-dependent disclosure of the nonlinear relationship between temperature and PM1 of this study holds for other urban and suburban areas.
This has implications for the PM concentrations in the context of climate change. The empirical perspective of the current study complement to the findings of various modelling studies (Dawson et al., 2007;Megaritis et al., 2013Megaritis et al., , 2014Sá et al., 475 2016; Doherty et al., 2017), the insights provided here from an empirical perspective could increase the confidence in air quality estimations under climate change. Furthermore, the impact of shallow MLHs on PM1 concentrations investigated here is comparable to results found in a previous, regional-scale study over central Europe that highlighted the dominant role of MLH on PM10 concentrations (Stirnberg et al., 2020). The importance of wind direction highlights the role of advected pollution by remote, highly polluted urban or industrial hotspots.. In general, the interpretation of pollution advection patterns requires 480 knowledge on source regions and terrain. Here, the Paris agglomeration is a major source of pollutants while the relatively flat terrain allows unimpeded advection of air masses. Urban areas in a more complex terrain would likely be affected by slightly different and possibly more complex mechanisms., such as terrain-and meteorology-dependent air stagnation events Wang et al. (2018) as well as orography driven wind and precipitation patterns (Rosenfeld et al., 2007). Still, given the task of disentangling the impact of the various meteorological drivers on air quality is already a complex scientific subject, a continental, 485 flat terrain city such as Paris was chosen as the subject area precisely to exclude other factors (such as orographic flow, or sea breeze) that would add further complexity. Certainly, the methods developed here could be transferred to more urban areas in more complex settings in the framework of future studies.
Furthermore, the analysis of meteorological drivers could be extended in future studies, e.g. by including information on anthropogenic emissions or further stations down-and upwind of SIRTA, which would allow further analysis of dominant 490 advection patterns. Furthermore, information on emissions or meteorology in the source region of air masses e.g., using Figure A1. Scatterplot for MLH [m agl] measured at SIRTA vs. MLH measured at Charles de Gaulle airport .
satellite-based observations, might be helpful to better reproduce particle transport patterns. This could be complemented by incorporating synoptic variables, e.g., the North Atlantic Oscillation (NAO) index.
For policy makers, the presented approach could prove beneficial in multiple ways. Knowledge of meteorological conditions that exacerbate air pollution could be used to issue preventative warnings to the public if these conditions are forecasted.

495
Another potential future application could be the quantitative assessment of policy measures, e.g., traffic bans, by comparing an "expected" level of air pollution under given meteorological conditions to actual observations (e.g., Cermak and Knutti, 2009), Finally, the presented model framework could be combined with short-term weather forecasts, which would allow to provide an air quality forecast based on the predictions of the statistical models.

500
Appendix A: Comparison of mixed layer height (MLH) measured at SIRTA and Charles de Gaulle airport As mentioned in section 2.2, ca. 90 missing MLH values in 2016 were replaced with measurements conducted at the Charles de Gaulle airport (see Fig. 1). Figs A2 and A1 summarize MLH values for 2016 when measurements from both sites are available (afternoon period). As shown in Fig. A1, measurements at both sites generally agree well, except for some outliers. Spearman's rank coefficient is significant (p-value < 0.05) and has a value of 0.51.