Regional-scale modelling for the assessment of atmospheric particulate matter concentrations at rural background locations in Europe

The application of regional-scale air quality models is an important tool in air quality assessment and management. For this reason, the understanding of model abilities and performances is mandatory. The main objective of this research was to investigate the spatial and temporal variability of background particulate matter (PM) concentrations, to evaluate the regional air quality modelling performance in simulating PM concentrations during statically stable conditions and to investigate processes that contribute to regionally increased PM concentrations with a focus on eastern and central Europe. The temporal and spatial variability of observed PM was analysed at 310 rural background stations in Europe during 2011. Two different regional air quality modelling systems (offline coupled European Monitoring and Evaluation Programme, EMEP, and online coupled Weather Research and Forecasting with Chemistry) were applied to simulate the transport of pollutants and to further investigate the processes that contributed to increased concentrations during high-pollution episodes. Background PM measurements from rural background stations, wind speed, surface pressure and ambient temperature data from 920 meteorological stations across Europe, classified according to the elevation, were used for the evaluation of individual model performance. Among the sea-level stations (up to 200 m), the best modelling performance, in terms of meteorology and chemistry, was found for both models. The underestimated modelled PM concentrations in some cases indicated the importance of the accurate assessment of regional air pollution transport under statically stable atmospheric conditions and the necessity of further model improvements.

1 Introduction

The increased concentration of particulate matter (PM) in the ambient environment is associated with a significant impact on human health (Anderson, 2009; Heal et al., 2012; Peters et al., 2001; Pope et al., 2002; Samet et al., 2000; Samoli et al., 2005). Continuous exposure to PM is considered to be among the top 10 most significant risk factors for public health globally, including Europe (Prank et al., 2016). The elevated PM concentrations in the atmosphere have effects on the ecosystem (acidification, eutrophication) and visibility (e.g. Putaud et al., 2010). These also affect various meteorological processes such as cloud formation and radiation. Consequently, PM has been recognised as a strong climate forcer (e.g. Andreae et al., 2005; Jiang et al., 2013) that also has an influence on Earth's heat balance through the direct radiative effects and cloud processes (Prank et al., 2016). European aerosol phenomenology studies (Van Dingenen et al., 2004; Putaud et al., 2004, 2010) have shown that the annual background concentrations of PM with an aerodynamic diameter ≤ 2.5 µm (PM_2.5) and ≤ 10 µm (PM₁₀) for continental Europe are strongly affected by regional aerosol transport. For example, long-range transport has been attributed to contributing up to about three-quarters of the total urban PM_2.5 concentrations in Finland (Karppinen et al., 2004; Pakkanen et al., 2001). A large fraction of the urban population is exposed to levels of PM₁₀ in excess of the limit values set for the protection of human health by national and international bodies. There have been numerous recent policy initiatives that aim to control PM concentrations to protect human health (EEA, 2015); yet high levels are reported regularly in different parts of the world (Kumar et al., 2015, 2016). The main problem in the assessment of PM₁₀ is in its diverse chemical composition across Europe. Nitrate is a main contributor in northwest (NW) Europe, mineral dust in south (S) Europe, desert dust from Africa over the Mediterranean, carbon in central Europe and sea salt in coastal areas of Europe. The total residence time of PM in the atmosphere is highly dependent on precipitation, which influences the deposition processes. Conversely, wind speed plays an important role in both PM advection and the alteration of PM size and composition. PM₁₀ usually deposits at closer distances from its sources than smaller particles (e.g. Dimitriou and Kassomenos, 2014). On average, the residence time of fine particles (PM_2.5) is usually about 4–6 d as opposed to 1–2 d for coarser particles (PM_2.5−10). The typical distances for deposition from the sources are around 2000 to 3000 km for the fine particles and 500 to 1000 km for coarse particles (WHO, 2006). PM₁₀ can be emitted directly to the atmosphere from various natural and anthropogenic sources (primary PM₁₀) or can be produced through photochemical reactions in the atmosphere (secondary PM₁₀). In addition, wind-blown soil and resuspended street dust contribute largely to the coarse-particle fraction (Amato et al., 2009; Forsberg et al., 2005; Harrison and Jones, 2005; Jeričević et al., 2012; Kumar and Goel, 2016; Luhana et al., 2004; Putaud et al., 2004). The contribution to PM emissions can be relevant at spatial scales ranging from local to regional including long-range transport (e.g. Juda-Rezler et al., 2011; Querol et al., 2004).

Air quality models (AQMs) play a significant role in the assessment and management of air quality. These are widely used in public health cohort studies given that the measurements are expensive and usually represent limited and small areas, e.g. rural areas or mountains (Ritter, 2013). Previous research on PM mass modelling (e.g. Vautard et al., 2007) identified the general underestimation of PM mass from large-scale models (grid spacing ∼50 km) and the difficulties in capturing the observed seasonal variations in an urban location. The complexity of PM mass modelling was also introduced in Prank et al. (2016) where various modelling systems were compared – the Unified European Monitoring and Evaluation Programme, EMEP (e.g. Simpson et al., 2012), LOTOS (e.g. Schaap et al., 2008), SILAM (e.g. Sofiev et al., 2008), and CMAQ (Community Multi-Scale Air Quality; Environmental Protection Agency) – which showed similar underestimations of PM concentrations. Applications of the Weather Research and Forecasting model with Chemistry, the WRF-Chem model (Grell et al., 2005), showed a relatively good comparison with measurements of the total PM mass over Europe (Tuccella et al., 2012), but the model did not capture the trends of PM compounds. Other studies (e.g. Saide et al., 2011) also indicated challenges in the modelling of PM mass, especially during statically stable atmospheric conditions, due to the choice of vertical and horizontal resolution as well as the influence of vertical and horizontal diffusion coefficients during model set-up (Jeričević et al., 2010). Furthermore, the WRF-Chem model was extensively tested during the intensive evaluation of online coupled models of the second phase of the Air Quality Model Evaluation International Initiative (AQMEII, 2012). During the exercise, overall underestimation of PM concentrations for all the stations was found, due to a relatively coarse grid spacing (23 km) together with the overestimation of wind speed, which can result in fast removal of pollutants from urban sources and underprediction of secondary organic aerosol (SOA) and grid-scale precipitation (e.g. Baró et al., 2015; Forkel et al., 2015). The EMEP performance is evaluated through continuous yearly technical reports such as EMEP (2016). The most recent studies showed significant technical improvements with updated initial and boundary conditions as well as with newer model versions, which include various modifications in the chemistry modules. Throughout the performed extensive tests (Gauss et al., 2016), the model generally underestimated the observed annual mean PM₁₀ levels by 24 %. However, there was an overall relatively good agreement (correlation coefficient, r=0.74) between modelled and measured annual mean PM₁₀ concentrations. A number of AQMs are currently available for practical applications. These models can be broadly divided into two main groups: offline and online models. The offline models consider solving separately meteorological conditions prior to chemistry during the simulation runs. There exists a huge variety of offline models such as the Comprehensive Air Quality Model with Extensions, CAMx (EVIRON, 2010); the Community Multi-scale Air Quality, CMAQ (U.S. Environmental Protection Agency); EMEP; and LOTOS-EUROS (e.g. Solazzo et al., 2012). In contrast to offline models, the online models were developed to include the more consistent description of processes such as atmospheric turbulence and to use a more frequent update of the meteorological variables within the chemistry part of the model. There are other reasons for online coupling such as the ability to treat feedback processes between aerosols and airflows. Examples of online models include WRF-Chem, the Environment – High-Resolution Limited Area Model (Enviro-HIRLAM), the Consortium for Small-scale Modelling Aerosols and Reactive Trace gases (COSMO-ART), and the non-hydrostatic mesoscale atmospheric model with climate module (Meso-NH-C); see, e.g., Baklanov et al. (2014).

