Incorporation of pollen data in source maps is vital for pollen dispersion models

Information about distribution of pollen sources, i.e. their presence and abundance in a specific region, is important, especially when atmospheric transport models are applied to forecast pollen concentrations. The goal of this study is to evaluate three pollen source maps using an atmospheric transport model and study the effect on the model results by combining these source maps with pollen data. Here we evaluate three maps for the birch taxon: (1) a map derived by combining a land cover data and forest inventory, (2) a map obtained from land cover data and calibrated using model simulations and pollen observations, and (3) a statistical map resulting from analysis of forest inventory and forest plot data. The maps were introduced to the Enviro-HIRLAM (Environment – High Resolution Limited Area Model) as input data to simulate birch pollen concentrations over Europe for the birch pollen season 2006. A total of 18 model runs were performed using each of the selected maps in turn with and without calibration with observed pollen data from 2006. The model results were compared with the pollen observation data at 12 measurement sites located in Finland, Denmark, and Russia. We show that calibration of the maps using pollen observations significantly improved the model performance for all three maps. The findings also indicate the large sensitivity of the model results to the source maps and agree well with other studies on birch showing that pollen or hybrid-based source maps provide the best model performance. This study highlights the importance of including pollen data in the production of source maps for pollen dispersion modelling and for exposure studies.

1 Introduction

Aeroallergens are a specific type of atmospheric aerosols causing allergic reactions among people suffering from allergic rhinitis, which is often connected with asthma (Bachert et al., 2004). The number of allergic patients sensitive to pollen is assessed to be 20 % of the European population (WHO, 2003). There are substantial variations in sensitization levels towards specific aeroallergens in Europe (Heinzerling et al., 2009). In northern Europe, pollen from the Poaceae (grasses) and Betulaceae (e.g. hazel, alder, and birch) families show the highest sensitization among patients (Heinzerling et al., 2009). Birch pollen is generally the most abundant of those four pollen types (Skjøth et al., 2013b), often with large interannual fluctuations in the overall pollen integral (e.g. Piotrowska and Kaszewski, 2011). In particular, for northern Europe the pollen season of 2006 was characterized by abundant birch pollen concentrations that resulted in an increase in patient calls for medical assistance during the spring period (Sofiev et al., 2011). Such episodes can potentially be described and forecast using atmospheric dispersion models (Sofiev et al., 2006), helping sensitized individuals to minimize their symptoms. Several European countries, such as Denmark, have seen an increase in sensitized individuals over the past few decades (e.g. Rasmussen, 2002; Linneberg et al., 2007). Generally, skin prick testing (SPT) is performed to identify sensitization of patients to common pollen allergens (e.g. Heinzerling et al., 2009). Haahtela et al. (2015) showed that the percentage of the population with positive SPT was about 35 % and 20 % in Finnish and Russian Karelia, respectively. These facts demonstrate the importance of pollen studies using atmospheric dispersion models applied for research and operational purposes (Klein et al., 2012).

A variety of pollen types are considered to be biological components by several atmospheric dispersion models such as SILAM, COSMO-ART, WRF/CMAQ, Enviro-HIRLAM, and others (see e.g. Zink et al., 2012, 2013; Zhang et al., 2014; Sofiev et al., 2015a; Baklanov et al., 2017; Sofiev, 2017; Zink et al., 2017). Information about distribution of pollen sources, i.e. their presence and abundance in a region of study, plays a crucial role in aerobiological modelling and forecasting using atmospheric dispersion models (e.g. Zink et al., 2017). In other words, it is one of the key input data used to simulate or forecast pollen emission, its atmospheric transport, and depositions in these models. For this reason, recent studies (e.g. Pauling et al., 2012; Bonini et al., 2018) are devoted to developing pollen source inventories and maps. Two main approaches aiming to build/construct the pollen source inventories are suggested in the literature: bottom-up and top-down (Skjøth et al., 2013b). The bottom-up approach is based on statistical data of pollen source distribution combined with land cover data obtained from remote sensing (Skjøth et al., 2008a). The limitation of this approach is that statistical data are either unavailable or cover too large of a territory (e.g. the whole country) in some regions. The top-down approach uses pollen observations as a starting point to estimate the abundance of different pollen taxa and it can be combined with land use data and/or models for more detailed assessment of pollen sources on a regional scale. Poor geographical and temporal coverage of the pollen observation data can limit this approach in some regions.

