MACC regional multi-model ensemble simulations of birch pollen dispersion in Europe

. This paper presents the ﬁrst ensemble modelling experiment in relation to birch pollen in Europe. The seven-model European ensemble of MACC-ENS, tested in trial simulations over the ﬂowering season of 2010, was run through the ﬂowering season of 2013. The simulations have been compared with observations in 11 countries, all members of the European Aeroallergen Network, for both individual models and the ensemble mean and median. It is shown that the models successfully reproduced the timing of the very late season of 2013, generally within a couple of days from the observed start of the season. The end of the season was generally predicted later than observed, by 5 days or more, which is a known feature of the source term used in the study. Absolute pollen concentrations during the sea-Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
During the last 30 years, the prevalence of airborne allergy and asthma in Europe has increased 4-fold, reaching 15-40 % of the population. According to the European Federation of Allergy and Airway Diseases Patients Associations, 80 million (24.4 %) adults living in Europe are allergic. The allergy prevalence in children is 30-40 % and increasing (Laatikainen et al., 2011;Rönmark et al., 2009). Allergy to various types of pollen in the air, exacerbated by co-exposure to chemical pollutants and fine aerosols, is the number one chronic disease in Europe, overshadowing allergy to house dust mite and affecting over 20 % of the population (Bousquet et al., 2007).
Among the allergenic plants, grass and birch pollen affect about 40 and 25 % of all hay fever sufferers in Europe respectively (Heinzerling et al., 2009). Birch is a strong allergy-provoking tree with a population-wide sensitisation of approximately 15 % (WHO, 2003). The distribution of silver birch (Betula pendula Roth.) and downy birch (B. pubescens L.) extends from mountains in the temperate climate of southern Europe to Fennoscandia and Siberia (Atkinson, 1992;OECD, 2003).
It has long been known that the bulk of pollen is deposited near the source plant (Raynor et al., 1970;Tampieri et al., 1977;Wright, 1953Wright, , 1952. However, birches, as well as other species (Alnus, Carpinus, Corylus, Ostrya, Fagus, Quercus, Castanea) belonging to the order Fagales, are windpollinated trees generating vast amounts of pollen to ensure a sufficient level of fertilisation of female flowers over receptor regions. Their pollen grains are sufficiently small and light to facilitate the atmospheric transport of up to 1 % of the released material over thousands of kilometres when weather conditions are suitable (Sofiev et al., 2006a). This phenomenon was noticed in the middle of the twentieth century (Erdtman, 1937(Erdtman, , 1935(Erdtman, , 1931Gregory, 1961). Later, it was recognised that the long-range transported pollen can have a substantial health impact (Viander and Koivikko, 1978) and facilitate a large-scale redistribution of genetic material (Lindgren et al., 1995). Long-range transport of pollen is practically unpredictable with local observations or statistical models. However, up until the last 2 decades no practical instruments had been developed for its quantitative evaluation and forecasting.
The main tool for analysing the pollen distribution on regional and continental scales is numerical modelling that combines phenological models covering pollen maturation and presentation (the pollen source term) with the atmospheric dispersion model. Currently, there exist four comparatively independent formulations of the source terms for birch pollen. The European-scale source term used in the current study was developed for the SILAM model (http://silam. fmi.fi) by an international consortium within the scope of the POLLEN project of the Academy of Finland (Siljamo et al., 2012;Sofiev et al., 2012a). Various versions of the model have been used for forecasts of pollen distribution in Europe starting from 2005 (Sofiev et al., 2006a) and reanalysis of the flowering seasons back to 1997 (Siljamo et al., 2008c;Veriankaitė et al., 2010). The COSMO-ART birch module was developed at the University of Karlsruhe (Helbig et al., 2004;Vogel et al., 2008) and MeteoSwiss (Pauling et al., 2012;Zink et al., 2013) and is currently used for pollen forecasting for central and south-western Europe. Development is also going on in Denmark with the regional ENVIRO-HIRLAM system (Mahura et al., 2009) applied for forecasting over northern Europe. Finally, combining the COSMO-ART and SILAM source terms, Efstathiou et al. (2011) developed a regional-scale model for the USA and applied it to birch and ragweed.
MACC (Monitoring of Atmospheric Composition and Climate, http://www.gmes-atmosphere.eu) pollen simulations are based on the SILAM source term. Its formulations and input data have been shared among the seven regional modelling teams of MACC and, in co-operation with the European Aeroallergen Network (EAN), set into operational multi-model ensemble forecasting of birch pollen in Europe.
The goal of the current paper is to present and evaluate the results of the first ensemble modelling of birch pollen in Europe during the season of 2013.
The next section will present the models and setup of the simulations as well as the observation data used for evaluation of the model predictions. The results section will present the outcome of the simulations and the quality scores of the individual models and the ensemble. The discussion section will be dedicated to analysis of the results, considerations of the efficiency of the multi-model ensemble for pollen, and identification of the most pressing development needs.

Materials and methods
This section presents the regional models used in the study, outlines the birch pollen source term implemented in all of them, and introduces the pollen observations used for evaluation of the model predictions.

Dispersion models
The dispersion models used in the study comprise the MACC-II European ensemble, which is described in detail in Marécal et al. (2015). In the following, only the model features relevant for the pollen atmospheric transport calculations are described.
The ensemble consisted of seven models.
-CHIMERE (http://www.lmd.polytechnique.fr/chimere/) is an Eulerian regional-scale chemistry-transport model for gaseous and aerosol species (Menut et al., 2013). Pollen is implemented as a special aerosol with a prescribed species-specific size (currently birch or ragweed) between 20 and 22 µm and density of 800 and 1050 kg m −3 respectively. The resulting gravitational settling velocity is 1.2-1.3 cm s −1 . The transport processes affecting pollens, such as advection, turbulent mixing, and wet deposition, are implemented in the same way as for other aerosols. Dry deposition is described via gravitational settling only, which dominates for pollens (Sofiev et al., 2006a), whereas resuspension is parameterized following Helbig et al. (2004).
-EMEP model of EMEP/MSC-W (European Monitoring and Evaluation Programme/Meteorological Synthesizing Centre -West) is a chemical transport model developed at the Norwegian Meteorological Institute and described in Simpson et al. (2012). It is flexible with respect to the choice of projection and grid resolution. Dry deposition is handled in the lowest model layer. A resistance analogy formulation is used to describe dry deposition of gases, whereas for aerosols the mass-conservation equation is adopted from Venkatram (1978) with the dry deposition velocities dependent on the land use type. Wet scavenging is dependent on precipitation intensity and is treated differently within and below clouds. The below-cloud scavenging rates for particles are calculated based on Scott (1979). The rates are size dependent, growing for larger particles.
-EURAD-IM (http://www.eurad.uni-koeln.de) is an Eulerian mesoscale chemistry transport model involving advection, diffusion, chemical transformation, wet and dry deposition, and sedimentation of tropospheric trace gases and aerosols (Hass et al., 1995;Memmesheimer et al., 2004). It includes 3D-VAR and 4D-VAR chemical data assimilation (Elbern et al., 2007) and is able to run in nesting mode. The positive definite advection scheme of Bott (1989) is used to solve the advective transport and the aerosol sedimentation. An eddy diffusion approach is applied to parameterize the vertical sub-grid-scale turbulent transport (Holtslag and Nieuwstadt, 1986). Dry deposition of aerosol species is treated as size dependent using the resistance model of Petroff and Zhang (2010). Wet deposition of pollen is parameterized according to Baklanov and Sorensen (2001).
-MATCH (http://www.smhi.se/en/research/researchdepartments/air-quality/match-transport-andchemistry-model-1.6831) is an Eulerian multi-scale chemical transport model with mass-conservative transport and diffusion based on a Bott-type advection scheme (Langner et al., 1998;Robertson and Langner, 1999). For birch pollen, dry deposition is mainly treated by sedimentation and a simplified wet scavenging scheme is applied. The temperature sum from March onwards, driving the birch pollen emission, is determined outside the model and fed into the emission module.
-MOCAGE (http://www.cnrm.meteo.fr/gmgec-old/site_ engl/mocage/mocage_en.html) is a multi-scale dispersion model with grid-nesting capability (Josse et al., 2004;Martet et al., 2009). The semi-Lagrangian advection scheme of Williamson and Rasch (1989) is used for the grid-scale transport. The convective transport is based on the parameterization proposed by Bechtold et al. (2001) whereas the turbulent diffusion follows the parameterization of Louis (1979). Dry deposition including the sedimentation scheme follows Seinfeld and Pandis (1998). The wet deposition by the convective and stratiform precipitations is based on Giorgi and Chameides (1986).
The dry deposition scheme (Kouznetsov and Sofiev, 2012) is applicable for a wide range of particle sizes including coarse aerosols, which are primarily removed by sedimentation. The wet deposition parameterization distinguishes between sub-and in-cloud scavenging by both rain and snow (Sofiev et al., 2006b). For coarse particles, impaction scavenging is dominant below the cloud. The model is capable of 3D-and 4D-VAR data assimilation (Vira and Sofiev, 2012), also applicable to birch.
-ENSEMBLE models were generated by the arithmetic average and median calculated from seven model fields for each hour.

Birch pollen source term
All models of this study are equipped with the same birch pollen source term  verified for the season of 2006 by Siljamo et al. (2012). The formulations and input data are open at http://silam.fmi.fi/MACC. The main input data set is the birch habitat map compiled by Sofiev et al. (2006a) with a spatial resolution of 0.5 • × 0.5 • longitudelatitude. The birch productivity is assumed to be the same in all years and equal to 10 9 pollen m −2 season −1 .
The flowering description follows the concept of thermal time phenological models and, in particular, the doublethreshold air temperature sum approach of Linkosalo et al. (2010) modified by , which determines the flowering propagation during the whole spring season. Within that approach, the heat accumulation starts on a day in spring (1 March in the current setup) and continues throughout the season. Flowering starts when the accumulated heat reaches the starting threshold and continues until the heat reaches the ending threshold. The rate of heat accumulation is the main controlling parameter for pollen emission: the model establishes direct proportionality between the flowering stage and fraction of the heat sum accumulated to-date.
Apart from temperature, the pollen release rate is modulated by ambient humidity, precipitation, and wind speed. Following , higher relative humidity (RH) and rain reduce the release, completely stopping it for RH > 80 % and/or rain > 0.1 mm h −1 . Strong wind promotes it by up to 50 %. Atmospheric turbulence is taken into account via the turbulent velocity scale and thus becomes important only in cases close to free convection. In stable or neutral stratification and calm conditions the release is suppressed by 50 %.
Local-scale variability of the flowering results in the need to include probabilistic description of the flowering propagation (Siljamo et al., 2008b). In the simplest form, the probability of an individual tree entering the flowering stage can be considered via the uncertainty of the temperature sum threshold determining the start of flowering for the grid cell.
The end of the season is described via the open-pocket principle: the flowering continues until the initially available amount of pollen is completely released.

Pollen observations
The observations for the model evaluation in 2013 have been provided by the following 11 members of the EAN: Austria, Estonia, Germany, Finland, France, Latvia, Lithuania, Russia, Spain, Switzerland, and Ukraine. Additionally, the data for the initial model testing for the season of 2010 were provided by Austria, Finland, Germany, Latvia, Lithuania, Russia, Spain, Switzerland, and Ukraine. In total, information from 165 sites in 2010 and 186 sites in 2013 was made available to the modelling teams. Among these, 21 stations in mountain valleys of the Alps and Pyrenees were flagged as not representative on the regional scale and excluded from the analysis (see the Discussion section). The analysis below concentrates on the season of 2013 as the data for 2010 were mainly used for setting up and verifying the pollen source term implementations. However, a comparison of these years is used to illustrate the variability of pollen seasons and the ability of the models to reproduce it.
Pollen monitoring was performed with Burkard 7-day and Lanzoni 2000 pollen traps based on the Hirst design (Hirst, 1954). The pollen grains were collected at an airflow rate of 10 L min −1 . The observations covered the period from March until September, with some variations between the countries. Daily observations were used. Following the EAN recommendations (Galán et al., 2014;Jäger et al., 1995), most samplers were located at heights of between 10 and 30 m on the roofs of suitable buildings. The places were frequently in the cities' downtown areas; i.e. they largely represent the urbanbackground conditions (although not always). With regard to microscopic analysis, the EAN recommendation is to count at least 10 % of the sample using horizontal or vertical strips (Galán et al., 2014). The actual procedures vary between the countries but generally comply. The counting in 2013 was performed along 12 vertical strips (in most countries), or two to four horizontal traverses (Switzerland, Spain), using a bihourly stratified random sampling (Finland) (Mandrioli and Comtois, 1998). In all cases, the data were expressed as mean daily concentrations (pollen m −3 ).

Setup of the simulations
Simulations followed the standards of MACC regional ensemble (Marécal et al., 2015). The domain spanned from 25 • W to 45 • E and from 30 to 70 • N. Each of the seven models was run with its own horizontal and vertical resolutions, which varied from 0.1 to 0.25 • of the horizontal grid cell size, and had from 3 up to 52 vertical layers within the troposphere (Table 1). This range of resolutions is not designed to reproduce local aspects of pollen distribution, instead covering the whole continent and describing the large-scale trans- port events. The limited number of vertical dispersion layers is a compromise allowing for high horizontal resolution. Thick layers are not a major limitation as long as the full vertical resolution of the input meteorological data is used for the evaluation of dispersion parameters (Sofiev, 2002).

Atmos
In the forecasting regime during the spring of 2013, the time range of the simulations was 96 h from 00:00 UTC on day 0 (D0) with hourly output on eight vertical levels (surface, 50, 250, 500, 1000, 2000, 3000, and 5000 m above the surface). After the end of the season, it was reanalysed by most of the models to correct technical problems experienced in the forecasting regime. For the reanalysis simulations (discussed in this paper), the models were run through the whole period without separation into individual forecast cycles. For those models that were not rerun, the first 24 h of each forecast were used. In all cases, only near-surface concentrations were analysed.
All models considered pollen as an inert water-insoluble particle 22 µm in diameter and with density of 800 kg m −3 (Bassett et al., 1978;Bucher and Kofler, 2015;Sofiev et al., 2006a) 3 Results for the flowering season of 2013

Observed peculiarities of the season
The season of 2013 had three major specifics, which distinguished it from "typical" pollen seasons and, in particular, from the training year of 2010: -A cold spring resulted in late flowering. In central Europe, the flowering started up to 2-3 weeks later than usual. For instance, 2013 in Switzerland was among the latest years since 1993 (the latest at five stations): 9 days later than in 2010. In Moscow, the cold start of the spring was compensated by its faster progression, so that the early-flowering alder was shifted by about 2 weeks but the birch season was delayed by only a few days. In Lithuania, however, the observed highconcentration time period started 10 days earlier than in 2010, almost simultaneously with France, which was probably caused by early long-range transport events.
-The duration of the season was up to 1 week shorter than usual. Thus, in Switzerland it lasted for ∼ 30 days (22-35 for different stations) as compared the long-term average of 37 days. In Finland, the difference between the season length in 2010 and in 2013 reached a factor of 2.4. -An anomalously low pollen season was recorded in northern Europe and Russia. The seasonal pollen index (SPI: the sum of daily pollen concentrations over the whole season) was 10-1000 times lower than in 2012 and about 10 times weaker than in 2010 (that year was comparatively usual). The SPI in central Europe was moderate, which resulted in an inverse load pattern: the SPI in the north was several tens of times lower than that in the central regions.
These peculiarities presented substantial challenges to the models. The phenological model of the source term has a mechanism that accounts for the season shift, but it still went beyond the verified range. The season strength, however, is currently not a predicted quantity, which made it impossible to capture the anomalously low season in the north.

Model results
All models predicted a quite standard load pattern for the SPI (Fig. 1). Its maximum is located over central Russia and Fennoscandia and the SPI gradually decreases towards the south-west. In central Europe, there is a substantial inhomogeneity of the SPI, which reflects the patchy birch habitat in the region.
Comparison with the observed SPI shows the challenge mentioned above: the very low observed concentrations in the north that were not reproduced by the models (Fig. 2, model predictions are overlaid by the observations -circles coloured following the same palette). This is in contrast with the previous years, particularly 2010, when the observed and modelled patterns were both typical and agreed very well (see example for SILAM, Fig. 3).
An example of the hourly concentrations at noon on 20 April 2013 is shown in Fig. 4. It depicts the middle of the season in central and eastern Europe. The models also showed the long-range transport of pollen to the south that reached Africa in most predictions.
The progress of the season is illustrated in Fig. 5, which depicts 5-day mean concentrations predicted by the ensemble median for four episodes: 1-6, 20-25 April, and 10-15 May and 1-6 June. The season progress in 2013 was quite usual: from the south-west to the north-east of the continent though delayed by up to 2 weeks due to the cold slow spring. The models successfully reproduced this development.
The primary parameters describing the season are its start and end days, often defined as the dates when 5 and 95 % of the cumulative seasonal pollen counts are reached respectively (Fig. 6). Outside the main source areas, the timing of the season is almost completely dictated by the episodes of long-range pollen transport, similar to the southern transport  episode shown in Fig. 4. The "fingerprints" of the plumes bringing the first pollen to the regions and those concluding the season are clearly seen in southern Europe, where Spanish stations show a presence of pollen almost as long as in Finland.
With regards to the model-measurement comparison statistics, one has to bear in mind that the time series of pollen concentration represents strongly non-stationary and non-ergodic processes, i.e. the usual statistics (bias, RMSE, correlation, etc.) that all rely on the process stationarity and ergodicity can be computed only within the main season and even then have to be taken with care. A series of such "standard" statistics was computed for the ensemble (Fig. 7), as well as for all individual models (Fig. 8). For the ensemble quality assessment, we used a discrete rank histogram, which is a simple and efficient way to understand the basic properties of the obtained ensemble. The theoretical basis and examples of this histogram and more sophisticated approaches can be found in Candille and Talagrand (2005) and Potempski and Galmarini (2009).
All statistics were calculated for the whole season: 15 March-24 April inclusive. Daily-mean values were used for RMSE and correlation coefficient.

Discussion
Within this section, the following issues are considered: (i) the ability of the model ensemble to predict the key features of the 2013 birch pollen season, (ii) main uncertainties of the current ensemble, (iii) specific features of the individual ensemble members, and (iv) parameters of the season, for which the use of the ensemble predictions is more beneficial than single-model simulations.

Model predictions for the key season parameters
The most important parameter for the users of pollen forecasts is the season start. Analysis shows (see Figs. 6 and 7) that the ensemble captured the season onset over the majority of central and western Europe with an error of just a couple of days, which is very small. This is in agreement with the source term evaluation by Siljamo et al. (2012).
Both in the north (Finland and Baltic states) and in Spain the pattern is very irregular: the error of the season start (the date when 5 % of the seasonal total is reached) at stations located a couple hundred kilometres from each other can differ by more than a week (Fig. 6). The main reason for such inhomogeneity is that the season start over these areas was largely influenced by remote sources in central Europe and long-range transport. A single episode affecting or pollen cloud passing by the station can result in a few weeks of the apparent-season shift.
The end of the season is more uncertain: the concentrations usually fall more slowly at the season end than they grow at its start, with substantial small-scale variability unresolved by the large-scale simulations. As seen from Fig. 6, the error usually stays within some 5 days but can also reach several weeks, especially in the mountainous regions (Pyrenees, Alps). Fortunately, this parameter is also less important for practical applications.
Representation of the absolute concentrations strongly varies over the European continent (Fig. 7). In its central part Atmos. Chem. Phys., 15, 8115-8130, 2015 www.atmos-chem-phys.net/15/8115/2015/ a) 5% obs , b) 95% obs c) 5%, mdl d) 95%, mdl  (Germany, Austria, part of France), there is a slight underestimation. In southern France and Spain, it gradually turns to a slight overprediction, suggesting a somewhat too long transport distance in the majority of the models. Finally, in the north all models strongly overpredict. These tenden-cies are practically not dependent on the longitude: available observation points in the east follow the same pattern of very slight overprediction in Ukraine and substantial overstatement in central Russia. The RMSE largely follows the bias field (Fig. 7). Correlation coefficient (Fig. 7) should be taken with care due to the evident non-stationarity of the process. However, it also highlights central and western Europe (as well as part of northern Europe) as the best-predicted areas. Mountains are the most difficult regions, along with the areas with few birch stands (southern Europe), where the habitat map is highly uncertain. Northern Europe is usually a well-predicted area but not in 2013: as seen from Fig. 7, correlation of time series in the Baltic states and southern part of Finland is quite low. It is high only in Moscow and northern Finland.

Main uncertainties of the ensemble
From the above analysis, one can deduce the main sources of uncertainties of the presented multi-model ensemble: (i) missing interannual variability of the birch productivity, (ii) errors in the mountainous regions, and (iii) birch distribution map.
Currently, there is no model for year-to-year variation of the birch productivity. A few studies reported in literature -e.g. Masaka (2001), Ranta et al. (2008Ranta et al. ( , 2005) -concentrate on predicting the SPI, which is a different quantity significantly affected by the current-year meteorological conditions. The total amount of pollen stored in catkins, in contrast, is decided by the previous-year summer and, to some extent, the following winter conditions. The second complication is that the existing studies are based on a limited number of observation points, which makes it difficult to generalize them to the continental scale. The work is ongoing, but so far the only way to obtain a reliable absolute level of concentrations is via data assimilation performed retrospectively.
The large uncertainties in the mountains originate from the insufficient resolution of both meteorological and dispersion models. As an illustration, the time series for Zams station in Austria (one of the 21 sites excluded from the comparison) shows that all models have shifted the season by several weeks: mid-June instead of late April to early May (Fig. 9). Some models also predicted peaks shortly before the season but these were the pollen plumes from remote sources. This error is exacerbated by a strong underestimation of the absolute values. The reason for the poor performance of all models is the complex-terrain environment with a characteristic width of the valley of barely 2 km. Continental-scale dispersion models, as well as the global meteorological model, all have a resolution 10-20 times coarser than that (Table 1). As a result, the grid-cell-scale temperature is not representative of the valley bottom (it is biased low), which leads to a late predicted start of the season. Moreover, pollen released at the bottom of the narrow valley is usually trapped inside it -in reality -whereas in the models it is mixed over the whole grid cell, which leads to strong underestimation even though the total released amount of pollen is reasonable.

Behaviour of individual models
Following the MACC standards, the setup of the ensemble members was largely harmonised. The emission term was implemented in all models with minimal differences from the one described in . The only known deviation was made in SILAM, where the pollen maturation and discharge were separated into two processes controlled by different environmental parameters, similar to Prank et al. (2013) and Zink et al. (2013). However, the impact of this variation at daily averaging is bound to be small. The other parameter that depended on the models was the thickness of the emission injection layer. The recommended layer was from the surface up to a height of 50 m but the model's geometry affected it. Meteorological input data were the same -the IFS forecasts. Meteorological data pre-processing was based on comparatively simple diagnostic procedures embedded in all models except for EURAD, which used the WRF model nested into IFS. The tasks of the preprocessors were (i) to derive the boundary layer characteristics missing from the IFS and (ii) for some models to re-diagnose the vertical wind component or refine the 3-D wind fields to ensure satisfaction of the continuity equation.
In light of this, the differences between the model predictions visible in Figs. 1, 4, and 8 should be mainly attributed to (i) model treatment of the 3-D pollen transport, (ii) vertical mixing, and (iii) removal mechanisms.
In general, the model results are quite similar and the main features of the pollen distribution are clearly visible in the individual-model patterns. On a closer look, one can see a few tendencies, such as (i) higher-than-average concentrations predicted by EURAD, (ii) lower-than-average values of MOCAGE, (iii) longer lifetime and farther atmospheric transport of CHIMERE, (iv) the shortest transport distance indicated by SILAM (Fig. 1), (v) about 10 % lowerthan-others correlation coefficient of EMEP, EURAD, and LOTOS-EUROS (Fig. 8), and (vi) a general tendency of a couple of days early start of the season for most models except for CHIMERE (6 days too early) and for EMEP and SILAM (1-2 days late season) (Fig. 8).
Interpretation of these tendencies is not unequivocal but some of them are connected. For instance, the season start is largely controlled by the possible long-range transport episodes, so that the model reporting the longest pollen transport (CHIMERE) predicts a too-early season start, especially in remote regions. Conversely, SILAM, with its shortest travelling distance, would report fewer such events, leading to later season onset. The same is true for EMEP, which also reported quite a short transport distance. The high predictions of EURAD lead to somewhat higher RMSE and low predictions of MOCAGE lead to lower RMSE, owing to the high bias of all the models in the north.
In a few cases, available observational information enables evaluation of these features. In particular, a quick north-tosouth reduction of the observed SPI in Spain (Fig. 2) suggests that the transport distance of pollen is indeed short (see also Fig. 3). This would also reduce the early bias of the season start shown by most models. Secondly, the approximately 10 % higher correlation coefficient of SILAM might possibly be attributed to the more articulated impact of local birch since the long-range effects are of lower importance for that model, owing to its short transport distance. These observations might also be projected to the resuspension process suggested by Helbig et al. (2004) but never explicitly verified or confirmed, as pointed out by . Among the seven models, only CHIMERE included it, which may be one of the reasons for its longer transport distance. The above comparison raises further concerns regarding this process. Evaluation of absolute concentrations is hardly conclusive due to the very specific season and overall uncertainty of this parameter.

Ensemble added value
Compilation of the multi-model ensemble out of individual models has proven to be beneficial. As seen from Fig. 8, the ensemble median, together with the ensemble mean and SILAM, has the highest correlation coefficient. The ensemble also showed among the shortest shifts of the season start and among the smallest RMSEs. Its robustness to outliers also turned out to be a very strong asset for pollen forecasting, which is a comparatively new area of modelling and has many unknowns in the governing processes.
However, the ensemble can only be as good as the majority of the individual models. As follows from the rank histogram (the Talagrand diagram, Fig. 10, left panel), the current ensemble is not perfect. The diagram shows the generally underestimating ensemble (tendency towards higher ranks of the observations) simultaneously with a fraction of observations being substantially overestimated (zeroth rank is also frequent). The latter feature is due not only to the northern sites (albeit few) but also to southern Spanish stations, where the overestimation is also systematic. Moderate underestimation takes place mainly in central Europe (Fig. 2), where most of the stations are located. However, even this imperfect ensemble outperformed most of the individual models.
Construction of the rank histogram for pollen faced a methodological problem: pollen concentrations long before and long after the flowering season are 0. One should also bear in mind the quite high detection limit for the microscopic analysis (about 0.5 pollen m −3 for daily mean), which increases the frequency of observed zero concentrations. In 2013, about 15 % of observations and 10-15 % of model predictions were below this limit. For the all-zero cases (zero observations and below-detection-limit bulk/all of the models), determination of the observation rank is meaningless: the ensemble variance collapses and its mean matches the observations perfectly. Neither of these cases can be ignored: they have clear physical meaning and manifest excellent ensemble behaviour. Therefore, for such cases the observation rank in Fig. 10 was picked as a random number from 0 to 7, i.e. it corresponded to a "perfect" ensemble. To illustrate the ensemble behaviour when the observations were below the detection limit, the histogram of the corresponding model values is shown in Fig. 10, right panel. As one can see, in over 80 % of zero observations models also showed the concentrations below 0.5 pollen m −3 . This peculiarity has to be kept in mind for future pollen-related ensembles.

Summary
The first-ever multi-model ensemble was created for predicting the concentrations of birch pollen in Europe. The ensemble was constructed from seven European atmospheric chemistry transport models of MACC and follows the main rules of the MACC regional ensemble: the models share the same source term, use the same meteorological input, and cover the MACC domain with similar resolution.
The ensemble was evaluated against the observational data of European Aeroallergen Network for the year 2013 and the basic statistical indicators were computed for each individual ensemble member, as well as for the ensemble mean and median.
Atmos. Chem. Phys., 15, 8115-8130, 2015 www.atmos-chem-phys.net/15/8115/2015/ The ensemble demonstrated good skills in predicting several main characteristics of the pollen season of 2013: the season start and propagation, the pollen distribution pattern in central, eastern, and southern Europe, and characteristic concentrations over these regions. The season timing was captured despite the anomalously late flowering due to the cold spring of 2013.
Representation of the pollen concentrations in northern Europe, the Baltic states, and central Russia was affected by the anomalously low flowering intensity in 2013. As a result, all models had strongly overestimated pollen levels there. This was in contrast to the usual pollen distribution pattern in Europe, such as the one of the quite typical year of 2010, which was reproduced much better.
The experiment showed the high added value of the ensemble. For most of the participating models this was the first experience of pollen simulations, which affected the reliability of their results. The ensemble median proved to be robust to the outliers, finally showing among the highest correlation coefficients and one of the smallest errors in the season timing and RMSEs.
The main areas of improvement referred to the interannual variation of the birch productivity as well as to the representation of the flowering timing in the complex-terrain conditions.