The main objective of this research was to investigate the spatial and temporal variability in background PM concentrations using 1 year of observed data, to evaluate the regional AQM performance in simulating PM concentrations during the colder part of the year and to analyse and evaluate the episodes of regionally increased PM concentrations that occurred in November 2011 in eastern and central Europe (the Pannonian Basin) during statically stable atmospheric conditions followed by drought periods. In this particular case, the pollution problems appeared to be of considerable concern in Hungary; e.g. smog alerts were issued in Budapest and eastern Hungary, various cars with high environmental impact were banned from the roads, and a ban was also issued on procedures such as burning leaves and garden debris (https://thecontrarianhungarian.wordpress.com/2011/11/08/hungarian-news-digest-nov-7-2011/, last access: 12 May 2020). Based on the analysis in Spinoni et al. (2015), the Pannonian Basin was characterised as an area with increased drought frequency per decade during the period from 1950 to 2012. This can have a strong effect on air quality problems, e.g. a dust-bowl effect (Stahl et al., 2016). Further assessment is conducted by applying two regional models: the offline Unified EMEP and the online coupled WRF-Chem in the simulation of PM mass transport. Model results are compared against observed concentrations at rural background sea-level, elevated and mountain stations in Europe.

Throughout the analysis, the indication of problems is given in the application of both regional models in simulating PM concentrations at different elevations (sea-level, elevated and mountain stations). We provide an individual validation of widely used different set-ups of the modelling systems without harmonisation of emission and meteorological input fields. This is a different approach than in, e.g., AQMEII exercises and enables an essential scientific baseline for choosing the appropriate model for future needs in terms of resolution, physical parameterisation, emission dataset and the complexity of orography representation in practical applications. Due to the complexity of air quality problems regarding PM, this work aims at filling the gaps in knowledge of regional modelling of PM over eastern Europe in terms of less information about PM concentrations (EEA, 2013) and therefore low accuracy in the PM emission inventory, and it fits in with the problems addressed in most of the air quality plans in Europe (Miranda et al., 2015).

2 Methodology

2.1 Measurements

The measurements of PM₁₀ from the rural background stations were taken from two available air quality databases. These were AirBase, the European air quality database maintained by the European Environmental Agency (https://www.eea.europa.eu/data-and-maps/data/airbase-the-european-air-quality-database-7, last access: 12 May 2020), and the database developed under the EU-funded PHARE 2006 project Establishment of Air Quality Monitoring and Management System, where 12 new rural stations were established in Croatia for PM measurements in 2011. For this study, PM₁₀ concentrations were available from six rural background stations in Croatia. The monitoring stations were divided into three categories based on their elevation to examine the spatial variability of pollution and to test the model performance at different levels: (i) sea-level (altitude from 0 to 200 m), (ii) elevated (from 200 to 500 m) and (iii) mountain stations (> 500 m). The differentiation of stations with respect to their elevation is important when dealing with station representativeness in models. According to current knowledge, it is found that numerical models perform differently at higher altitudes. This is mostly related to the vertical resolution of the model within the boundary layer (Bernier and Bélair, 2012). With respect to the elevation, the total numbers of stations used for further analysis (Sect. 3.1) and model validation (Sect. 3.3) are shown in Table 1. When interpreting average observed yearly, seasonal and episode PM₁₀ concentrations, it is important to note that the majority of the surface stations are in northern and western Europe, while the elevated and mountain stations are generally situated in central and eastern Europe. The density of rural background stations varies geographically with a significantly greater number of stations in western and northern Europe compared to central and eastern Europe. The PM₁₀ measurements were acquired with different approaches: the gravimetric method (EN12341) using high-volume samplers (HVSs) and low-volume samplers (LVSs), β-attenuation monitoring (e.g. Willeke and Baron, 1993), TEOMs (tapered element oscillating microbalances) measurements (e.g. Patashnick and Rupprecht, 1980), and the optical particle counters of the GRIMM 180 instrument. The comparison of the PM₁₀ concentration data obtained by different measurement methods is still considered to be a complicated issue. The standard gravimetric method (EN12341) is a classic method of weighing the mass deposited on a filter. It is accepted as a standard reference method against which all other measurement methods are validated (Noble et al., 2001; EC, 2010). Although this is a standard method used for compliance reasons in the EU, there are numerous studies showing that chemical reactions between the air and the deposited particles, as well as within the aerosol mass, also compromise these measurements. The ambient temperature and relative humidity greatly influence the actual mass loaded on the filter (Allen et al., 1997; Eisner and Wiener, 2002; Pang et al., 2002). For example, aerosol particles can contain up to 30 % water at 50 % relative humidity (Putaud et al., 2004). Conversely, calibration, temperature and humidity issues can create artefacts that must be taken into account for TEOMs and β-attenuation monitoring (Allen et al., 1997; Hauck et al., 2004). Lacey and Faulkner (2015) addressed three objectives for the treatment of uncertainties gained with PM measurements: estimate the uncertainty, identify the measurement with the greatest impact on uncertainty and finally determine the sensitivity of total uncertainty to all measured parameters. As is common in these types of studies, the authors did not consider the uncertainty of measurements in further analysis.

Table 1 The number of stations used in the analysis.

2.2 Statistical analysis

The evaluation of model performance is a comprehensive task. In order to evaluate and validate modelling performance, various statistical measures such as bias (BIAS), index of agreement (IOA), correlation coefficient (r), root mean square error (RMSE), normalised mean square error (NMSE), and systematic (NMSE_sys) and unsystematic (NMSE_unsys) normalised mean square error were used (Chang and Hanna, 2004):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{BIAS}=\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{M}-\stackrel{\mathrm{‾}}{O}}{\stackrel{\mathrm{‾}}{O}}}\right)\times \mathrm{100}\mathit{\%},\text{(2)}& {\displaystyle }\mathrm{BIAS}=\stackrel{\mathrm{‾}}{M}-\stackrel{\mathrm{‾}}{O},\text{(3)}& {\displaystyle }\mathrm{IOA}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(O_i-M_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N\left(\mathrm{abs}\right(M_i-\stackrel{\mathrm{‾}}{O})+\mathrm{abs}(O_i-\stackrel{\mathrm{‾}}{O}){)}^{\mathrm{2}}}},\text{(4)}& {\displaystyle }r={\displaystyle \frac{n\sum O_i{M}_i-\sum O_i\sum M_i}{\sqrt{n\sum O_i^{\mathrm{2}}-(\sum O_i{)}^{\mathrm{2}}}\sqrt{n\sum M_i^{\mathrm{2}}-(\sum M_i{)}^{\mathrm{2}}}}},\text{(5)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^N(M_i-O_i{)}^{\mathrm{2}}},\text{(6)}& {\displaystyle }\mathrm{NMSE}={\displaystyle \frac{\stackrel{\mathrm{‾}}{(O-M{)}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{O}\stackrel{\mathrm{‾}}{M}}},\text{(7)}& {\displaystyle }{\mathrm{NMSE}}_{\mathrm{sys}}={\displaystyle \frac{\mathrm{4}{\left(\frac{\left(\stackrel{\mathrm{‾}}{O}-\stackrel{\mathrm{‾}}{M}\right)}{O.5\left(\stackrel{\mathrm{‾}}{O}+\stackrel{\mathrm{‾}}{M}\right)}\right)}^{\mathrm{2}}}{\mathrm{4}-{\left(\frac{\left(\stackrel{\mathrm{‾}}{O}-\stackrel{\mathrm{‾}}{M}\right)}{O.5\left(\stackrel{\mathrm{‾}}{O}+\stackrel{\mathrm{‾}}{M}\right)}\right)}^{\mathrm{2}}}},\text{(8)}& {\displaystyle }\mathrm{NMSE}={\mathrm{NMSE}}_{\mathrm{sys}}+{\mathrm{NMSE}}_{\mathrm{unsys}},\end{array}$$

where M stands for model predictions and O for observations.

As there is no single best modelling performance measure, it is recommended by Chang and Hanna (2004) that a suite of different performance measures needs to be applied. Results should be carefully interpreted by taking into account advantages and disadvantages of all applied statistical measures and ensuring that those are complementary to each other and lead to the same conclusion on the certain ability of the model performance. Therefore as already previously noted in this section, a set of different statistical measures is used in order to understand the ability of the model to properly estimate high-pollution episodes of PM₁₀ concentrations and to evaluate the relations between chemical and meteorological parameters. BIAS refers to the arithmetic difference between M and O indicating the model's general overestimation or underestimation of analysed parameters. It is known that a model whose predictions are completely out of phase with observations can still have a BIAS=0 because of compensating errors. Different BIAS was used: for evaluating model performance regarding PM₁₀ we used BIAS under Eq. (1) as opposed to meteorological parameters under Eq. (2). r and IOA are dimensionless measures of model accuracy. r is sensitive to good agreement of extreme data pairs, and a scatter plot might show generally poor agreement, but the presence of good agreement for a few extreme pairs will greatly improve r. The IOA is the ratio of the mean square error and the potential error and then subtracted from 1 (Willmott, 1984). The IOA varies from 0 to 1, with higher index values indicating that M have better agreement with the O. Although the IOA provides some improvement over r, it is still sensitive to extreme values due to the square differences in the mean square error in the numerator. RMSE gives information on the spread of the residuals from the regression line; it depends greatly on the magnitude of the parameter on which RMSE is applied and therefore it cannot be compared with RMSE of some other parameter. NMSE_sys is a measure with which NMSE_unsys provides information on systematic and unsystematic (random) errors in the model.

2.3 Boundary layer height determination

One of the widely used methods for deriving boundary layer height from numerical models is based on the assumption that turbulence collapses to laminar flow when the bulk Richardson number Ri_B, exceeds values of a critical Ri_B (∼0.25 and larger), and the height at which this occurs can be considered as a boundary layer height (Jeričević et al., 2010). Using sounding and modelled data, Ri_B was calculated based on the following expression and shown in Sect. 3.3.3:

$$\begin{array}{}\text{(9)}& {\mathit{Ri}}_{\mathrm{B}}={\displaystyle \frac{g(z-z_{\mathrm{0}})}{\stackrel{\mathrm{‾}}{\mathit{\theta }\left(z\right)}}}{\displaystyle \frac{\mathit{\theta }\left(z\right)-\mathit{\theta }\left(z_{\mathrm{0}}\right)}{\left(u\right(z){)}^{\mathrm{2}}+(v\left(z\right){)}^{\mathrm{2}}}},\end{array}$$

where z is the height of the particular model level, z₀ is the height of the first level in the model, θ(z) is the potential temperature at the height z, θ(z₀) is the potential temperature at the height z₀, and $\stackrel{\mathrm{‾}}{\mathit{\theta }\left(z\right)}$ is the averaged potential temperature between the first level (z₀) and particular level (z). u(z) and v(z) are the wind components on particular levels.

Comparison of estimated planetary boundary layer height (PBLH) was carried out using Eq. (9) rather than comparing the direct output of model-derived PBLH values as each model uses a different method for the calculation of the PBLH. By using the same methodology for PBLH determination, uncertainties are reduced and the more realistic evaluation of two modelled PBLH values is ensured.

2.4 Air quality models

This work is based on the intensive tests performed in Gašparac et al. (2016), where the WRF-Chem, Unified EMEP and WRF-CAMx models were evaluated against the surface measurement stations over Croatia under different atmospheric static stability conditions. Here, both EMEP and WRF-Chem AQMs are used to determine the spatial and temporal distribution of PM₁₀ concentrations and possible transboundary transport and to evaluate the performance of the individual model systems during a 1-month period at the sea-level, elevated and mountain rural background stations in Europe. The WRF-Chem simulation is performed from 29 October to 30 November and EMEP from 1 October to 30 November. As all statistical analysis was done for dates after 1 November the simulation length was long enough to overcome the effects of spin-up time.

2.4.1 The EMEP model

The EMEP chemical transport model (Simpson et al., 2012), developed by the Meteorological Synthesizing Centre-West (MSC-W) was used to perform calculations of PM₁₀ concentrations (http://www.emep.int/, last access: 12 May 2020). The model domain encompassed all of Europe with a horizontal grid spacing of 50×50 km², extending vertically from surface level (first model level height around 42 m) to the tropopause at 100 hPa, as seen in the Supplement Fig. S1. The basic physical formulation of the EMEP model is derived from Berge and Jakobsen (1998). The model derives its horizontal and vertical grid from the input meteorological data. The daily meteorological input data used for the EMEP/MSC-W model for 2011 were based on experimental forecast runs with the Integrated Forecast System (IFS), a global operational forecasting model from the European Centre for Medium-Range Weather Forecasts (ECMWF). Vertically, the 60 eta levels of the IFS model were interpolated onto the 37 EMEP sigma levels. The emission input for the EMEP/MSC-W model, with a horizontal grid spacing of 50×50 km², consists of gridded annual national emissions based on emission data reported every year to EMEP/MSC-W (until 2005) and to the Centre on Emission Inventories and Projections (from 2006) by each participating country. The standard emissions input required by the EMEP model consists of gridded annual national emissions of sulfur dioxide (SO₂), nitrogen oxides (${\mathrm{NO}}_x=\mathrm{NO}+{\mathrm{NO}}_{\mathrm{2}}$), ammonia (NH₃), non-methane volatile organic compounds (NMVOCs), carbon monoxide (CO) and particulates (PM_2.5, and PM_2.5−10). The PM categories can be further divided into elemental carbon, organic matter and other compounds as required. Emissions can be set from anthropogenic sources such as the burning of fossil and biomass-based fuels and solvent release or from natural sources such as foliar volatile organic compound (VOC) emissions or volcanoes. Several sources are challenging to categorise into anthropogenic versus natural categories (Winiwarter and Simpson, 1999), for example, emissions of NO from microbes in soils being promoted by N-deposition and fertiliser usage. The anthropogenic emissions are categorised into 11 SNAP (Selected Nomenclature for sources of Air Pollution) sectors based on their sources. Emission integration during simulation is distributed vertically, based on the SNAP sectors and plume-rise calculations performed for different types of emission sources and, temporally, based upon time factors (i.e. monthly, daily, day of week, weekly, hourly).

Regarding the planetary boundary layer parameterisations under statically stable atmospheric conditions, EMEP includes a non-local vertical diffusion scheme based on a linear exponential profile with coefficients calculated from large eddy simulation (LES) data and boundary layer height determined using the bulk Richardson number method (Jeričević and Večenaj, 2009; Jeričević et al., 2010; Simpson et al., 2012). Other mechanism used in this work (e.g. chemical scheme: EmChem09; chemical preprocessor: GenChem) are described in Simpson et al. (2012).

The above set-up of the EMEP model with the IFS meteorology as an initial and boundary meteorological condition is later on referred to and used as the “EMEP model”. Any further comparison of meteorological conditions obtained in EMEP simulations is related to the IFS model and PM₁₀ to the choice of EMEP chemistry parameterisation.

2.4.2 The WRF-Chem model

The WRF-Chem model is the WRF (Weather Research and Forecasting) model (http://www.wrf-model.org, last access: 12 May 2020) coupled with chemistry. It is a state-of-the-art air quality model (Grell et al., 2005), in which the chemistry (emission, transport, mixing, and chemical transformation of trace gases and aerosols) is simultaneously simulated with meteorology (online coupling). The WRF is a mesoscale numerical weather prediction system designed for operational forecasting needs and atmospheric research (Skamarock et al., 2008). The model set-up was based on earlier research (Gašparac et al., 2016; Grgurić et al., 2013; Jeričević et al., 2017) where the results were evaluated against measurements at meteorological stations in Croatia. In this paper, we used the WRF-Chem version 3.5.1. A Mercator projection was used in a one-domain run on 170 points in the east–west direction and 145 points in the north–south direction, with a cell size of 18×18 km² (Fig. S1) and a vertical grid spacing encompassing the atmosphere from surface level (first model level height around 22 m) to the height of ∼23 km in 50 unequally sorted sigma levels that were more densely distributed near the ground level. Initial and boundary meteorological conditions were provided by NCEP (National Centers for Environmental Prediction) Final Analysis (FNL ds083.2) with 1^∘ of horizontal resolution and a time step of 6 h. They were selected based on previous research and other conducted studies with the WRF or WRF-Chem model (Gašparac et al., 2016; Grgurić et al., 2013; Jeričević et al., 2017; Syrakov et al., 2015). FNL analyses are a product of the Global Data Assimilation System (GDAS), which continuously makes multiple analyses of collected observational data from the Global Telecommunications System (GTS) and various other sources. The whole analysis is available at 26 pressure levels from the surface to a height of ∼28 km. The input emissions were prepared via the PREP-CHEM Sources tool (Freitas et al., 2011) with the EDGAR (version 4.3.1., Emissions Database for Global Atmospheric Research) inventory for the year 2011. Biogenic emissions were calculated from MEGAN (Model of Emissions of Gases and Aerosols from Nature; Guenther et al., 2006) and lateral boundary and initial conditions were created from the global chemistry model MOZART (Emmons et al., 2010). The detailed WRF-Chem set-up is shown in Table 2.

Table 2 Details of the WRF-Chem parameterisations.

It is worth pointing out that the results of statistical analysis and model evaluation further on in the text will not describe the performance of the model itself but rather will describe the performance of a set of selected parameterisations and chemical and meteorological initial and boundary conditions used in the WRF-Chem model. Following this, when referring to the “WRF-Chem model” in the text, the authors are referring to the WRF-Chem model with the above-described set-up (Table 2).

3 Results

Available daily averaged rural background PM₁₀ concentrations ($\left(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}\right)$ over Europe (Table 1) were analysed in the following sections with annual temporal variations and the episodes of very high $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations that occurred during November 2011.

3.1 Analysis of PM measurements

We analysed the $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ measurements from 310 stations over a period of 1 year during 2011. Following the air quality report in Europe (EEA, 2013), $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ limit values (2008/50/EC Directive, LV =50 µg m⁻³) were exceeded at both urban and rural sites in Europe during 2011. These “hotspots”, locations with exceedances of the LV, were in south Poland, the Czech Republic, the Po Valley, the Balkan Peninsula, Portugal and Turkey. In this work, we focused on the area and rural background stations shown in Fig. 1. The analysis of measurements from 310 rural background stations showed that observed $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ exceeded LV 5456 times during 2011 and were mainly located in the hotspot areas (Fig. 1). The seasonal variation in $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ during 2011 was significant at the 5 % level (based on analysis of variance, ANOVA; p=0). The applied ANOVA is calculated via scipy python package. This particular one-way ANOVA tests the null hypothesis that two or more groups have the same population mean. The p value is common variable used in hypothesis testing, the smaller the p value, the stronger is the evidence that the hypothesis needs to be rejected (Heiman, 2001).

Figure 1 Number of days exceeding the daily limit PM₁₀ value (LV) at rural background stations during the year 2011 in the domain of the research. Stations marked with a grey circle represent less than or equal to 35 permitted exceedances during the year (2008/50/EC Directive).

Spatially averaged seasonal values of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ were 21.62, 21.74, 14.96 and 20.87 µg m⁻³ for DJF, MAM, JJA and SON, respectively. Only during summer (JJA) was a decrease found with respect to other seasons over Europe. However, it should be noted that significant differences in PM levels across Europe are recognised (Putaud et al., 2004) and a deeper analysis of spatial and temporal variations in background PM₁₀ concentration is needed. Fig. 2 presents individual $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values for each rural background station (panel b), spatially averaged $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ over the all stations (green line, panel a) and the maximum $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values among all rural background stations (red line, panel a) during 2011. The time series of these $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations indicate the increase in concentrations at all rural background stations (Fig. 2) during DJF and SON seasons (i.e. the colder part of the year). During these seasons, $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values at all rural background stations were relatively high, reaching 40 µg m⁻³. During the colder part of the year, most of the stations recorded $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values above the permitted LV which is mainly due to increased emissions from domestic heating and industrial activities (EEA, 2013). Moreover, according to, e.g., EEA (2013) and Saarikoski et al. (2008), aside from the primary sources (natural and anthropogenic), the secondary inorganic aerosols (SIAs) and SOAs vary substantially across Europe from season to season, which indicates the presence of various PM₁₀ sources. SIA contributions are mostly related to SON–DJF, domestic heating and large combustion plants, while SOA contribution is instead related to MAM–JJA seasons, e.g. emissions from vegetation. This can explain the relatively high daily concentrations in the MAM season (Fig. 2).

Figure 2 The spatial average (a) over all the rural background stations (the green line, corresponding to the right green y axis) and the maximum of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ for all rural background stations (the red line, corresponding to the left red y axis) and $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ (b) during 2011. The values above 50 µg m⁻³ (red colour) represent values above the daily limit values for PM₁₀ under the 2008/50/EC Directive.

3.2 Analysis of PM measurements and meteorological conditions during episodes in November 2011

Further analysis of the observed and modelled $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values is focused on November 2011, as the highest $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations were present during the colder part of the year and prevailing meteorological conditions enabled the accumulation of the pollutants in the lower layers of the atmosphere over Europe. According to Blunden et al. (2012), a strong high-pressure field encompasses the area over central and southern Europe during November 2011. Moreover, this month was the coldest in 2011 and extremely dry; it was the driest month in Bulgaria and Serbia with less than 25 % of the national total averaged precipitation. During the SON season in 2011, anticyclonic conditions prevailed and below-average precipitation conditions were recorded. Following Cindrić et al. (2016), the drought was present in the continental part of Croatia, encompassing the Pannonian Basin and surrounding countries, and was characterised by an extremely long duration. It started in February 2011 and reached the most intense extremely dry conditions in November, when an increase in $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ was recorded at the majority of the analysed rural background stations (Fig. 2). In western Europe, the autumn season temperature was above average normal (1961–1990) and was characterised by a prevailing high-pressure field. This was observed particularly in November, during which monthly average temperature records were exceeded (e.g. UK, France and Switzerland reported their second warmest autumn in last 100 years). Contrary to the western Europe, the increased nocturnal cooling decreased temperatures in southeastern Europe. The dominating high-pressure field resulted in a decrease in precipitation in some western and central Europe countries, e.g. southern France, the Alpine region, Germany, Austria, the Czech Republic, Slovakia and Hungary. All those countries reported the driest November in more than the last 100 years (Blunden et al., 2012). In order to identify the episodes and the areas of enhanced $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values, differences (DFs) between the $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ and annually averaged PM₁₀ ($\left(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{a}}\right)$ at rural background stations were used, defined as

$$\begin{array}{}\text{(10)}& \mathrm{DF}={\displaystyle \frac{(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}-(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{a}}}{(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{a}}}}\times \mathrm{100}\mathit{\%}.\end{array}$$

The spatial distribution of DF values in percentage is shown in Fig. S2. The significant increase in $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ is defined as an increase in DF of more than 100 % with respect to the annual mean. If a significant increase in DF was detected and lasted at least 2 consequent days, the area was identified as an area experiencing a high-pollution episode. During November 2011, a significant increase in $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ occurred generally over the addressed hotspots within the domain, and two high-pollution episodes (DF > 100 %) were found (Figs. 3–4). During both episodes identified, the highest peaks (9 November in the first episode, Fig.  3; 14 November in the second episode, Fig. 4) occurred in the area of central Europe and coastal part of western Europe with a DF above 200 %. Later on, observed meteorological conditions (daily averaged sea-level pressure field ($(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$), daily averaged surface temperature ($(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$), daily averaged relative humidity ($(\stackrel{\mathrm{‾}}{\mathrm{rh}}{)}_{\mathrm{d}}$), Figs. 3–4; daily averaged surface wind speed ($(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$) and direction ($(\stackrel{\mathrm{‾}}{\mathrm{wd}}{)}_{\mathrm{d}}$), Fig. 5) along with the DF (Figs. 3–4) were analysed to determine the mechanisms and relationships between the meteorology and the high-pollution episodes.

Figure 3 DF as “Conc ∼ year” and measurements from synoptic stations (relative humidity (Rel hum), ambient temperature at 2 m (Temp) and surface pressure) from Ogimet (https://www.ogimet.com/, last access: 12 May 2020) during the first large-scale episode (5 to 9 November). Stations with a temperature between 0 and 5 ^∘C are marked with little grey dots due to better representativeness on the map.

Figure 4 Same as Fig. 3, but during the second large-scale episode (12 to 16 November).

At the beginning of November, values of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ were mainly at (or lower than) the mean monthly average values over most of the analysed stations while an increase in DF ranging from 50 % to 100 % above annual averages was found over hotspots areas (south Poland, the Czech Republic, Po Valley, Balkan Peninsula; Fig. S2). On 3 November, cyclone Rolf in the Gulf of Genoa generated intense rainfall in northern Italy (not shown). These conditions were followed by high S to SE winds over the Adriatic Sea and nearby countries in the following days (Blunden et al., 2012). The characteristic meteorological conditions during or following Genoa low cyclones are a strong flow aloft (Sirocco wind over the Adriatic Sea and Italy), rainfall in mid-central Europe (Austria, the Czech Republic and Poland) and the formation of high-$(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ fields over eastern Europe (Blunden et al., 2012). From 5 November, a first large-scale episode (DF > 100 %, Fig. S2) started in central and northern Europe. The onset of the event was in Poland and northeastern Germany and encompassed the coastal areas of northern Europe, the Benelux countries and northern France in the following days until 9 November. During the first episode, a high-$(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ field (Fig. 3) formed over continental Europe, first affecting the east of Europe and gradually spreading to western Europe. Over the affected area (DF > 100 %, Fig. 3), the wind speed was generally reduced below 3 m s⁻¹, except at some isolated stations (Fig. 5, left). Moderate to strong NE wind (5–6 m s⁻¹) started to blow in coastal and northern Europe from 7 November until the end of the first episode, when they turned into an ESE direction (Fig. 5, left). Over the mountainous region in central Europe (the Czech Republic, Slovakia and south Germany), the wind speed was persistent during the episode, with relatively high magnitude (above 7 m s⁻¹) and generally in a SSE direction. Over the area with increased concentrations (DF > 100 %, Fig. 3), a gradual moderate decrease in $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ from east to west from the beginning to the end of the first episode was found (i.e. Poland < 0 ^∘C, Germany, the Czech Republic and Slovakia 0–5 ^∘C). On 10 November the wind speed was lower than 3 m s⁻¹ over all of Europe (not shown); values were reduced and comparable to $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{a}}$ (Fig. S2).

Figure 5 Daily averaged wind speed and directions during two high-pollution episodes: from 5 to 9 November (episode 1, left) and from 12 to 16 November (episode 2, right). Stations with a measured wind speed below 3 m s⁻¹ are marked with grey dots. Source of measurements: Ogimet (https://www.ogimet.com/, last access: 12 May 2020).

A building up of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ started again from 12 November (Figs. S2 and 4), mainly affecting stations in central and coastal western Europe. The observed concentrations exceeded the annual averages by up to 100 % (DF) affecting the areas with hotspots (southern Poland, the Czech Republic, Benelux countries) and up to 200 % in central Germany and Slovakia (Fig. 4). In the following days, from 13 to 16 November, increased concentrations (DF > 100 %) encompassing the area from central Europe in the northwest direction through coastal areas in Germany, the UK and Ireland were present in the southeastern direction across the Czech Republic, Austria, Slovenia, western Hungary and Croatia. During this, second episode, a high-$(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ field again influenced the weather conditions (Fig. 4). Low $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ (<3 m s⁻¹; Fig. 5, right) and a decrease in $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ were found with the lowest $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ measured in eastern and central Europe (Fig. 4, below −5 ^∘C). Previously mentioned persistent conditions influenced the formation of statically stable atmospheric conditions during this episode (see Sect. 3.3.3). Over particular areas with highly increased concentrations (DF > 200 %, Poland, Germany, Slovakia, the Czech Republic; Fig. 4), an increase in $(\stackrel{\mathrm{‾}}{\mathrm{rh}}{)}_{\mathrm{d}}$ was found, except in the Pannonian Basin (Fig. S3) where relatively lower $(\stackrel{\mathrm{‾}}{\mathrm{rh}}{)}_{\mathrm{d}}$ and higher $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ values of up to 20 % and 5 ^∘C, respectively, were recorded in comparison with the surrounding areas. Moreover, within the areas of the Pannonian Basin, a high $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ and low wind speed conditions prevailed for 1 d longer (Figs. 4–5, right) in comparison with the surrounding areas. On 19 November a large-scale decrease in $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ was detected and values of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ were generally reduced to those of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{a}}$ at all stations (Figs. S2, S6).

3.3 Model evaluation

Numerical simulations using the EMEP (with a grid spacing of 50×50 km²) and WRF-Chem (with a grid spacing of 18×18 km²; Fig. S1) models were provided for November 2011 to evaluate the performances of the individual, state-of-the-art models during November 2011 and to further investigate the processes contributing to the increased concentrations during the high-pollution episodes. It is worth noting that differences between the emission databases used were found in the spatial variability of PM₁₀ emissions and in the gridded input emission fields above the entire domains of EMEP and WRF-Chem. Notable differences in emissions were found over the coastal areas and eastern part of the domain, particularly over Bosnia and Herzegovina, Serbia and Hungary which are crucial for the case studies analysed here. Aside from this, the difference in vertical resolution (first model level height – EMEP at 46 m, WRF-Chem at 22 m) can have a strong impact on surface concentrations and thus can be related to the differences in surface PM₁₀ concentrations obtained from the two models used.

3.3.1 Evaluation of model performances during November 2011

Meteorological conditions

Vertical wind profile plays an important role in the dispersion of particulate matter. Hence, a validation of the modelled wind speed against measurements using mast-mounted instruments (Fig. 6; Cabauw, Netherlands, 51.97^∘ N, 4.95^∘ E, and Karlsruhe, in the western part of Germany, 48.98^∘ N, 8.39^∘ E) was performed. During November there was no significant difference between modelled vertical profiles of wind speed below 75 m (Fig. 6) for both sites. Modelled vertical wind profiles were close to measurements at the Cabauw site (up to 75 m), while at Karlsruhe the models underestimated the observed wind speed values in the first 180 m for the WRF-Chem model and much higher above ground level for EMEP. The relatively coarse horizontal resolutions of the models have a great impact on wind values (e.g. Jeričević et al., 2012), which is why the modelled values correspond better to the observed wind values at the Cabauw site, situated in flat terrain, than to the values observed over the moderately complex terrain at the Karlsruhe site. Above 100 m, a change in the slope of the vertical wind speed profile for WRF-Chem was found. The difference in model performance above the surface layer was previously attributed to the proper choice of boundary layer parameterisation in Boadh et al. (2016).

Figure 6 Vertical profiles of measured and modelled wind speeds at Karlsruhe (a, measurement source: Institute of Meteorology and Climate Research, Atmospheric Environmental Research, Karlsruhe Institute of Technology) and the Cabauw mast station (b, measurement source: Cesar Observatory, http://www.cesar-observatory.nl, last access: 12 May 2020) during November 2011.

The modelled $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ and $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ were compared to measurements from 920 synoptic stations within the domain taking into account the elevation of the station. A detailed statistical evaluation of the two individual model performances was conducted by calculation and analyses of six different statistical measures (Fig. 7): BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$, IOA$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$, $r\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$, RMSE$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$, NMSE${}_{\mathrm{sys}}\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$ and NMSE${}_{\mathrm{unsys}}\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, and $\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$. In the following Fig. 7, individual scales for each analysed meteorological parameter are given as their magnitudes differs greatly. Statistic measures calculated for wind speed are given in units of metres per second, temperature in degrees Celsius and pressure in hectopascal. This is important for the interpretation of model scores in simulating different meteorological parameters as, e.g., RMSE or NMSE depend on their magnitude. Furthermore, the results from Fig. 7 should be viewed as individual model performance rather than inter-comparison of two different model performances. According to BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}})$, the WRF-Chem model generally overestimated the observed $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, which is in accordance with other similar studies (e.g. Solazzo et al., 2012). The median of overestimation of $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ increases with the station altitude; BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}})$ was 1.8 m s⁻¹ at sea-level, 1.9 m s⁻¹ at elevated and 2.8 m s⁻¹ at mountain stations. WRF-Chem successfully predicted $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ and $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ as BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}},\left(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}\right)$ values were very low at sea level and elevated stations while small to moderate (BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ )∼1.2 hPa, BIAS$\left(\right(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}})\pm \mathrm{1}$ ^∘C) at mountain stations. The BIAS ($\left(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}\right)$ increases with height for both models. On elevated stations, a median of BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}})$ decreased up to 1 hPa for both models. On mountain stations the spread of BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}})$ is higher with respect to lower altitudes (−5 to 2.6 hPa for WRF-Chem and −4.9 to 4 hPa for the EMEP model). The median of BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}})$ is the same as for elevated stations for the WRF-Chem model and 0.5 hPa for the EMEP model. The EMEP model predicted $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ and $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ well with low BIAS values at sea-level and elevated stations, while for surface $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ values, underestimation was found (BIAS$\left(\right(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}})$ $\sim -\mathrm{2}$, 3 and 4 ^∘C at sea-level, elevated and mountain stations, respectively). The median IOA$\left(\right(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}})$ was relatively high for both models, while for IOA$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}})$ it was to a small extent lower. For both parameters a decrease in performance with height was found. This indicates problems in simulations with regional models over complex terrain, which is confirmed by the values of r that were consistent for both models. As a result of small BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}})$ over sea-level and elevated stations, the IOA$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}})$ was close to 1. However, over the mountain stations a high spread of values was found as the formulation of IOA is very sensitive to the extreme values. The models did not show any substantial unsystematic and systematic errors for $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$. The range of both systematic and unsystematic errors increased with height for $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$; the median values of NMSE_sys and NMSE${}_{\mathrm{unsys}}\left(\right(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}})$ for the EMEP model were the highest for elevated stations. In the case of the WRF-Chem model, NMSE${}_{\mathrm{sys}}\left(\right(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}})$ increases with height, while for the EMEP model, the highest NMSE${}_{\mathrm{unsys}}\left(\right(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}})$ median was found at elevated stations.