Skjøth et al. (2008a) used the bottom-up approach to create the Tree Species Inventory (TSI) for 39 species (including birch, oak, alder, etc.) redistributed to a 50 km grid to be used in atmospheric dispersion models. The top-down approach has been used to produce ragweed pollen inventories for the Pannonian Basin (Skjøth et al., 2010), France (Thibaudon et al., 2014), Austria (Karrer et al., 2015), and Italy (Bonini et al., 2018). The inventories are based on a combination of pollen observation data, knowledge on ragweed ecology, and detailed land cover information. The last (i.e. Bonini et al., 2018) also includes influence of the Ophraella communa beetle on the plants. Prank et al. (2013) as well as Hamaoui-Laguel et al. (2015) have both used combined occupancy and climatic habitat quality maps for ragweed from an ecological model. The map covering Europe has been introduced to the atmospheric models as input data, run for several years and calibrated using pollen observation data. A birch pollen source map has been derived by Pauling et al. (2012). The methodology includes using a forest inventory combined with local and global land use datasets. The obtained map has been calibrated by observational pollen data and it covers Europe. An urban-scale grass pollen source inventory has been obtained by Skjøth et al. (2013a) using GIS and remote sensing data. A high-resolution (1 km) pollen source inventory based on 12 taxa of trees, grass, and weeds has been produced by McInnes et al. (2017) for the UK to be mainly used in exposure studies as well as other possible applications. Many of these different pollen source maps show substantial variation, exemplified by the birch maps over the UK produced by McInnes et al. (2017) and Pauling et al. (2012) or the European scale by comparing the maps by Pauling et al. (2012) with those by Siljamo et al. (2013). This suggests substantial disagreement or uncertainty in the production of the pollen source maps – the key input data for pollen dispersion models. An additional problem is the large variability of the pollen production per plant in each specific year. Several approaches have been proposed to address this issue (e.g. Masaka and Maguchi, 2001; Ranta et al., 2005; Ritenberga et al., 2018).

Zink et al. (2017) applied an atmospheric dispersion model to compare different ragweed pollen source maps using the COSMO-ART model. The results show that the pollen source map resulting from combining land cover and pollen observation data provides the best model performance. The advantage of including a dispersion model is that, in principle, such models take into account atmospheric transport and its effect on pollen concentrations during a specific year (e.g. Prank et al., 2013). A recent study showed high potential of extended four-dimensional variational data assimilation as a rigorous procedure for obtaining the source maps' calibrations (Sofiev, 2019).

The literature review presented in this section shows there are multiple pollen source maps produced for different pollen types and regions and using different approaches and data. Despite this, they are rarely applied to and intercompared using atmospheric dispersion models. This study aims to evaluate three pollen source maps for the birch taxon using an atmospheric dispersion model and study the effect on the model results by combining these source maps with pollen data. The maps and model are described in the Methods section below.

2 Methods

2.1 Birch pollen source maps

Three birch pollen source maps have been chosen in this study: (1) a map derived by combining land cover data and forest inventory, (2) a map obtained from land cover data and calibrated using model simulations and pollen observations, and (3) a statistical map resulting from analysis of forest inventory and forest plot data (see Fig. A1 in the Appendix section). The maps have been preprocessed for two modelling domains: P15 and T15. P15 covers part of Spain in the south to part of Scandinavia in the north, with Denmark near the centre of the model domain (Fig. 1a). T15 covers Finland, parts of Sweden, northwest Russia, the Baltic states, and Belarus (Fig. 1b).

Figure 1 Location of the selected birch pollen observation sites with boundaries of the modelling domains: P15 (a) and T15 (b).