Figure 7 Intercomparison of the applied statistical measures (BIAS, IOA, r, RMSE, NMSEsys, NMSEunsys) between modelled (WRF-Chem – red boxes; EMEP – blue boxes) and measured (from 920 meteorological stations across all of Europe, source: Ogimet (https://www.ogimet.com/, last access: 12 May 2020.)) wind speed (WS, //), temperature (Temp, ^∘∘) and surface pressure (Pressure, $\left|\right|$) during November 2011 for sea-level (S), elevated (E) and mountain (M) stations. The units of selected meteorological parameters are metres per second for wind speed, degrees Celsius for temperature and hectopascal for surface pressure.

Overall, during a 1-month period of simulation, EMEP had the lowest systematic errors for $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, while WRF-Chem had the lowest systematic errors for $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$. Based on the statistics given, overall model performance regarding meteorological parameters was in accordance with similar modelling studies. For example, negative BIAS and high r for $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ was found in, e.g., Skjøth et al. (2015) and Qu et al. (2014). Positive BIAS for $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ was already addressed as an issue in related studies such as, e.g., Baró et al. (2015) and Forkel et al. (2015), while results for $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ for sea-level and/or elevated stations are in accordance with, e.g., Qu et al. (2014).

Chemistry

The modelled $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values were compared with the available corresponding measurements (Table 1) with respect to height by applying statistical measures (Fig. 8, Table S1). Although the number of stations varies within altitude groups (Table 1), the overall model performance can be inferred from the above-mentioned figure and table. The underestimation of concentrations was found at sea-level (the median of −44 % and −26 % for the WRF-Chem and EMEP models, respectively) and elevated stations (−55 % and −29 % for the WRF-Chem and EMEP models, respectively; Figs. S4–S5). At mountain stations, EMEP had good agreement of ∼13 %, while underestimation with respect to WRF-Chem is still present (∼33 %). According to Figs. S4–S5, the BIAS$\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ in both model simulations showed a similar distribution with respect to the height of the station, i.e. moving from underestimation towards overestimation. IOA$\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ was generally equally persistent with height for both models (Fig. 8) with a somewhat higher score for simulations with the EMEP model except for the sea-level stations where the median of both models had equal values (0.9, Table S1). The highest $r\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ values were above 0.87 for both models; however, the overall performance in terms of r (median, Table S1) for both models was relatively low, particularly for the elevated and mountain stations. The average values over the domain for the WRF-Chem and EMEP models were 0.39, 0.21 and 0.19 and 0.48, 0.28 and 0.24 for sea-level, elevated and mountain stations, respectively. High variability in r values over the domain for both models is found (Figs. S4–S5). As r is a measure of linearity and is highly dependent on the estimation of peak values and trends, the low values at all stations are attributed to a mismatch between modelled and measured peak values during the period of analysis. Even a small discrepancy between measured and modelled $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ can lead to a decrease in r. RMSE$\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ decreases with height, and the highest median RMSE values were found over sea-level stations (20.7 for the WRF-Chem model and 17.3 for the EMEP model; Fig. 8, Table S1). It should be noted that RMSE depends greatly on the concentration magnitudes. Higher values of RMSE for both models generally correspond to the stations with low r values (the hotspot areas: southern Poland, the Czech Republic, Po Valley; Figs. S4–S5). Figure 8 shows that the trends of systematic errors differ between the models. The lowest errors in the WRF-Chem model were found over sea-level stations, while the highest were found over elevated stations. The errors in the EMEP model were comparable at all altitudes; however, the range of errors increased with height. A similar performance was found for the EMEP model and unsystematic errors. The median values were comparable at all altitudes, while the range slightly increased with height. In the case of the WRF-Chem model, a moderate increase in the median and the range of unsystematic errors with height was found. The areas affected with increased NMSE_unsys($\left(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}\right)$ were the hotspot areas (Po Valley and southern Poland) in the EMEP model, while in the WRF-Chem model, the increase in NMSE${}_{\mathrm{unsys}}\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ is found at almost all stations, particularly at the mountain level (Figs. S4–S5). It must be pointed out that both NMSE_sys and NMSE_unsys of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ in the EMEP model were substantially smaller at all altitudes with respect to the WRF-Chem model.

Figure 8 Intercomparison of the applied statistical measures (BIAS, IOA, r, RMSE, NMSEsys, NMSEunsys) between measured $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ (310 rural background stations from AirBase, https://www.eea.europa.eu/data-and-maps/data/airbase-the-european-air-quality-database-7, and the EU-PHARE project) and modelled $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ with the WRF-Chem (red boxes) and EMEP (blue boxes) models during November 2011 with respect to the station height.

The overall performance of the models regarding $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ is in agreement with similar modelling studies (e.g. Werner et al., 2015; Baró et al., 2015; Forkel et al., 2015; Gauss et al., 2016). Due to the coarser grid resolutions, differences in terrain height could lead to a problem in station representativeness in regional models. Generally, from the analysis given, it can be concluded that the performance of both models varies with height. There is moderate agreement in all of the analysed meteorological parameters and $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$, which shows a trend in the decrease in performance with height. This can be seen in Figs. 7–8. The better modelling performance was found for $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ using the WRF-Chem model, while it was found for $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_d$ in the case with the EMEP model. Both systematic and unsystematic errors for $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ were the lowest for sea-level stations and at comparable levels between models. Values of $r\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ and RMSE $\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ decreased with height for both models. A substantial number of elevated stations are located in the vicinity of hotspot areas (south Poland, the Czech Republic, etc.; Figs. S4, S5) and are therefore subject to strong influence from high emissions. This can explain the relatively lower performance (e.g. NMSE${}_{\mathrm{sys}}\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ for the WRF-Chem model; RMSE$\left(\right(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}})$ for both applied models) of a number of stations at an elevated level with respect to other altitudes in this area.

3.3.2 Analysis of model performance during the large-scale episodes

Here we focus on the analysis of spatial and temporal variations in the mean surface daily fields ($(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$, and $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ with $(\stackrel{\mathrm{‾}}{\mathrm{wd}}{)}_{\mathrm{d}}$) between the two applied models in order to investigate the mechanisms behind the high-pollution episodes.

Figure 9 Modelled $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ as “Conc”, and $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ as pressure, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ as “Temp”, $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ as PBLH, and $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ with $(\stackrel{\mathrm{‾}}{\mathrm{wd}}{)}_{\mathrm{d}}$ as WS for 2 typical days during the first (9 November 2011) and second (14 November 2011) high-pollution episodes from the WRF-Chem and EMEP models.

In Fig. 9, the modelled surface $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ together with $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$, and $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ with $(\stackrel{\mathrm{‾}}{\mathrm{wd}}{)}_{\mathrm{d}}$ for the 2 d with peak $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations (9 and 14 November 2011) during the two high-pollution episodes obtained with the EMEP and WRF-Chem models are shown. The distribution of $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ for both selected days was generally equal over the entire domain for both models. The $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ tends to have lower values (< 100 m) in the WRF-Chem simulation, and gradients in the pressure fields are much higher in comparison with the EMEP model. Values of $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ were generally higher within the domain for the WRF-Chem simulation. However, both models indicated the same areas with lowered wind speed, which is in accordance with the measurements (Fig. 5). Generally, both models correctly indicated areas affected by high-pollution episodes (DF > 100 %, Figs. 3–4). Over areas with $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ below 100 m, peaks of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ were found, reaching measured $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values (Fig. S6). For both peak days the models are consistent, showing prevailing high $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ fields, relatively cold areas with low $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ (more evident in the case of the WRF-Chem model) and low $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ conditions (more evident in the EMEP model) over the affected areas with $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations (Figs. 3–4). The Tables S1–S2 show the minimum, maximum and median values of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$, $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ and $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ over the domain (Fig. 1) for both models during the episodes. Minimum, maximum and median values of $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ between models were similar. The average minimum $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ over domain was 1004.77 and 1005.55 hPa, the average maximum was 1031.93 and 1031.44 hPa, and the average median was 1021.18 and 1020.33 hPa for the WRF-Chem and EMEP models, respectively. The average minimum $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ for WRF-Chem ($\sim -\mathrm{5.54}$ ^∘C) was lower with respect to the EMEP model ($\sim -\mathrm{2.31}$ ^∘C); however, the average maximum $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ (∼20 ^∘C) and median $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ (∼10 ^∘C) values were the same for both models. $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ in the WRF-Chem model varied from an average minimum value of 38.97 m to an average maximum value of 1612.29 m, while EMEP had a much higher average minimum value of 137.62 m (due to a coarser vertical resolution of the EMEP model) and a somewhat lower average maximum value of ∼ 1585.81 m (Tables S1–S2). $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ was more variable over the domain for WRF-Chem with respect to the EMEP model. During both episodes, minimum $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ in WRF-Chem was in the range of 0 to 0.11 m s⁻¹, while the maximum varied from 19.77 up to 36.34 m s⁻¹; the average median $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ was 5.00 m s⁻¹. For the EMEP model, minimum $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ was similar to WRF-Chem and in the range of 0.01 to 0.18 m s⁻¹, while maximum $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ was lower than that obtained with the WRF-Chem simulation: in the range of 12.74 to 16.77 m s⁻¹. The same applied to the average median $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, which was lower than in the WRF-Chem simulation: 3.60 m s⁻¹. The average $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations were generally higher in the EMEP model. The average minimum $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations were between 0.19 and 1.51 µg m⁻³, average maximum $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ was 62.04 and 84.45 µg m⁻³, and average median $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ values were between 6.91 and 13.46 µg m⁻³ for the WRF-Chem and EMEP models, respectively, during both episodes. The absolute maximum concentration obtained with the WRF-Chem model was 63.55 and 81.32 µg m⁻³, while for the EMEP model, it was 110.09 and 97.84 µg m⁻³ during the first and second episode, respectively.

During the first episode, the presence of cyclone Rolf in the Gulf of Genoa was evident in both models (Figs. S7–S8). Stronger surface winds occurred in the WRF-Chem simulation over Europe compared to the EMEP simulation, which consequently resulted in different dynamics within the boundary layer (Fig. S9). In the EMEP model, the onset of the high-pollution event was in central Europe as shown in the measurements, but with lower concentrations with respect to the measurements (Figs. 3 and S7). With NE winds over the coastal areas of northern Europe, the pollution gradually spread to western Europe. In the WRF-Chem model, the higher surface wind speed over central Europe was well estimated (Figs. 5, S7) and surface wind speeds over coastal areas in northern Europe were well-represented in the second part of the episode, leading to a good estimation of potential transport of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ to western Europe (Fig. S8). This agrees with similar studies where the dependence of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ on BIAS$\left(\right(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}})$ was identified (e.g. Solazzo et al., 2012). During both episodes, the $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ on the synoptic-scale was correctly predicted by both models over the domain (Figs. 3, 4, S7, S8, S10, S11). Aside from the $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$, notable differences between models performances were found in $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ (up to 200 m) and $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ (up to 5 ^∘C), which had an impact on the distribution and magnitude of the estimated high $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations in both episodes. In simulations with the WRF-Chem model (Fig. S11), the onset of the second episode was delayed by up to 1.5 d in comparison with the measurements (Fig. 4). Moreover, in the second episode, over areas with increased concentrations in central Europe, the decrease in $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ followed by a weak wind speed was found to be in accordance with the measurements (Figs. 5, right; S12). Recognised statically stable conditions (elaborated in Sect. 3.3.3.) with the presence of colder days prevailed over all of Europe. This favoured the build-up of concentrations in northwest and central Europe affecting all of central Europe (Figs. S8–S9). The representation of meteorological conditions over the affected areas (DF > 100 %, Figs. 3–4) agreed well with measurements during both episodes (Figs. 3–5, Figs. S7–S12). Although differences in $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ were found between the models (Figs. S9–S12), the areas with increased $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ were appropriately similar. However, as previously pointed out, the models underestimated the measured surface concentrations (Figs. S6–S8, S10, S11).

3.3.3 Intercomparison of modelled PBLH against radio soundings

More detailed analyses of model results and the influence of meteorological parameters during the second episode were made against measurements within the area of the Pannonian Basin (Fig. S3). The Pannonian Basin endured high-pollution events during the second high-pollution episode that were mainly found at urban stations (not shown) due to the lack of rural background measurements. In the analysed period, increased values of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ can be depicted only on one available rural background station in the area (e.g. Fig. S2, Fig. 1). The increased concentrations can also be observed from modelling results (Figs. 9, S2). The area of increased concentrations is in accordance with the area of weak wind conditions (Fig. 5) and low $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ values and can be described as an area of potentially statically stable conditions. The mean modelled vertical profiles during episodes at all available sounding stations within the area of interest agreed well with the measurements, except at the Belgrade station, where both models underestimated wind speed by up to 10 m s⁻¹ in the first 2000 m (Fig. S13).

Figure 10 Time series of the vertical profile of the bulk Richardson number (Eq. 9, the colour bar on the right) for the Zagreb, Budapest, Szeged and Belgrade sites from the WRF-Chem and EMEP models and sounding measurements (source: http://www.weather.uwyo.edu/upperair/sounding.html; last access: 12 May 2020) before/after and during second pollution episode (from 11 to 21 November). The black line indicates the boundary layer height (Ri_B=0.25).

Using sounding measurements and Eq. (9), the Ri_B and boundary layer height (H_bl) were calculated for four sites within the Pannonian Basin (Fig. S3) and are shown in Fig. 10. The same parameters were calculated from the WRF-Chem and EMEP modelling data. It must be pointed out that available sounding measurements were instantaneous values at 00:00 UTC only, while the time step in the WRF-Chem model was 1 h and in EMEP 3 h. The Ri_B values calculated from soundings and modelled data shown in Fig. 10 are represented with the same time step as input data: 24 h for measurements, 1 h for WRF-Chem and 3 h for the EMEP model. According to Fig. 10, the models were consistent in Ri_B and in estimating H_bl. The development of the atmospheric boundary layer started early in the morning with sunrise and reached values up to 350–400 m around 14:00 (local time), except between 17 and 21 November when a decrease in H_bl was found. During this period the peak values of H_bl reached 200 m and the statically stable conditions (Ri_B > 0.25) were dominant (light blue to dark blue colour up to value of 2; above that: white colour). As a strong increase in statically stable conditions occurred at all four stations, which are spread out within the basin, it can be concluded that statically stable conditions prevailed over the Pannonian Basin during this particular event. A similar conclusion comes from values of Ri_B and H_bl calculated from soundings (Fig. 10). Due to the coarse vertical resolution in some periods and the low time step (24 h), the contours are rough and the effect of sunrise on the development of H_bl cannot be seen. However, at all four measurement stations with prevailing statically stable conditions during the second high-pollution episode were indicated, which is in accordance with the modelling results.

4 Summary and conclusions

Numerical modelling of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ with different AQMs is still challenging (Baró et al., 2015; Prank et al., 2016; Laurent et al., 2016). It is therefore important to further analyse the different performances of regional models that have been widely used in practical applications. The main task of the current work was to investigate one of the weakest model capabilities, i.e. the simulations of AQMs under statically stable boundary conditions (e.g. Gašparac et al., 2016; Grisogono and Belušić, 2008) focusing on dynamic model aspects during episodes of elevated $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations over central and eastern Europe. Here, two different regional AQMs, namely, EMEP and WRF-Chem, were applied to evaluate their individual state-of-the-art performance and to investigate the processes that contributed to a high $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentration during pollution episodes that occurred in Europe. Other model intercomparison research studies over Europe and north America were done within the AQMEII project (e.g. Im et al., 2015; Solazzo et al., 2012; Rao et al., 2011). However, with respect to those large exercises with harmonised input data (same meteorology, emissions, boundary and initial conditions), the focus of this research was on the specific meteorological situations when statically stable atmospheric conditions prevailed accompanied by the occurrence of high $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations. The offline EMEP and online WRF-Chem modelling systems were used with the available input data that are usually implemented in practical applications (e.g. environmental assessment studies). The added value here is in the individual statistical evaluation of such modelling systems using data from the large number of meteorological and air quality stations in eastern Europe that have been less represented in other similar exercises. The analysed and modelled meteorological parameters were validated using surface measurements from 920 synoptic stations, soundings within the Pannonian region and mast-mounted instrument measurements. The $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations were validated against surface measurements from 310 rural background stations. During the colder part of the year, when usually higher PM concentrations are observed, the following model features are established.

• According to the low systematic errors a very good model performance is found in simulating $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ over sea-level and elevated stations, while there is a moderate performance due to the high spread of BIAS over mountain stations.

• Good performance in modelling $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ in EMEP and $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ in the WRF-Chem model is found, while, by contrast, the EMEP model highly overestimated $(\stackrel{\mathrm{‾}}{t_{\mathrm{2}\mathrm{m}}}{)}_{\mathrm{d}}$ and WRF-Chem overestimated $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$.

• The differences in boundary layer dynamics were found in models through the analysis of vertical wind profiles.

• Based on calculated values of Ri_B, the evaluation of modelled $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ agreed well with the measurements for both models. However, according to the spatial $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ fields, derived directly from model, the WRF-Chem model generally tends to estimate lower $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ with respect to the EMEP model over areas affected by high pollution (DF > 100 %).

• From the results of the simulation of a 1-month period that encompassed various meteorological conditions and different terrain types, we found the following:

  • Strong influence of meteorological conditions on increased background $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ and correct estimation of the $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ is recognised as one of the main factors in the dispersion of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$.

  • General underestimation of background $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations with both models, except with EMEP for mountain stations (slight overestimation).

  • Statistical analysis with respect to the terrain type shows the best modelling performance of $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ and meteorology over sea-level stations (flat terrain). Both models tend to agree on decrease in performance with height, indicating problems in regional model simulations over complex terrain.

• From the analysis of the high-pollution episodes, we can conclude the following:

  • During the first high-pollution episode, a high $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$ in the WRF-Chem model resulted in a decrease in surface $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$, while favourable conditions prevailed for the build-up of concentration in central Europe over hotspot areas with a decrease in surface $(\stackrel{\mathrm{‾}}{\mathrm{ws}}{)}_{\mathrm{d}}$.

  • Low wind speed conditions during the entire second episode, followed by high $(\stackrel{\mathrm{‾}}{\mathrm{mslp}}{)}_{\mathrm{d}}$ and low $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$, prevailed over the affected area (DF > 100 %).

  • Statically stable conditions were recognised as the main mechanism for the build-up of concentrations during the second episode. Both models produced low values of $(\stackrel{\mathrm{‾}}{\mathrm{pblh}}{)}_{\mathrm{d}}$ (< 100 m in WRF-Chem and 100–200 m in EMEP) over areas where stations recorded $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations > 200 % (DF) with respect to the annual mean (Figs. 3–4, S7, S8, S10, S11).

  • Underestimation of background $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations with regional models is in accordance with other modelling studies (Gauss et al., 2016; Forkel et al., 2015).

• Reasons for the underestimation of modelled $(\stackrel{\mathrm{‾}}{{\mathrm{PM}}_{\mathrm{10}}}{)}_{\mathrm{d}}$ concentrations were attributed to the uncertainty of associated and inadequate treatments of formation processes that usually omit some components of atmospheric aerosols (e.g. SOA, SIA) and thus fail to estimate the total PM budget properly.

Dynamic model properties are very important: horizontal and vertical model resolutions and the boundary layer parameterisations in statically stable atmospheric conditions should be selected carefully. Furthermore, model simulations using a more accurate emission inventory and larger (nested) domains with a finer resolution are necessary for further improving the model predictions. Future work using longer periods of simulations for both models, including other pollutants (e.g. NO_x, SO_x, PM compounds and O₃) is recommended to make comparisons under various meteorological conditions.