Map 1 (see Fig. 2a, b) has been obtained by combining three datasets: Global Land Cover Characteristics (GLCC, https://doi.org/10.5066/F7GB230D) version 2 (Belward et al., 1999), European Forest Institute tree cover (EFI; Päivinen et al., 2001), and Tree Species Inventory (TSI; Skjøth et al. (2008a)). GLCC has global coverage and includes geographical distribution of land cover, soil type, and other various properties with 1 km horizontal resolution. EFI covers Europe, has 1 km horizontal resolution, and contains spatial distribution of forest with subclasses, i.e. broadleaved forest and coniferous forest. TSI is presented by 39 types of tree species (birch, oak, alder, etc.) and covers Europe with 50 km horizontal resolution and it is based on national forest inventories and national statistics. The following layers have been extracted from these databases and combined to derive the birch pollen source: land mask, forest, vegetation, urban territories (from GLCC), broadleaved forest (from EFI), and birch forest fraction in broadleaved forest (from TSI). For details of the combining procedure see Kurganskiy et al. (2015) and Kurganskiy (2017).

Figure 2 Birch pollen source maps processed for two different modelling domains P15 (a, c, e) and T15 (b, d, f) with 15 km grid resolution. Panels (a, b) correspond to Map 1, (c, d) to Map 2, and (e, f) to Map 3. The maps are shown as a percentage – the percentage cover of birch in each 15 km×15 km grid square.

Map 2 (see Fig. 2c, d) is a product of two datasets and a calibration using model simulations. The two datasets are EFI tree cover and ECOCLIMAP (Masson et al., 2003). ECOCLIMAP is a global dataset with 1 km resolution and it includes 215 ecosystems. The model simulations have been performed with the System for Integrated modeLing of Atmospheric coMposition (SILAM, http://silam.fmi.fi, last access: 23 January 2020, Sofiev et al., 2015b) and pollen observation data for several years to ensure average unbiased representation of the Seasonal Pollen Integral (SPIn) with the SILAM model, where SPIn is the integral of pollen concentrations over time (Galán et al., 2017). The calibrated map is calculated for a 1 km grid covering the whole of Europe with nearby regions. Details of the map calibration procedure are described by Prank et al. (2013) using ragweed, and its application to birch in the SILAM model can be found in Sofiev (2017).

Map 3 (see Fig. 2e, f) is the result of combining ICP-Level-I forest plot data (http://www.icp-forests.org, last access: 23 January 2020) and National Forest Inventory (NFI) plot data as well as the NFI statistics (Brus et al., 2012). Two techniques were used to derive the map: compositional kriging in areas with NFI plot data and regression modelling outside of the areas. The National Forest Inventory statistics were used to scale the regression model results. The map has 1 km horizontal resolution and covers most of Europe. However, information about birch pollen sources in Russia is missing in Map 3. Therefore, the Russian sources were extracted from Map 1 and combined with Map 3 to fill the gap.

The maps were introduced to Enviro-HIRLAM (Environment – High Resolution Limited Area Model) as input data to simulate birch pollen concentrations for the modelling domains.

2.2 Model description

Enviro-HIRLAM is an online coupled meteorology–chemistry model allowing us to take into account spatial–temporal evolution of atmospheric chemical (Korsholm, 2009; Baklanov et al., 2017) and biological (Kurganskiy, 2017) aerosols simulated inside the meteorological HIRLAM model (Undén et al., 2002). The current version of the Enviro-HIRLAM model is based on HIRLAM v7.2 (https://hirlam.org). The detailed model description of all components can be found in Korsholm (2009), Kurganskiy (2017), and Baklanov et al. (2017) while this article is focused on description of the Enviro-HIRLAM features relevant to pollen. In order to simulate birch pollen concentrations, Enviro-HIRLAM requires external information/data on the spatial pollen source distribution and phenology, i.e. start of birch flowering, pollen production, characteristics of pollen release (emissions) within the given study region, and where the spatial pollen source distribution is represented by the maps described above in Sect. 2.1. Atmospheric transport and dispersion as well as dry and wet depositions are calculated within Enviro-HIRLAM at each time step. In numerical models, the start of birch pollen flowering is often determined by applying simple crop growth models (e.g. McMaster and Wilhelm, 1997); here the growing degree day (GDD) method is implemented in Enviro-HIRLAM. GDD is based on accumulation of 2 m daily mean air temperatures above a given threshold (cut-off temperature) as this is the WMO definition of GDD. It is assumed here that birch flowering starts as soon as the accumulated temperature (GDD) reaches a certain value depending on geographical location. These temperature sum threshold maps have been produced by Sofiev et al. (2013) and have been implemented and used as input data for the GDD parameterization in Enviro-HIRLAM. The data were obtained by the fitting procedure applied using the observed leaf unfolding dates (Siljamo et al., 2008) and modelled ones by the GDD parameterization and ECMWF’s ERA-40 reanalysis air temperature data (Uppala et al., 2005) utilized as input for GDD. It was shown in Kurganskiy (2017) that the temperature sum threshold maps are applicable for use to simulate the birch flowering start in Enviro-HIRLAM, even though another model (ERA-40) had been originally used to obtain the maps. GDD is an integral part of Enviro-HIRLAM, and GDD calculations are performed for the same dates as the model is run over. Birch pollen emission is meteorology-dependent and uses parameterizations proposed by Sofiev et al. (2013). It is based on dimensionless functions correcting the seasonal pollen productivity, here chosen to be 3.7×10⁸ $\mathrm{pollen}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{season}}^{-\mathrm{1}}$ as a default value obtained for SILAM (Sofiev et al., 2013). The main meteorological parameters affecting the emissions are air temperature, relative humidity, wind speed, and accumulated precipitation. The emitted birch pollen particles are subject to atmospheric transport, which is parameterized in the same way as for aerosols using the locally mass-conserving semi-Lagrangian (LMCSL) scheme (Kaas, 2008; Sørensen et al., 2013). The scheme provides proper mass conservation, shape preservation, and multi-tracer efficiency required from a numerical point of view. Dry deposition of birch pollen particles is determined by gravitational settling (Seinfeld and Pandis, 2006). The gravitational settling parameterization is based on calculation of settling velocity according to Stokes law and taking into account density of birch pollen particles (800 kg m⁻³) and an estimated size of 22 µm (e.g. Mäkelä, 1996) with near spherical shape (Sofiev et al., 2006). Wet deposition distinguishes between in-cloud (Stier et al., 2005) and below-cloud (Baklanov and Sørensen, 2001) scavenging, i.e. removal of pollen particles from the atmosphere by precipitation. Both scavenging types are parameterized through the scavenging coefficients used for aerosol particles with r≥10 µm.

2.3 Model configuration

The Enviro-HIRLAM model is used to calculate pollen concentrations in two different modelling domains shown in Fig. 1. Both domains have about 15 km horizontal resolution and 40 hybrid vertical levels. Input meteorological initial and boundary conditions were taken from the ECMWF IFS model (Persson, 2011) with 15 km resolution and 6 h intervals. The birch pollen emission module has the same settings as described in Sofiev et al. (2013) and Siljamo (2013). The Enviro-HIRLAM model was run from 1 March until 15 June 2006 with the ECMWF input data covering the same period. The two model domains have been simulated with three types of calculations: (1) a standard calculation, (2) a calculation that includes a correction for 2 m air temperature (T_2 m) bias between assimilated and simulated T_2 m (termed COR), and (3) a calculation that includes both the correction for temperature and a grid-based scaling factor for pollen emissions (termed SCF). These calculations were performed using each of the selected pollen source maps, hence nine model calculations for each of the two domains. The model results from model simulations 2 and 3 are presented in the Results section. The correction for T_2 m bias is done in two steps. Firstly, the model output from the simulations of type 1 is used to calculate the biases (differences between simulated and assimilated T_2 m) for each assimilation window (i.e. 6 h). Secondly, the calculated biases are introduced in the model simulations of type 2 (COR) where the simulated T_2 m is corrected at each model time step using the average bias between the two nearest/closest assimilation windows. Further details of this procedure can be found in Kurganskiy (2017). The scaling factor for pollen emissions is based on the ratio between simulated and observed seasonal pollen integral (SPIn) (see next section) from the simulations of type 2 (COR). This provides 15 point values that have been interpolated to the entire model domain, thereby providing a grid-based scaling factor for each of the tree source maps. The interpolation procedure takes into account the weighted distance between each observation point and grid cell with a constant radius of influence (1 km in this study). The procedure also ensures the scaling factors are equal to around 1 in the areas located far away from the observation points. This is especially visible in P15 domain, e.g. in France (Fig. 3a, c, e). Further details of the interpolation procedure can be found in Kurganskiy (2017). The grid-based scaling factor has then been extracted for each model domain, thus providing six grid-based scaling factors (Fig. 3) applied for the SCF simulations for all maps. The principal difference between the runs is that the COR runs include a constant value for birch pollen productivity (3.7×10⁸ $\mathrm{pollen}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{season}}^{-\mathrm{1}}$) in each grid cell, whereas the productivity is calibrated using the grid-based scaling factor (obtained from the COR runs and shown in Fig. 3) in the SCF runs. The model results were compared with the pollen observation data at 12 measurement sites located in Finland, Denmark, and Russia (see Fig. 1) for all model runs.

Figure 3 SPIn ratio used to calibrate pollen productivity in the SCF runs for two different domains P15 (a, c, e) and T15 (b, d, f). Panels (a, b) correspond to Map 1, (c, d) to Map 2, and (e, f) to Map 3.

2.4 Used observations and statistical evaluation of model results

Daily concentrations of birch pollen were used from 12 sites (Fig. 1a, b), and furthermore annual SPIn values were obtained from three additional sites in Lithuania, published by Veriankaite (2010), in order to strengthen the calculation of the grid-based scaling factor (Sect. 2.3). The observational data were obtained using 7 d volumetric samplers (Hirst, 1952) and analysed according to standard methods in aerobiology (Galán et al., 2014). Station names and coordinates are presented in the Appendix section along with statistical summaries for each station (see Tables A1–A4), while overall results are presented in the Results section. The statistical evaluation includes standard and threshold-based statistical approaches, e.g. previously used by Siljamo et al. (2013) and Sofiev et al. (2015a) and calculated in this paper for each map (1, 2, 3) and run (COR, SCF) type. The start and end of the birch pollen season are calculated using the retrospective method according to Nilsson and Persson (1981), here defined as the dates when the total sums of birch pollen concentrations (SPIn) reach 5 % and 95 %, correspondingly. The criteria are chosen to provide filtration of LRT (long-range-transport) episodes which have a significant impact in northern Europe, especially in the beginning of the birch pollen season (Ranta et al., 2006; Skjøth et al., 2007; Mahura et al., 2007; Veriankaite, 2010). The standard statistical approach includes calculations of correlation coefficient (R), coefficient of determination (R²), mean bias (MB), and root-mean-square error (RMSE), and the evaluation primarily concerns correlation and RMSE. Two additional metrics, normalized mean bias factor (NMBF) and normalized mean absolute error factor (NMAEF), are also calculated according to Yu et al. (2006) since RMSE is highly sensitive to outliers. The R² values are shown on the time series plots for each station, map, and run type as well as on the global scatter plots. The standard statistical metrics have been calculated for the “main” pollen season (30 April–15 June 2006) due to non-ergodicity (i.e. non-similarity in space and time) and non-stationarity of the birch pollen time series (Sofiev et al., 2015a). The threshold-based statistical approach (Siljamo et al., 2013) comprises calculations of the model accuracy (MA), probability of detection (POD), false alarm ratio (FAR), probability of false detection (POFD), and odds ratio (OR) using the threshold value for birch pollen concentrations C_th=50 pollen m⁻³. The threshold value C_th is chosen since most of the population with pollen allergies might start suffering from allergic reactions when daily mean birch pollen concentration C≥C_th in the air (Jantunen et al., 2012). Additionally, the modelled and observed birch pollen concentrations have been sorted into classes in order to estimate the number of cases for low (1–10 pollen m⁻³), moderate (10–100 pollen m⁻³), high (100–1000 pollen m⁻³), and very high (>1000 pollen m⁻³) birch pollen concentrations. The sorted data plotted with corresponding standard errors are shown in the Results section. Hit rates (HRs) are calculated for the classes, and their results are shown in tables in the Results section. The results of the statistical significance test (chi-squared test) performed for the HR values are also discussed. For summary of the thresholds, classes, and metrics used to calculate the threshold-based statistics, the reader is referred to Tables A5–A6 in the Appendix section.

3 Results

The simulated time series of birch pollen concentrations have been compared to the selected pollen observation sites by determining the start and end of the birch pollen season and applying the standard statistics as well as threshold-based statistics described in the Methods section. The modelled and observed time series are also shown in Figs. 4–5 for COR and SCF runs, respectively. The start of the birch pollen season is simulated quite well for both runs and all maps. Station-wise the bias does not exceed 2–3 d at all stations and indicates 2 d too early flowering on average for all maps. Simulation of the last flowering day turns out to be the most challenging for Enviro-HIRLAM. The bias shows too late flowering at almost all stations. On average, it varies from 15 d (Map 3) to 9 d (Map 2) according to the COR run. Rescaling of the pollen productivity (SCF run) reduced the bias by 2–3 d for each map. The largest biases for the last flowering day are found at five Finnish sites (Helsinki, Joutseno, Kangasala, Vaasa, and Kuopio) and two Russian sites (Moscow and Smolensk), and they are related to overestimation of the observed birch pollen concentrations at the end of the season for both COR and SCF runs (see Figs. 4 and 5). It can be addressed to large variations in the year-to-year pollen production in that region which, in turn, introduce a large bias in the model simulations.

Figure 4 Time series of the modelled and observed birch pollen daily mean concentrations for the selected pollen observation sites. The modelled concentrations correspond to each simulated map according to the COR run. A base-10 log scale is used for the y axis.

Figure 5 Time series of the modelled and observed birch pollen daily mean concentrations for the selected pollen observation sites. The modelled concentrations correspond to each simulated map according to the SCF run. A base-10 log scale is used for the y axis.

Generally, the time series of the modelled birch pollen concentrations (Figs. 4 and  5) show good agreement with observations for all maps at most of the selected stations. The COR run indicates that the model behaves similarly with respect to describing variance in observations (see R² values in Fig. 4), except for some stations (e.g. Moscow, Helsinki) where Map 2 provides higher R² values in comparison with Map 1 and Map 3. The SCF run did not affect the R² values significantly, but it reduced the mean bias (MB, NMBF) at the stations (see Tables A1–A4). The highest R² values are found at two Danish sites (Copenhagen and Viborg) where the model explains more than 70 % of variance found in observations for all considered maps and model runs. It should be noted here that the R² values decrease slightly at the Danish sites according to the SCF runs. The reason for this is that the simulated concentrations in the SCF runs are increased almost every single day, which has the effect of reducing the underestimation found in the central season (e.g. Fig. 5a) and increasing LRT, most likely originating from sources found in Germany and Poland at the beginning of the season and from Scandinavia at the end of the season. For more detailed statistical summaries at the stations, the reader is referred to tables in the Appendix section.

The global standard statistical metrics for all stations are presented in Table 1. According to the COR run the model performs slightly better with Map 2. It is supported by a higher correlation coefficient (R = 0.59) and lower mean bias (MB =69.08 pollen m⁻³, NMBF =0.09) in comparison with two other maps. One can assume that this result should be expected since pollen observations are already included in the methodology used to obtain Map 2. Rescaling of the pollen productivity (SCF run) provides improvement of the Enviro-HIRLAM model performance by increasing the R values (R=0.72 for Maps 1 and 2 and R=0.68 for Map 3) and decreasing the MB, NMBF, and NMAEF values. However, RMSE is large (i.e. >1000 pollen m⁻³) for all maps and both model runs. It can be explained by sensitivity of RMSE to large differences between modelled and observed values. The large bias for the last flowering day simulation contributes significantly to the RMSE values.

Table 1 Global “standard” statistical metrics obtained by comparing the modelled ($\stackrel{\mathrm{‾}}{C_{\mathrm{m}}}$) and observed ($\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}$) daily mean birch pollen concentrations. The metrics were calculated using all available station data for both simulated domains, over the 30 April–15 June 2006 period. $\stackrel{\mathrm{‾}}{C_{\mathrm{m}}}$, $\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}$, MB, and RMSE are given (pollen m⁻³). The results were statistically significant for all maps (p value <0.01 according to a Student's t test).

The global scatter plots (Fig. 6) also indicate improvements in the modelled–observed concentration agreement expressed as an increase in the R² for the SCF model setup. As seen, the coefficients of determination do not exceed 0.34 (Map 2) for all maps in the COR run. The SCF run provides R² values around 0.5 for all maps. This means that the model can explain about 50 % of variance found in observations after the rescaling using any of the maps.

Figure 6 Scatter plots for birch pollen concentrations when comparing COR (a–c) and SCF (d–f) model runs with observations for Map 1 (a, d), Map 2 (b, e), and Map 3 (c, f).

The threshold-based statistical metrics calculated for each map and model run are shown in Table 2. The results indicate similar behaviour of the model when comparing performance of the maps. The values of the statistical metrics are slightly different between the maps, but the difference does not exceed 7 %. The model accuracy (MA) – the fraction of correct simulations – is higher than 70 % for all model simulations. The POD values show that the Enviro-HIRLAM model simulates the high concentrations (>50 pollen m⁻³) correctly in more than 70 % of cases. The fraction of incorrect simulations for the high concentrations (FAR values) is less than 30 % for the model setups. The POFD values show that the fraction of the cases when concentrations are simulated as high but observed as low (i.e. <50 pollen m⁻³) is more than 30 % for the COR runs, and it decreases below 30 % after rescaling of the pollen productivity (SCF run). As expected the odds ratio (OD) values increase in the SCF run for all maps. The modelled and observed birch pollen concentrations have been grouped in order to estimate the number of cases for low, moderate, high, and very high birch pollen concentrations (see Fig. 7). As can be seen, the Enviro-HIRLAM model simulates the low and very high pollen concentrations quite well with a difference less than ±5 % between the modelled and observed number of cases for both COR (Fig. 7a) and SCF (Fig. 7b) runs and all maps. The largest differences between the modelled and observed number of cases are found for the moderate and high pollen concentrations, but the differences are within the error bars (see Fig. 7). Map 2 underestimates the number of cases for the moderate concentrations by 10 %, whereas Map 3 overestimates the high ones by more than 10 % (COR run). The SCF run slightly improves the results for the moderate and high number of cases.

Table 2 The threshold-based statistical metrics obtained by comparing model simulations with observed daily mean concentrations. MA, POD, FAR, and POFD are given as a percentage.

Figure 7 Histogram of the observed and modelled birch pollen concentrations sorted into low (1–10 pollen m⁻³), moderate (10–100) pollen m⁻³, high (100–1000 pollen m⁻³), and very high (>1000 pollen m⁻³) classes. Panel (a) corresponds to COR and (b) to the SCF model simulations.

The HR calculations (see Tables 3–4) show that simulations of the moderate and high birch pollen concentrations are the most challenging for all model setups. The SCF run (Table 4) indicates improvements of the model score for the moderate and high birch pollen concentrations. However, the hit rates still do not exceed 50 % for all maps. This feature is also visible on the scatter plots in Fig. 6. The HR values for the low pollen concentrations are quite high and exceed 90 % for COR and 80 % for SCF runs. The very high pollen concentrations are simulated slightly better by using Map 3 with HR = 70 %. However, according to the chi-square test the differences between the HR values for different maps are not statistically significant for both COR and SCF runs.

Table 3 Hit rates of the Enviro-HIRLAM birch pollen simulations (COR run) for the classes
zero or low, moderate, high, and very high. The classes are based on daily mean pollen concentrations (pollen m⁻³).

Table 4 Hit rates of the Enviro-HIRLAM birch pollen simulations (SCF run) for the classes zero or low, moderate, high, and very high. The classes are based on daily mean pollen concentrations (pollen m⁻³).

4 Discussion

Analysis of the start of birch pollen season showed good model performance for all maps with the bias not exceeding 2 d on average. Such error is acceptable by state-of-the-art atmospheric dispersion models (e.g. Sofiev et al., 2015a). The analysis did not reveal significant dependency of the start of the birch pollen season on the underlying pollen source map and/or pollen productivity. This could be explained by the fact that the start of the birch pollen season is calculated using meteorological parameters (i.e. temperate sum thresholds) employed in GDD parameterization in the model. However, long-distance transport of birch pollen has often been observed in that region (e.g. Hjelmroos, 1992; Mahura et al., 2007), which can systematically affect the start of the pollen season (Skjøth et al., 2007). It could therefore be assumed that changing the source term in remote regions such as Poland could affect the start of the calculated pollen season in Denmark and Sweden. This phenomenon was not observed, which may be explained by the fact that the scaling factor for Poland was near unity for all three maps, while the scaling factor for Copenhagen and the close surroundings was substantially higher. This observation therefore complements the suggestion from other studies that long-distance transport of pollen is an episodic phenomenon (e.g. Smith et al., 2008; Fernández-Rodríguez et al., 2014) and that the overall contribution to locally observed birch pollen concentration is due to sources found within the region (Sofiev, 2017). Add to this that urban areas have previously been identified as areas that may enhance pollen production of flowering plants (Ziska et al., 2003). This may be particularly relevant for Copenhagen as a substantial fraction of the birch pollen observed is expected to originate from a source found in the city itself (Skjøth et al., 2008b). The last flowering day simulation was the most challenging for Enviro-HIRLAM. Map 2 showed the lowest bias estimated as 6 d (delay of end of flowering) on average for the stations according to the SCF runs. However, such bias is a known feature for state-of-the-art atmospheric models (e.g. Sofiev et al., 2015a). Introducing more pollen observation stations and performing more iterations for rescaling of pollen emissions (i.e. repeating SCF runs several times) would strengthen the current approach, and it could potentially improve the modelling of the last flowering day. Station-wise the Finnish and Russian sites had bias ≥ 2 weeks in all model runs. One of the possible reasons for this is that the large pollen emission area found in Russia, north of Moscow and Smolensk, is overestimated for this particular year. The source is found in all three source maps and the overestimation is further supported by the fact that the model simulates large amounts of birch pollen concentration very late in the season at Moscow and Smolensk (Fig. 4g and j), most likely originating from more northerly regions. This area will occasionally also affect the Finnish stations as it is well known that this region receives LRT from Russia (e.g. Sofiev et al., 2006). The study therefore complements a multi-model ensemble that demonstrated generally good results for Europe and also highlighted that large variations in the year-to-year pollen production in that region can introduce a large bias in the model simulations (Sofiev et al., 2015a). Large annual variations in the birch pollen integral are well known for a number of European countries (Latałowa et al., 2002; Dahl et al., 2013; Piotrowska and Kaszewski, 2011), as these annual fluctuations are generally not accounted for in the atmospheric models (e.g. Ritenberga et al., 2018).

The standard statistical metrics indicated a slightly better performance of the Enviro-HIRLAM with Map 2 on average for the stations according to the COR run. However, even for this map, the model explained less than 35 % of variance in observations. It is therefore plausible to attribute the remaining part of the disagreement to small-scale patterns of the birch habitation, which are neither captured by existing calibration methods nor satisfactorily reproduced by the land cover datasets. In this study, a scaling of the pollen productivity improved the model results significantly for all source terms and provided near-identical results for Map 1 and Map 2. The small-scale pattern in birch habitation will be difficult to capture by the land cover datasets, which will have the effect that a small but diffuse emission source found in many locations is not present, causing the model to underestimate medium concentration levels. Conversely, the large emission source found in parts of Russia will cause the model to overestimate periods when long-distance transport is present. This effect of over- and underestimation will vary from year to year, depending on pollen productivity, which is known to vary from year to year and between regions (Ranta et al., 2005). This suggests that accurate exposure calculations that use dispersion models should preferably use data fusion that combines a detailed inventory-based source term (e.g. Skjøth et al., 2008a; Pekkarinen et al., 2009; Sofiev, 2017) along with observed pollen concentrations covering the period of interest. Other possible sources of error (e.g. timing of pollen release, quality of meteorological parameters, simulation of atmospheric transport and removal processes) could also be taken into account. However, their effect is less important in comparison with the quality of emission (source) maps – the main uncertainty in air pollution and pollen dispersion modelling.

The threshold-based statistics calculated relative to one threshold value (C_th=50 pollen m⁻³) showed a good model agreement with observations independently of the map used. The threshold-based analysis was also carried out for low (1–10 pollen m⁻³), moderate (10–100 pollen m⁻³), high (100–1000 pollen m⁻³), and very high (>1000 pollen m⁻³) birch pollen concentrations. According to the analysis, simulations of the moderate and high birch pollen concentrations appeared to be the most challenging for Enviro-HIRLAM for all maps. The hit rates were less than 50 % for all model setups. This reveals a need for improvement of the rescaling of birch pollen productivity by introducing more observational points, in particular in areas sparse with data or regions with high emissions such as parts of Russia, and/or performing more iterations. The differences of the HR values appeared to be not significant during inter-comparison of the maps. This study concerning birch pollen complements other studies on ragweed (Prank et al., 2013; Hamaoui-Laguel et al., 2015), which also demonstrated the need for recalibration of the source term. However, it has since been shown by Zink et al. (2017) that source terms combining pollen data from several years with detailed land cover data can outperform other approaches minimizing the need for local calibrations. This approach is also limited by the fact that only a few regions have a sufficiently high number of pollen observations (Zink et al., 2017). Furthermore, abrupt annual changes in the underlying pollen integral can be caused by changes in the pollen emission (Bonini et al., 2015, 2018).

The insufficient density of the observational network also became the methodological limitation of the current study: we have used the same data from the same stations and for the same years for both calibration and evaluation. However, we separated them frequency-wise: calibration uses the annual or multi-annual average whereas the evaluation primarily concerned correlation and RMSE. It should be noted that an independent evaluation that does not use pollen data is shown by analysis of the COR runs in the study. Another approach could be to carry out a cross-validation procedure using the so-called leave-one-out procedure as this tests the sensitivity of individual data points. However, such a sensitivity study requires 15 additional model simulations for each map and each modelling domain, hence 90 simulations. The increased computational cost makes this exercise unfeasible, suggesting that other approaches need to be developed.

5 Conclusions

The aim of the study was to evaluate three birch pollen source maps using the atmospheric dispersion model Enviro-HIRLAM. The evaluation has been performed for 12 pollen observation sites located in Denmark, Finland, and Russia. The modelled and observed time series of birch pollen concentrations have been analysed for the start and end of flowering as well as by calculating the standard and threshold-based statistical metrics for two sets of the Enviro-HIRLAM model runs: without and with calibration using pollen observations.

The analysis did not reveal significant dependency of the start–end of the birch pollen season on the underlying pollen source map. The statistical analysis showed a good model agreement with the observed birch pollen concentrations in the region studied. The model can explain up to 50 % of variance found in observations after the calibration with all maps. It was shown that the remaining part of the disagreement should be attributed to small-scale patterns of the birch habitation, which are neither captured by existing calibration methods nor satisfactorily reproduced by the land cover datasets. The analysis also revealed that the insufficient density of the observational network was the methodological limitation of the study.

Generally, it was shown that calibration of the maps using pollen observations covering the investigation year significantly improved the model performance for all three maps. The findings also indicate the large sensitivity of the model results to the source maps and agree well with other studies on birch showing that pollen or hybrid-based source maps provide the best model performance. This study highlights the importance of including pollen data in the production of source maps for pollen dispersion modelling and for exposure studies.

Appendix A

Figure A1 Flowcharts outlining the data and methods used to obtain Map 1 (a), Map 2 (b), and Map 3 (c).

Table A1 Birch pollen statistical metrics for the selected observation sites: COR run. $\stackrel{\mathrm{‾}}{C_{\mathrm{m}}}$, $\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}$, MB, and RMSE are in pollen per cubic metre; p value <0.05.

Table A2 Birch pollen statistical metrics for the selected observation sites: COR run. $\stackrel{\mathrm{‾}}{C_{\mathrm{m}}}$, $\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}$, MB, and RMSE are in pollen per cubic metre; p value <0.05.

Table A3 Birch pollen statistical metrics for the selected observation sites: SCF run. $\stackrel{\mathrm{‾}}{C_{\mathrm{m}}}$, $\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}$, MB, and RMSE are in pollen per cubic metre; p value <0.05.

Table A4 Birch pollen statistical metrics for the selected observation sites: SCF run. $\stackrel{\mathrm{‾}}{C_{\mathrm{m}}}$, $\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}$, MB, and RMSE are in pollen per cubic metre; p value <0.05.

Table A5 Thresholds and classes used to calculate the threshold-based statistical metrics for 12 selected pollen observational sites.

Table A6 Definitions of the threshold-based statistical metrics used to compare the modelled and observed daily mean birch pollen concentrations. See Siljamo et al. (2013) for more details.