# European pollen reanalysis, 1980-2022, for alder, birch, and olive

Mikhail Sofiev, Julia Palamarchuk, Rostislav Kouznetsov, Tamuna

Abramidze, Beverley Adams-Groom, Célia M. Antunes, Arturo H. Ariño,

Maximilian Bastl, Jordina Belmonte, Uwe E. Berger, et al.

► To cite this version:

Mikhail Sofiev, Julia Palamarchuk, Rostislav Kouznetsov, Tamuna Abramidze, Beverley Adams-Groom, et al.. European pollen reanalysis, 1980-2022, for alder, birch, and olive. Scientific Data, 2024, Scientific Data, 11, pp.1082. 10.1038/s41597-024-03686-2 . hal-04780449

HAL Id: hal-04780449

https://lilloa.hal.science/hal-04780449v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

DATA DESCRIPTOR

# European pollen reanalysis, 1980–2022, for alder, birch, and olive

Mikhail Sofiev et al. #

The dataset presents a 43 year-long reanalysis of pollen seasons for three major allergenic genera of trees in Europe: alder (Alnus), birch (Betula), and olive (Olea). Driven by the meteorological reanalysis ERA5, the atmospheric composition model SILAM predicted the flowering period and calculated the Europe-wide dispersion pattern of pollen for the years 1980–2022. The model applied an extended 4-dimensional variational data assimilation of in-situ observations of aerobiological networks in 34 European countries to reproduce the inter-annual variability and trends of pollen production and distribution. The control variable of the assimilation procedure was the total pollen release during each flowering season, implemented as an annual correction factor to the mean pollen production. The dataset was designed as an input to studies on climate-induced and anthropogenically driven changes in the European vegetation, biodiversity monitoring, bioaerosol modelling and assessment, as well as, in combination with intra-seasonal observations, for health-related applications.

## Background & Summary

Airborne pollen released by anemophilous plants during their flowering season can cause significant allergic manifestations impairing public health, especially in combination with other air pollutants and/or weather phenomena¹⁻⁵. The health impact of pollen has exhibited an upward trend, with only a very low fraction of the population being concerned in the 1960s. In the late-1990s – early-2000s the prevalence of allergic rhinitis (AR) in Europe was estimated at ~15%⁶; with a range of 11–32% deduced by a study of 10 European countries⁷. In the 2010-s, the European Community Respiratory Health Survey established the prevalence of AR to be from 4% to 32%⁸. Recently, a thorough review indicated the range of 1% - 40% and reported a gradual increase of prevalence in already heavily affected countries⁹. The rise of AR cases coincided with an onset of asthma epidemics. The topic got particular attention due to a recent outbreak of thunderstorm-related asthma attacks¹⁰.

The increase of the prevalence of pollinosis and other types of allergy is commonly attributed to a “western lifestyle”, following the hygiene hypothesis suggested in the 1980s¹¹. However, since then it has been clarified that a significant contribution stems from other factors¹²,¹³. From the environmental standpoint, land use and land cover transformations, on-going climate change, rising CO₂ levels and temperature can contribute to changes in exposure to pollen and to subsequent changes in allergy prevalence. Therefore, a retrospective assessment of concentrations of allergenic pollen and their trends is important for understanding the epidemiology of pollen allergy and constructing long-term exposure-response models¹⁴,¹⁵.

Apart from the public health-related motivation, information on biological particles in the air can shed light on large-scale changes in biodiversity, ecosystem services, species migration, habitat degradation, planting preferences, etc¹⁶. Pollen seasonal abundance could serve as one of potential markers of the ecosystem state and composition.

Today, long-term retrospective assessments of pollen concentrations are based on a limited number of long *in-situ* observation time series¹⁷⁻¹⁹, whose representativeness for other regions is often unclear. The only pollen-related modelling study covering several decades was made for birch and grass without assimilating the observations²⁰. It aimed at evaluation of meteorology-driven trends and variability, thus leaving the land use and plant productivity changes out of scope.

Model-based assimilation of observed concentrations of atmospheric tracers can potentially improve both air quality forecasts and retrospective analysis²¹⁻²⁷. However, practical large-scale and operational applications of data assimilation (DA) technique for air quality are comparatively rare, partly because the classical DA brings limited improvement of model performance, whereas the extended formulations are complicated

<footnote>*A full list of authors and their affiliations appears at the end of the paper.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1 Scheme of the European Pollen Reanalysis.</figcaption></figure>

and costly<sup>21-23,28</sup>. The extended methods are regularly used only for CO<sub>2</sub> source apportionment, inversion, and attribution<sup>21,29-32</sup>. The CO<sub>2</sub> emission monitoring, initially mostly based on an inversion of transport matrices<sup>32</sup>, now involves sophisticated variational DA methods, Ensemble Kalman Filter and Smoother, Bayesian inversion, etc<sup>33-35</sup>. An important feature of most CO<sub>2</sub>-related studies is a long temporal scale (window) of assimilation. It reduces the impact of random fluctuations of an otherwise quite homogeneous field of CO<sub>2</sub> concentration, thus increasing sensitivity of the method to small but systematic sources and sinks.

The goals of the current study are: (i) to perform model-based reanalysis of pollen dispersion for alder (1980–2022), birch (1980–2022), and olive (1985–2022) over Europe, (ii) to evaluate the extended 4-dimensional variational data assimilation (4D-VAR) technique with a long assimilation window for reproducing the pollen season strength, and, finally, (iii) to produce a publicly available long-term Europe-wide pollen reanalysis dataset. The main components and procedures used for generation of the pollen reanalysis are summarized in Fig. 1.

## Methods

**Pollen modelling and data assimilation technology.** Modern pollen dispersion models demonstrate quite good skills in representing timing of pollen seasons of trees, as well as variations of pollen concentrations during the season<sup>36-40</sup>. Comparatively simple but efficient parameterizations for pollen production and release were also described long ago<sup>39-42</sup>. The absolute levels of pollen concentration are reproduced much worse. To-date, there is no Europe-wide model for predicting the intensity of the pollen season (e.g., Seasonal Pollen Integral, SPIn<sup>43</sup>). Practically applicable models have been made only for Northern Europe<sup>44</sup> and Bavaria<sup>45</sup>. Improving the SPIn representation via DA is therefore a high priority task with an important practical outcome.

To-date, there are few practical examples of data assimilation applied to pollen dispersion, whether in operational forecasting or in retrospective analysis. Individual experiments cover a range of approaches, such as pre-processing the model setup using real-time pollen observations, model-measurement fusion and other post-processing, and up to full-scale pollen DA studies<sup>28,38,46</sup>. The first demonstration of 4D-VAR for pollen source inversion was performed with the SILAM model (System for Integrated modeLling of Atmospheric coMposition, http://silam.fmi.fi<sup>47</sup>) within the MACC project (Monitoring Atmospheric Composition and Climate). A more detailed study by Sofiev<sup>28</sup>, hereinafter referred to as S19, expanded the methodology of Vira & Sofiev<sup>27</sup>, hereinafter VS12, and compared efficiency of the extended 4D-VAR and classical DA methods.

**Problem formulation.** Addressing arguably the most-uncertain parameter of pollen dispersion models, the reanalysis took the total Seasonal Pollen Integral (SPIn<sup>43</sup>), a time-wise integral of concentration [pollen day m<sup>-3</sup>], as the target parameter for assimilation. For many tree taxa, the total amount of pollen released during a season is related to the previous-season flowering intensity and to conditions during the inflorescence buds formation, which typically occurs during the spring-summer months. The conditions of the current season can only change the timing of the release and, in rare cases, reduce its intensity by destroying the tree inflorescences by a late frost or long intense rain. Under this assumption, a convenient target parameter of the SILAM pollen module is a map of total pollen production E(i, j, yr): the number of pollen grains released during the whole season per unit area of surface covered by the corresponding tree taxon. Here, i, j are grid indices in the corresponding maps, and yr is year. This parameter was suggested by the S19 study<sup>28</sup>.

**Input data.** The input to the reanalysis includes meteorological data, land use information, plant distribution areas, and pollen observations.

**Meteorological input.** The pollen reanalysis is driven by the meteorological reanalysis ERA5 of the European Centre for Medium-Range Weather Forecasting, ECMWF (http://www.ecmwf.int)<sup>48</sup>. The dataset has been downloaded from the ECMWF archive in a regular lon-lat grid with resolution of 0.25° × 0.25°, which corresponds to the high-resolution products available in the Climate Data Store, CDS (https://cds.climate.copernicus.eu)<sup>48</sup>. The vertical structure of ERA5 consists of 137 hybrid levels, of which SILAM used levels 61–137 (from ~104 hPa down to the surface). The temporal resolution of ERA5 is one hour.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Emission distribution maps of (a) alder, (b) birch, and (c) olive (map from 2020) used in the reanalysis. Unit: see colour bars. Each map shows the multi-annual mean release of pollen grains of the corresponding type per 1 m² of the grid cell area.</figcaption></figure>

Inconsistencies in wind fields due to reprojection from the reduced-Gaussian grid are handled with the SILAM meteorological preprocessor⁴⁹, which also produced additional quantities characterising the atmospheric boundary layer⁵⁰.

*Distribution maps of alder, birch, and olive.* There is no homogenised dataset quantifying the distributions of alder, birch, and olive. Even the most-comprehensive maps of the European Forest Institute (EFI, https://efi.int/knowledge/maps/treespecies, visited 9 Dec 2023), European Atlas of Forest Tree species (https://forest.jrc.ec.europa.eu/en/european-atlas/, visited 9 Dec 2023)⁵¹, or the COoRdinate Information on the Environment (CORINE, https://land.copernicus.eu/pan-european/corine-land-cover, visited 13 Mar 2024) do not cover the Russian & Belarusian territory and tend to contain artefacts in areas with poor local data. Birch and alder are not distinguishable in satellite images, where they are just labeled as a mixture of broadleaf deciduous vegetation. Therefore, the distribution of birch was compiled from several sources following the procedure set by Sofiev *et al.*⁴². Areas with good-quality EFI maps were used as-is and extrapolated latitude-wise to poorly described regions using satellite-observed broadleaf forest maps as proxies. The resulting map had a 0.01° × 0.01° resolution. For alder, the tree-specific information is even scarcer than for birch, therefore a generic land use dataset ECOCLIMAP⁵² was used with a fraction of appropriate land use types assigned to alder.

For olive, a both naturally occurring and cultivated tree, with plantations recognizable from high-resolution satellite images, the CORINE database served as the best proxy. This inventory also has a time dimension: the first map of olive plantations was made in 1990 with collection period covering 1980s. After that, the inventory was updated every 6–10 years, thus allowing for a time-resolving source term. However, CORINE covers only a part of Europe and has no data about Africa. Therefore, Africa had to be filled-in from the ECOCLIMAP general-forest maps, with fractions of olive trees roughly determined by visual inspection of satellite images. Last but not least, olive is a popular ornamental plant in Southern Europe (e.g., https://www.gardenista.com/posts/simple-landscaping-ideas-10-genius-gardens-with-an-olive-tree, visited 6 Jun 2024), which was accounted for by assigning a small fraction (<1%) of urban and suburban land use type to be a source of olive pollen. The urban areas were also taken from CORINE, thus accounting for their evolution in time.

For all species, the final step was a multi-annual calibration of the maps as described by Prank *et al.*⁵³ aiming at unbiased representation of the SPIn, all over the modelling domain. The result is a distribution of mean annual pollen emission per unit area (Fig. 2).

The difficulties in construction of the distribution maps and the related uncertainties resulted in a forced decision to keep the distribution maps of alder and birch constant throughout the reanalysis, thus relying on DA to reproduce both the short-term and long-term changes in the plant abundance and pollen production.

*Pollen observations.* The standard device for pollen observations in Europe through decades has been the Hirst-type pollen trap⁵⁴, which is one of the longest-used devices in atmospheric composition observations. It is an impactor: the particles are sucked into the inlet orifice at ~10 l min⁻¹ flow rate and, due to their inertia, hit a sticky tape located 0.7 mm behind the orifice. The flow rate and the tape position are selected so that >50% of particles of approximately 5 µm of aerodynamic diameter hit the tape and get stuck there. This 5 µm size is considered to be a cut-off of the trap: smaller particles avoid colliding with the tape, whereas practically all particles notably larger than 5 µm get captured. The tape is fixed on a drum rotating with a speed of one revolution per week, which corresponds to about two mm per hour of the tape movement. Upon completion of a weekly cycle, the tape is cut into daily segments and mounted on microscope slides and then particles stuck on the tape are manually counted using a light microscope. The device is simple in construction and affordable, which made it a de-facto standard in aerobiology for over 70 years. The observational procedures and good practices have been homogenised⁴³,⁵⁵,⁵⁶ and recently brought to the European standard CEN 16868⁵⁷.

A limitation of the Hirst-type devices is the quite tedious manual microscopic analysis of slides collected usually during the previous week, or even earlier. A typical resolution of the produced data is one day, but at a price of additional counting time, 2-hourly data can be generated as well. A higher temporal resolution is not possible without altering the drum rotation speed: the width of the air jet exposing the moving tape with particles approximately corresponds to two hours of the tape motion. The device also has other inherent uncertainties and, sometimes, reliability issues⁵⁸,⁵⁹.
<figure><img src="image_5.png" /><figcaption>Fig. 3 The number of stations used for assimilation and validation for each species over the reanalysis period.</figcaption></figure>

In Europe, most observational groups participate in a voluntary European Aeroallergen Network (EAN), which maintains a common database of observations since 1974 (https://ean.polleninfo.eu/Ean, visited 30 Nov 2023). This reanalysis brought together teams from 34 European countries, thus representing the largest joint effort of the European aerobiological community to date. A full list of participating networks and organisations is provided in the Acknowledgement section. The number of stations available for the reanalysis varied from barely a dozen of sites in the 1970s up to over 300 in 2020 s (Fig. 3), i.e., the number of monitoring sites has increased almost 30-fold compared to the earlier years.

A limitation of the existing observations is the network sparsity in several Eastern European countries and Russia. In addition, the Pollen Information Dienst (PID) network in Germany decided not to participate in the reanalysis. Therefore, in the currently used dataset Germany is represented by the Electronic Pollen Information Network ePIN⁶⁰, which has a comparatively short history of observations, i.e., 5 years for most of the stations, the site in Munich with the time series starting in 1989, and a dense network of 24 stations operated in Bavaria in 2015 within a preparatory project towards ePIN. The ePIN network operates the automatic impactor BAA-500⁶¹,⁶², which was among the devices that compared favourably with the Hirst-type impactors within a recent intercomparison campaign⁶³. All other sites included in the reanalysis operated Hirst-type devices throughout the considered period.

**SILAM atmospheric composition model.** The reanalysis was performed with the SILAM model, v.5.9.1. It is an offline global-to-local chemistry transport model developed for evaluating atmospheric composition and air quality⁴⁷, emergency decision support applications⁶⁴, source inversion problems²⁶,²⁷, and analysis of observations⁶⁵,⁶⁶. SILAM is equipped with a variety of source modules, of which this study used the pollen release module³⁸,³⁹,⁶⁷. Pollen is transported in the air as a chemically inert aerosol, whose mixing and deposition followed the standard SILAM descriptions⁵⁰,⁶⁸,⁶⁹.

SILAM AQ predictions have been evaluated in many regional and global studies showing robust performance⁴⁹,⁷⁰⁻⁷⁴. Operational evaluations of its global and regional atmospheric composition and AQ forecasts are available at http://silam.fmi.fi (visited 20 May 2024) for a global comparison against satellite-derived aerosol optical depth, at https://cams2-83.aeroval.met.no/evaluation.php (visited 20 May 2024) for a European evaluation of 6 pollutants against *in-situ* EIONET data, at https://hpfx.collab.science.gc.ca/~svfs000/na-aq-mm-fe/dist/ (visited 20 Sep 2023) for a North-American evaluation of O₃, NO₂, PM₂.₅ against *in-situ* sites in Canada and US, and at https://dust.aemet.es/ (visited 20 Sep 2023) for a North-African, Middle East and South-European evaluation of dust against satellite and *in-situ* observations.

The SILAM bioaerosol source module currently includes twelve types of bioaerosols: pollen of alder, birch, grass, olive, five groups of *Artemisia* genera, ragweed, hazel, and a recent extension to nonclassical bioaerosols, such as aphids and ladybirds³⁸,³⁹,⁵³,⁶⁷,⁷⁵,⁷⁶.

SILAM is an open-code system (https://github.com/fmidev/silam-model, accessed 8 Dec 2023). The reanalysis was performed with the latest operational release of SILAM v.5.9.1 (DOI 10.5281/zenodo.10351493) available from the above repository.
<table><caption>**Table 1.** Parameters of the SILAM pollen source term for trees.</caption><thead><tr><th></th><th>Alder</th><th>Birch</th><th>Olive</th></tr></thead><tbody><tr><td colspan="4"><strong>Pollen development model</strong></td></tr><tr><td>Heat-sum <i>H</i> type</td><td>Hourly mean-T</td><td>Daily mean-T</td><td>Daily mean-T</td></tr><tr><td>Uncertainty of <i>H</i>, [%]</td><td>10</td><td>10</td><td>10</td></tr><tr><td>Start day of heat accumulation</td><td>1 January</td><td>1 March</td><td>1 January</td></tr><tr><td>Cut-off temperature, [°C]</td><td>4</td><td>3.5</td><td>0</td></tr><tr><td>Standard pollen release, <i>N</i><sub>tot</sub> [pollen grains m<sup>-2</sup> yr<sup>-1</sup>]</td><td>10<sup>8</sup></td><td>10<sup>8</sup></td><td>3 × 10<sup>8</sup></td></tr><tr><td>Uncertainty of <i>N</i><sub>tot</sub> [%]</td><td>10</td><td>10</td><td>10</td></tr><tr><td colspan="4"><strong>Pollen release model</strong></td></tr><tr><td>Shortest release time, <i>τ</i>, [hour]</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Low-humidity threshold, <i>q</i><sub>low</sub>, [%]</td><td>50</td><td>50</td><td>50</td></tr><tr><td>High-humidity threshold, <i>q</i><sub>high</sub>, [%]</td><td>80</td><td>90</td><td>80</td></tr><tr><td>Precipitation_threshold, <i>P</i><sub>high</sub> [mm hour<sup>-1</sup>]</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>Wind speed saturation level, [m sec<sup>-1</sup>]</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Wind speed max scaling</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>Emission injection height range, [m]</td><td>1–50</td><td>1–50</td><td>2–50</td></tr><tr><td>Pollen size, [µm]</td><td>22</td><td>22</td><td>28</td></tr><tr><td>Pollen density, [kg m<sup>-3</sup>]</td><td>800</td><td>800</td><td>800</td></tr></tbody></table>

**Tree pollen source term.** All tree pollen sources in SILAM are described by using the same double-threshold heat-sum approach³⁹,⁴¹ with tree-specific distribution maps and coefficients of the governing equations. It distinguishes between two stages of pollen emission: development of ready-to-fly pollen and its release into the air. The development stage is controlled exclusively by accumulated heat, whereas the release is a function of actual meteorological conditions³⁹. The increment of ready-to-fly pollen *p*<sub>rdy</sub> is described as³⁹:

$$ \frac{dp_{rdy}}{dt} = \begin{cases} 0: & H < H_{fs} * (1 - \delta_H) \\ S\phi N_{tot} \frac{T - T_{co}}{\Delta H} h(T - T_{co}) p_{fs} \left(\frac{H}{H_{fs}}\right) p_{fe} (R): & H > H_{fs}, R < 1 + \delta_N \\ 0: & R > 1 + \delta_N \end{cases} \quad (1) $$

Here *S(i,j)* is area [m²] of grid cell (i,j), *φ(i,j)* is fraction of area covered by the corresponding tree crowns, *N<sub>tot</sub>* is amount of pollen released during the full season per unit area fully covered by crowns of the pollen-producing trees [pollen grains] (Fig. 2 shows the product (*φ* *N<sub>tot</sub>)), *T* is temperature [K], either the hourly or the daily mean, depending on the taxon, *T<sub>co</sub>* is cut-off temperature [K], *ΔH* is difference between the heat-sum of the end and the start of the flowering season (*H<sub>fe</sub>*, *H<sub>fs</sub>*), respectively [K day] or [K hour], *h(.)* is Heaviside step function (=0 if its argument is negative and =1 otherwise), *p<sub>fs</sub>* and *p<sub>fe</sub>* are probabilities of a single tree to start or end the flowering, *R* is fraction of *N<sub>tot</sub>* released until the current moment, *δ<sub>H</sub>* is uncertainty of heat-sum threshold, and *p<sub>rdy</sub>* is amount of pollen available for release in the grid cell [pollen grains].

Whenever $p_{rdy} > 0$, the release rate is³⁹:

$$ E(i, j) = p_{rdy}(i, j) \left(1 - e^{-\frac{\Delta t}{\tau}}\right) f_{wind}(U, w^*) f_{thr}(q, q_{low}, q_{high}) f_{thr}(P, P_{low}, P_{high}) \quad (2) $$

Here *E* is emission rate in the grid cell [pollen sec<sup>-1</sup>], *τ* is time constant for the exponential decay of the ready-to-fly pollen under normal conditions [sec], *f<sub>wind</sub>* is correction for wind speed *U* [m sec<sup>-1</sup>] and convective velocity scale *w*<sup>*</sup> [m sec<sup>-1</sup>], *f<sub>thr</sub>* is threshold correction function applied to relative humidity *q* (thresholds *q*<sub>low</sub> and *q*<sub>high</sub>) [%] and precipitation rate *P* (thresholds *P*<sub>low</sub> = 0 and *P*<sub>high</sub>) [mm sec<sup>-1</sup>]. Formulas for the components of Equations (1) and (2) are presented in Sofiev *et al.*³⁹, coefficients for the specific genera of trees are shown in Table 1, and maps of heat sum thresholds for the start and the end of flowering periods are presented in Fig. 4.

**Setup of the 4D-VAR assimilation procedure.** Control variable. Similar to S19, this reanalysis uses the extended 4D-VAR technique²⁷. We denote the target parameter, or a control variable, as *ξ*, and define the model operator *M* mapping this variable to a unique phase-space trajectory of the model state *x* through *x* = *Mξ*. The mapping is defined over some finite time interval referred to as assimilation window. The vector of observations *y* corresponds to *x* via the observation operator *O*: *y* = *O(x)* + *ε*, where *ε* is the observation error. Assuming *ε* to be Gaussian, the maximum likelihood of *ξ* corresponds to a minimum of the cost function:

$$ J(\xi) = \frac{1}{2}(y - Ox)^T \Phi^{-1}(y - Ox) + \frac{1}{2}(\xi - \xi_b)^T B^{-1}(\xi - \xi_b) \quad (3) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4 Heat-sum start/end flowering thresholds for alder, birch and olive. No birch trees were assumed south of 40 N.</figcaption></figure>

Here $\Phi$ and $B$ are the observational and background covariance matrices, respectively. The first term is a penalty for deviation from observations whereas the second one penalizes the deviation from the initial guess of the control variable.

The parameter deciding the total seasonal pollen release, $N_{tot}$, has a climatological value specific for each tree type (Table 1). The dependence of $N_{tot}$ on year and location has been disregarded in the original model³⁹, but it is straightforward to introduce it via a correction map. The control variable $\xi$: of Eq. (3) can be set equal to this map and defined as:

$$N(i, j, yr) = \xi(i, j, yr)N_{tot} \qquad (4)$$

This correction map is unitless and it's a-priori value is $\xi_b = 1$.

There is a physically meaningful restriction, $\xi \ge 0$, which is inconvenient from a practical standpoint because requires conditional minimization procedure. Therefore, in SILAM 4D-VAR, this variable is used in a log-transformed form: $\xi = \exp(\zeta)$, $\zeta_b = 0$. This transformation is monotonic and does not change the features of Eq. (3), but it eliminates the conditional minimization procedure since the range of meaningful values of $\zeta$ is $(-\infty, \infty)$. An evident downside of this transformation is that the optimization can no longer completely switch off the source in a grid cell: the requirement $\xi \ge 0$ is enforced as $\xi > 0$. which may be too limiting in some cases. In practice, however, false emission source can be made negligibly low, so that its impact does not affect the results.

Covariance matrices and measures against over-fitting. We employed the background covariance matrices of VS12 with a few adjustments.

Matrix $B$ incorporates standard deviation of the forecast error and spatial correlation between the grid cells. The former is assumed to be constant, while the latter is set through a correlation distance $\rho$ [m]:

$$corr(\mathbf{r}, \mathbf{r}') = \exp\left(-\frac{\|\mathbf{r}, \mathbf{r}'\|^2}{\rho^2}\right) \qquad (5)$$

Here, $\mathbf{r}, \mathbf{r}'$ are 2-D vectors representing two locations within the 2-D computation grid, $\|\cdot\|$ is the $L_2$ norm of such a vector. Compared to VS12, the following simplifications were made: (i) it is assumed that most of pollen is released and transported within the boundary layer, which is well-mixed, thus eliminating the vertical dimension of the DA problem, (ii) correlation distances latitude- and longitude-wise are assumed to be the same and equal to $\rho$. The first assumption is supported by the near-surface location of pollen sources and the large size of the grains, which keeps them within the boundary layer most of the time. The second one is supported by the VS12 analysis, which showed a barely 10% difference between the meridional and latitudinal correlation distances for sulphur oxides (SOx). For pollen the difference is expected to be even smaller due to the short pollen transport distance and a wide distribution of the sources.

The actual value of the correlation distance $\rho$ depends on the transport features and source distribution, but also on temporal averaging. For this reanalysis, $\rho$ was determined experimentally, by running the assimilation over several years with different $\rho$ values and evaluating the improvement of the model-measurement agreement
for non-assimilated stations. The best improvement was found for $\rho \sim 250$ km, but the sensitivity turned out to be quite low. Compared to $\rho \sim 80$ km for SOx emission determined by VS12, the 250 km correlation distance is reasonable: SOx sources are mostly isolated, e.g., power plants, whereas tree pollen production is controlled by the interplay of synoptic-scale meteorological processes, distribution of sources, and major topographic elements⁴⁴.

Specific absolute values of $B$ and $\Phi$ are not relevant: the location of the minimum of the cost function (3) does not depend on the magnitude of the covariances, but rather on their ratio. Therefore, the observational covariance was assumed constant and diagonal, equal to 6 pollen grains m⁻³, close to the Hirst-device detection limit⁷⁷, and the background correlation $B$, (5), was multiplied with 10, thus allowing an order of magnitude variation of $\xi$ with minor penalty, i.e., $B = 1000 \text{ corr}(r, r')$. Both values are at the low end of the corresponding uncertainties, but their ratio is of the correct order of magnitude. Computations showed that with such values, the observation-related term of the cost function is at least 100 times larger than the background, i.e. the minimization is sensitive to even small (1–10%) improvements of the model-measurement distance but disregards the low level noise present in the term (~0.1%).

As a final measure against over-fitting, the iterative minimization procedure of 4D-VAR was monitored with a truncated-iterations watchdog⁷⁸. After every successful adjoint-and-forward cycle, a so-called L-curve⁷⁹ was calculated to check for signs of over-fitting. This method proved efficient in identifying the optimal iteration, after which further minimization iterations were prone to over-fitting. The result of the data assimilation process was set to be equal the value of $\xi$ at the optimal iteration step, thus truncating the minimization cycle. The method of truncated iterations was used by Vira et al.²⁶ and S19, with positive outcome.

**Assimilation window.** The selection of the time-independent control variable suggests that the assimilation window should cover the whole flowering season of each taxon. In Europe, the season gradually propagates from south-west to north-east following the increase of temperature and the melting of snow. Among the considered taxa, alder is the first one to flower, followed by birch and then olives. Time windows that cover their seasons over Europe in all considered years are:

* Alder: 5 January – 31 May

* Birch: 10 March – 1 July

* Olive: 1 April – 31 July

The simulations started on the first day of the heat sum accumulation (Table 1), but the cost function (3) was calculated and minimised only within the assimilation window.

**Observational data averaging and filtering.** The averaging period of the observations for the time-independent control variable could be as long as the assimilation window. However, the mere existence of the pollen season, with the start and the end, makes the problem non-stationary and non-ergodic, i.e., renders averaging inapplicable. Moreover, season-long averaging would require nearly 100% data availability because any value missing from the observations during the season would affect the SPIn and potentially disqualify the time series. However, some averaging was still possible.

The strongest argument in favour of the long averaging period is that it removes high-frequency variations of concentrations, which are not needed for the SPIn assimilation but penalised by RMSE, i.e., the first term of Eq. (3). To demonstrate this, the Eq. (3) is simplified by assuming the observation operator $O$ to be an identity operator and the observational covariance matrix $\Phi$ to be a scalar equal to 1. Let's split both observations and model predictions into their seasonal mean value and high-frequency normalised anomalies: $x = \bar{x}(1 + \mu)$, $y = \bar{y}(1 + \nu)$, where overbar denotes averaging and $\bar{\mu} = \bar{\nu} = 0$, $\bar{x}^2 = \bar{x}^2$ and $\bar{y}^2 = \bar{y}^2$. Then the RMSE term is transformed to:

$$ J = (y - Ox)^T \Phi^{-1} (y - Ox) = \overline{(y - x)^2} = \overline{(\bar{y}(1 + \nu) - \bar{x}(1 + \mu))^2} \\ = \bar{x}^2(1 + \overline{\mu^2}) + \bar{y}^2(1 + \overline{\nu^2}) + 2\bar{y}\bar{x}(1 + \overline{\mu\nu}) \quad (6) $$

Introducing the control variable of assimilation from the Eq. (4): $x = \xi z$, $\bar{x} = \xi \bar{z}$, $z = \bar{z}(1 + \mu)$, one can find the minimum of the error with regard to $\xi$ analytically. Indeed, by substituting $x \to z$ into the Eq. (6), one obtains a quadratic form with regard to $\xi$, which should be minimised:

$$ J = \xi^2 \bar{z}^2 (1 + \overline{\mu^2}) + \bar{y}^2 (1 + \overline{\nu^2}) - 2\xi \bar{z} \quad \bar{y}(1 + \overline{\mu\nu}) \to \min_{\xi} \quad (7) $$

Its minimum is:

$$ \xi_{opt} = \frac{\bar{y}}{\bar{z}} * \frac{1 + \overline{\mu\nu}}{1 + \overline{\mu^2}} \quad (8) $$

The term $\overline{\mu\nu}$ in the numerator of the ratio (8) is the normalised temporal covariance of observations and model predictions, i.e., the correlation coefficient. Should the time series correlate strongly, their normalised fluctuations would be equal, $\xi_{opt}$ would be just a ratio of the mean values, and assimilation would lead to an unbiased solution. But if correlation is below unity, the scaling will be smaller, and the assimilation will
<figure><img src="image_7.png" /><figcaption>Fig. 5 Histogram of correlation coefficients for 1-day (blue) and 2-day (green) mean observations and model predictions. All stations satisfying the completeness requirements (>30 days of data in the given year) are included.</figcaption></figure>

under-estimate the observed concentrations. The degree of the under-estimation is proportional to the standard deviation of the model output.

As a compromising solution, an averaging period of two days was selected, and only stations reproduced with a temporal correlation coefficient better than 0.3 were used. Averaging over two days resulted in a very significant increase of temporal correlation (Fig. 5), thus both reducing the bias (8) and increasing the number of stations passing the 0.3 inclusion threshold. Also, at least 30 days of valid observations were required for a particular year for the station time series to be included, and at least 5 of them were required to be non-zero.

The accepted stations were split into assimilation (80%) and evaluation (20%) subsets. The split was semi-random: it was required that spatially isolated stations were always used for assimilation. Therefore, in each year and for each taxon, a list of mandatorily assimilated stations was manually prescribed, after which the remaining stations were split randomly. The distributions were subsequently checked and, if all stations in some regions appeared in the evaluation list, the split was rerun. This requirement helped maintaining the acceptable coverage of the assimilated stations in the pre-1990 period when the network consisted of few stations (Fig. 3).

**Bias reduction: calibration of the mean release strength.** To correct the regular low bias of the assimilation output originated from the imperfect model-measurement correlation, after the assimilation, the mean pollen source term strength was calibrated against the whole-period observed SPIn following the procedure of Prank *et al.*⁵³.

The SILAM model was run through the whole period with the assimilated season strength and the total multi-annual SPIn was calculated from the co-located observations and model predictions for each station (SPIn<sub>obs</sub> and SPIn<sub>mdl</sub>, respectively). The computations were made only for stations and years with at least 30 daily observations available (no threshold for model-measurement correlation). The relative error r<sub>s</sub> was then computed:

$$r_s = \frac{SPIn_{mdl}}{SPIn_{obs}} \qquad (9)$$

The ratio r<sub>s</sub> was subsequently inter-/extrapolated over the whole domain using a Radial-Basis Function with a linear kernel and a smoothing parameter equal to 10. The obtained smooth gridded correction was used as a climatological correction factor to the map of seasonal pollen emission. Computations were performed independently for each tree taxa.

**Setup of the SILAM model runs.** Totally, three sets of the SILAM model runs were performed: the first guess run, the data assimilation run, and the final run.

**The first guess (reference) run** was made with the unconstrained model, through the whole period for all species not accounting for any year-to-year variability in pollen production. The run set the reference point for the reanalysis and showed the skill scores to outperform. The setup followed the SILAM operational pollen forecast. The horizontal grid was the same as in the European ensemble of Copernicus Atmosphere Monitoring Service (CAMS): 700 × 420 grid cells, resolution 0.1° × 0.1°, longitude range (25W-45E), and latitude range (30N-72N). The vertical structure consisted of 9 uneven stacked layers, up to 6725 m above the surface: 25 m, 50 m, 100 m, 200 m, 400 m, 750 m, 1200 m, 2000 m, and 2000 m thick. Output included hourly 3D concentrations and 2D dry and wet deposition.

**The DA run** was performed with assimilation of the pollen data as described above generating the set of annual pollen emission correction maps for each year and species. Due to high computational demand of 4D-VAR, the assimilation was performed with a coarser resolution of 0.25° × 0.25° and a vertical with 6 layers with thicknesses of 50 m, 100 m, 400 m, 1000 m, 2000m, and 3000 m. The domain was also reduced to cover the observational network with ~5° margin: horizontal grid of 200 × 168 cells, longitude range (10W-40E), latitude range (30N-72N).

The DA run produced two types of output: the annual emission correction maps for each year and species, and near-surface pollen concentrations. The latter was used to calculate the constant-in-time bias-reducing correction map (9). It did not require additional SILAM computations.
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Mean multi-annual SPIn, [pollen day m⁻³]. The Seasonal Pollen Integral for each pollen type is computed as a sum of daily-mean near-surface concentrations, separately for each year. The obtained annual integrals are averaged over the whole period 1980–2022.</figcaption></figure>

**The final run** used the annual emission intensity correction map from the DA run, extrapolated to the east and linearly downscaled from the DA grid to the source grid, and additionally scaled with the bias-reducing map (9). The rest of the setup was identical to the first guess unconstrained run, thus allowing for a direct comparison.

All simulations used zero lateral boundary conditions and a fully reflective top boundary.

Computations were performed using the SILAM v.5.9.1 (https://doi.org/10.5281/zenodo.10351493) with OMP+MPI parallelization (Open Multi-Processing + Message Passing Interface). Assimilation was set year-and species-wise, each using two 128-CPU (Central Processing Unit) nodes of a supercomputer. The runtime for both the DA and the bias elimination run were between 1 and 3 hours depending on the speed of the 4D-VAR convergence and the simulation grids.

## Data Records

The European Pollen Reanalysis v.1.1 is freely available for the public. The metadata and the link to the main archive⁸⁰ have been registered as https://doi.org/10.57707/fmi-b2share.85841086f9db46b882d750eaa9e42515 and PID: http://hdl.handle.net/11304/4766ae3f-f7cb-4967-b601-9b0479600e98 A direct link to the data is: https://european-pollen-reanalysis.lake.fmi.fi/index.html.

The **input** data that were used for constructing the reanalysis have been described in the Methods section.

The **output** data records of the Pollen Reanalysis v.1.1 are the following (for more details, see the above metadata DOI). The **output horizontal grid** is the same as that of the European ensemble of Copernicus Atmosphere Monitoring Service: 700 × 420 grid cells, Mercator projection, resolution 0.1° × 0.1°, longitude range (25W-45E), and latitude range (30N-72N). The **output vertical** consists of 9 stacked layers, up to 6725 m above the surface and being 25 m, 50 m, 100 m, 200 m, 400 m, 750 m, 1200 m, 2000m, and 2000m thick. Supplementary 2D fields are provided in the **data assimilation horizontal grid**, Mercator projection, resolution 0.25° × 0.25° and 200 × 168 grid cells, longitude range (10W-40E), latitude range (30N-72N).

All fields cover the full reanalysis period of 1980–2022.

The following variables are provided in the dataset:

* 2D near-surface pollen hourly concentrations in the output horizontal grid (example in Fig. 6).

* 3D pollen hourly concentrations in the output horizontal grid and output vertical

* 2D hourly dry and wet deposition fields in the output horizontal grid

* 2D annual pollen productivity correction fields, in the DA grid

* 2D seasonal footprint area from the assimilated stations in the DA grid

All files are in the netCDF4 format, closely following the CF-1.3 convention (https://cfconventions.org, visited 3 Dec 2023), tested for viewing with GrADS v.2.0, Python 3.7 netCDF4 library, and NASA PanoPly netCDF/HDF/GRIB data viewer (https://www.giss.nasa.gov/tools/panoply visited 3 Dec 2023).

## Technical Validation

**Validation of the DA procedure.** Online validation during the assimilation iterations is based on the truncated-iteration and L-curve technology⁷⁸,⁷⁹. It verifies the cost-function reduction and checks for potential over-fitting. A typical example is shown in Fig. 7 for alder assimilation for the year 1990. This diagnostic is primarily used for stopping the iterations not relying on absolute tolerances or other formal criteria of the minimization routines, but on more physically grounded considerations of regularization via truncated iterations⁷⁸,⁷⁹. These are based on the assumption that the large eigenvalues of the gradient matrix will be driving the initial quick reduction of the cost function along the corresponding phase dimension(s), whereas the smaller eigenvalues will start driving the process only after the fast decay is over (iteration 6 in Fig. 7). At that point or soon after it the iterations can be stopped. The actual criterion for stopping the iterations in the current reanalysis was the 20-fold reduction of the initial slope of the L-curve.

The strength of the over-fitting and representativeness of the assimilation and evaluation stations are visible from the comparison of dark and light blue curves from the left-hand panel of Fig. 7. In the provided example, 20× reduction of deviation from the assimilated stations corresponds to 10× reduction of error in the evaluation
<figure><img src="image_9.png" /></figure>

**Fig. 7** Diagnostics of the DA iterations: error reduction and L-curve-based stopping criterion for alder in 1990. Left-hand panel: changes in the squared deviation from observations (first term in Eq. (1), dark blue line), squared deviation from the background (second term in Eq. (1), brown line), norm of the gradient of the full cost function (1) (green line), squared deviation from the evaluation subset of the observations (light blue line), and cosine of an angle between two sequential iteration steps (black line at the lower panel), all as functions of the 4D-Var iterations (iteration 0 corresponds to the background state). The lower panel also shows the iterations used for fitting the analytical approximation of the L-curve (green dots at the bottom) – and those excluded from the fitting due to suspected instability of the optimiser. Right-hand panel: L-curve, a dependence of the squared model deviation from the assimilated observations on the deviation of the control variable from the background (the first and the second terms of the Eq. (1), respectively): magenta broken line is plotted iteration-wise, light blue curve is its smooth analytical approximation with a hyperbolic function. In both panels, the red star shows the stopping criterion determined from the analytical L-curve approximation as 5% threshold of the initial slope; the nearest iteration is then considered as the optimal one. The orange cross denotes the first estimate made during the calculations with a simplified approximation procedure, whereas the violet diamond shows the first iteration when the full cost function reaches a 1.1 level of its overall minimum during the optimization.

stations, with a marginal sign of error growth at the iteration 10. This confirms a good correspondence of the assimilation and evaluation stations and manifests negligible over-fitting in this case. A large distortion of the procedure after the iteration 9 reflects the complexity of the problem and a limited stability of the quasi-Newton optimizer: upon completion of the initial well-defined descend, the optimiser becomes sensitive to small errors and noise in the data. The ultimate task of the truncation procedure is to stop the iterations before it happens.

The summary for all years and tree species of the L-curve diagnostics is provided in Fig. 8. As one can see, the most challenging task was assimilation of the olive pollen data: the improvement of performance was limited or sometimes absent. On the upside, the over-fitting was insignificant: RMSE was not much worse for evaluation stations than for assimilated stations, i.e. the assimilation procedure remained consistent even when the improvement was limited. There are several reasons for that behaviour. Firstly, comparatively few stations observing olive pollen are distributed over a very long but narrow area along the Mediterranean coast. Atmospheric transport between these sites is limited, which complicates the selection of evaluation stations: every station is unique. Secondly, the inter-annual variation of olive pollen production is low due to human management of the agricultural crop: irrigation, fertilization, intensification, etc., i.e. the signal to assimilate is low⁸¹. Nevertheless, for about two thirds of the years, the improvement reached or exceeded ~20%, which is a good outcome for such a complicated case.

Assimilation of alder and birch was more efficient: the RMSE improvement was more than a factor of two for alder (down to 48% of initial value) and over a third for birch for the evaluation stations. The over-fitting was very limited: it exceeded a factor of 2 ± 0.1 for only two early years for birch (small green dots outside the range marked by grey lines) when the network density was not sufficient to constrain the model.

**Inter-annual SPIn variation.** The primary target of the assimilation process was to resolve the inter-annual variation of the SPIn, Eq. (4). The effect of the DA can be described in terms of correlation coefficient of the observed and predicted SPIn (Fig. 9). Comparison of the panels of Fig. 9 reveals three signals: (i) even an
<figure><img src="image_10.png" /></figure>

**Fig. 8** Summary of L-curve diagnostics for the whole reanalysis. X-axis: RMSE for assimilated stations of the optimal 4D-VAR iteration normalised with that of the 0th iteration (unconstrained run), y-axis: RMSE for evaluation stations at the optimal iteration normalised with that of the 0th iteration. On both axes, 1.0 means no improvement due to the data assimilation. Points above the 1:1 line show the years where improvement was higher for assimilated stations than for evaluation ones: large displacement shows an overfit. Two side dashed lines show the slope of 0.5 and 2 with an offset of -0.1 and 0.1, respectively. The size of the dots is the smallest for 1980 and gradually grows towards 2022. Legend also shows mean RMSE improvement (a fraction of RMSE left after DA).

<figure><img src="image_11.png" /></figure>

**Fig. 9** Histograms of correlation coefficient between the observed and predicted/assimilated inter-annual SPIn variations. Three cases are considered: unconstrained reference run, raw output of the DA run, and bias-corrected final run. All years and all stations with >3 years of valid observations are included.

unconstrained run (cyan-coloured bars) shows certain skill in reproducing the inter-annual SPIn variation, (ii) assimilated time series demonstrate much higher correlation with the observed SPIn variability, (iii) some stations remain poorly reproduced in all runs. The first observation reflects the effect of the transport conditions during the season on the SPIn, which appeared to be significant, in line with the earlier hindcast computations²⁰. The second one confirms efficiency of the selected DA procedure. The third one shows the areas of improvement: the most-probable reason for persistently poor representation of some stations is a deficiency of regional/local sources, which lead to representativeness issues for the corresponding stations.

The assimilated SPIn fields also exhibit excellent frequency distributions. The plots in Fig. 10 show practically perfect quantiles for all species.

**Efficiency of the post-DA bias reduction.** As shown in the Methods section, assimilation of the imperfectly correlating data leads to a negative bias in the results. This bias has been corrected by the post-DA recalibration of the emission maps following the procedure of Prank et al.⁵³ (Fig. 1). Its efficiency was verified by computing the ratio of the observed and predicted SPIns over the whole period of 43 years at each station:
<figure><img src="image_12.png" /><figcaption>Fig. 10 Quantile plots for assimilated SPIn. The plot is obtained by an independent sorting of the observations and the corresponding model predictions, thus disregarding their temporal co-location and only accounting for the relation between the distribution functions of the observed and predicted concentrations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11 Ratio of observed and modelled multi-annual SPIn.</figcaption></figure>

$$r_{st} = \frac{\sum x_{st}}{\sum y_{st}} \qquad (10)$$

The summation in Eq. (10) is performed over all valid observed daily data of the station, and the corresponding model predictions. Full co-location, in both space and time, ensured consistency and comparability of the multi-annual SPIn.

Figure 11 shows that the model predictions for all species are in-average unbiased: except for a few outliers, the vast majority of the stations report a ratio of modelled and observed SPIn close to unity. This is in agreement with the quantile plots of Fig. 10.

**Temporal correlation of daily concentrations.** The assimilation of the target variable of the reanalysis, SPIn, has a limited and non-equivocal effect on the day-to-day variations of concentrations within the season. Indeed, changing the regional pollen emission intensity alters the long-range transport patterns, thus affecting both source and receptor locations. As SILAM has been calibrated to an inter-annually constant release rate, this reanalysis opens new possibilities for model improvement. However, in the current study the impact of improved SPIn to intra-seasonal model performance was negligible (Fig. 12). The final fields, with the assimilated SPIn and the corrected bias, show similar skills as the initial run, thus confirming that the short-term features of the season remain essentially unconstrained.

**Areas constrained by the observations.** Due to varying number of stations (Fig. 3), area constrained by the data is different each year. Secondly, the PID network did not participate in the reanalysis, potentially reducing the reliability of the results in northern Germany. However, even the comparatively sparse network was efficiently used by the assimilation procedure constraining the emission not only in the vicinity of the stations, but also in the regions covered by footprints of the available stations (Fig. 13 upper row). For instance, even in 1980 when no German networks were operational, pollen release in north-western part of the country was constrained by the data from Belgian stations, whereas its southern part was visible for Austrian sites. Nevertheless, one has to be careful in interpretation of the obtained results during this decade because the well-observed regions in Europe were isolated from each other, and the current assimilation procedure provided only limited information about their pollen levels. With a gradual growth of the network, the coverage improved, so that after ~1990 most of Europe was covered by the station footprints and constrained in the reanalysis (Fig. 13 lower row).
<figure><img src="image_14.png" /><figcaption>**Fig. 12** Histograms of temporal correlation coefficient between observed and predicted 2-daily-mean pollen concentrations, 1980–2022. The coefficient is computed for each station and for each year.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13** Areas observed by the active stations in 1980 (upper row, 1985 for olives) and 2022 (lower row), [relative]. The presented variable is the sum of footprints of all active stations throughout the growth and flowering season 1 January–31 July.</figcaption></figure>

In view of the limited network coverage in some periods and regions, existing model evaluation exercises can be used to estimate uncertainty of the unconstrained model simulations. In particular, birch model was evaluated in a series of studies³⁷,³⁹,⁴²,⁶⁷,⁸²⁻⁸⁵ showing robust performance. An important addition to these studies is a recent model evaluation report of Copernicus Atmospheric Monitoring Service (CAMS)⁸⁶, which evaluated 11 European models, all equipped with the SILAM pollen emission module, against the EAN pollen monitoring data provided to CAMS for model evaluation purposes. The report confirms the small bias and very high correlation of the SILAM unconstrained forecasts for birch and olive in 2022 (Fig. 14).

## Usage Notes

The reanalysis output includes, apart from the near-surface pollen concentrations, the 3-D fields of hourly pollen concentrations, thus allowing for nested high-resolution computations over any part of Europe. Therefore, high-resolution model-based assessments can use the re-analysis as a boundary condition.

Long-term re-analysis is a convenient source for various trend estimates. However, a common problem (frequently under-estimated) is that the trend of availability of the observations can significantly affect the analysis. For pollen re-analysis, the number of monitoring stations varies greatly depending on the specific year (Fig. 3). Therefore, SILAM computations over several regions in 1980s and, to a less extent, 1990s are practically unconstrained. To facilitate the interpretation of the re-analysis and its application to long-term studies, it is accompanied with annual footprints of the available monitoring stations (Fig. 13). As shown by Sofiev et al.⁸⁷, footprints of monitoring stations provide valuable information on the network fidelity. Illustrating the phenomenon, the maps of Fig. 13 show the areas constrained by the data assimilation in 1980 and 2022. A full set of years is provided as a part of the re-analysis dataset. Areas not covered by the footprints cannot be constrained by the current DA technology with the available monitoring network during the particular year. Consequently, the pollen
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Bias (left column, panels a,c) and correlation coefficient (right column, panels b,d) of the SILAM birch (upper row, panels a,b) and olive (lower row, panels c, d) forecasts in 2022, as evaluated by Copernicus Atmospheric Monitoring Service⁸⁶. Adopted from CAMS283-2023 report, available for unrestricted use from the CAMS Web portal⁸⁶.</figcaption></figure>

release over these areas is equal to its regional long-term mean. Unconstrained areas should be considered with care and excluded, for instance, from a season severity trend analysis.

The reanalysis followed the protocol suggested by the VS12 and S19 studies, which also pointed out the possibility of operational forecasting with assimilation of real-time pollen measurements. To date, the scarcity of such observations in Europe is limiting its implementation at the continental scale, but regional applications can be already feasible.

While making use of the reanalysis, one should keep in mind that the main attention was paid to the overall severity of the pollen season, i.e., the SPIn. The reanalysis did not aim to resolve day-to-day concentrations, which remained the essentially unconstrained SILAM model predictions, with their strengths and weaknesses (Fig. 12). Therefore, usage of the reanalysis for tasks sensitive to daily or hourly values (e.g., clinical trials) can be recommended only in combination with the actual pollen monitoring data. Details of the monitoring data requesting procedure can be obtained from https://www.ean-net.org/, (visited 7 Dec 2023).

## Code availability

The reanalysis has been performed with the latest operational release of the SILAM model v.5.9.1, which is freely available from https://github.com/fmidev/silam-model. The version actually used for the reanalysis is released by Kouznetsov & Tyuryakov⁸⁸.

Received: 8 January 2024; Accepted: 26 July 2024;

Published online: 03 October 2024

## References

1. D'Amato, G. et al. Thunderstorm allergy and asthma: state of the art. *Multidis Res Med* **16** (2021).

2. D'Amato, G. et al. Allergenic pollen and pollen allergy in Europe. *Allergy* **976–990**, https://doi.org/10.1111/j.1398-9995.2007.01393.x (2007).

3. de Weger, L. et al. Impact of pollen. in *Allergenic pollen. A review of the production, release, distribution and health impacts* (eds. Sofiev, M. & Bergmann, K.-C. x + 247, https://doi.org/10.1007/978-94007-4881-1 (Springer Netherlands, Dordrecht, 2013).

4. *Allergenic Pollen. A Review of Production, Release, Distribution and Health Impact*. (Springer-Verlag Berlin, Heidelberg, 2013).

5. Beggs, P. J. Thunderstorm Asthma and Climate Change. *JAMA* **331**, 878 (2024).

6. WHO. *Phenology and Human Health: Allergic Disorders*. **55** (2003).

7. Dahl, R., Andersen, P. S., Chivato, T., Valovirta, E. & De Monchy, J. National prevalence of respiratory allergic disorders. *Respiratory Medicine* **98**, 398–403 (2004).

8. Pawankar, R., Canonica, G. W., Holgate, S. T., Lockey, R. F. & Blasis, M. S. *WAO White Book on Allergy: Update 2013*. (World Allergy Organization, Milwaukee, Wisconsin, USA, 2013).
9. Savouré, M. et al. Worldwide prevalence of rhinitis in adults: A review of definitions and temporal evolution. *Clinical & Translational All* **12** (2022).

10. Thien, F. et al. The Melbourne epidemic thunderstorm asthma event 2016: an investigation of environmental triggers, effect on health services, and patient risk factors. *The Lancet Planetary Health* **2**, e255–e263 (2018).

11. Strachan, D. P. Hay fever, hygiene, and household size. *BMJ* **299**, 1259–1260 (1989).

12. Van Tilburg Bernardes, E. & Arrieta, M.-C. Hygiene Hypothesis in Asthma Development: Is Hygiene to Blame? *Archives of Medical Research* **48**, 717–726 (2017).

13. Akdis, C. A. Does the epithelial barrier hypothesis explain the increase in allergy, autoimmunity and other chronic conditions? *Nat Rev Immunol* **21**, 739–751 (2021).

14. Idrose, N. S. et al. A Review of the Respiratory Health Burden Attributable to Short-Term Exposure to Pollen. *IJERPH* **19**, 7541 (2022).

15. Jaakkola, J. J. K. et al. Airborne pollen concentrations and daily mortality from respiratory and cardiovascular causes. *European Journal of Public Health* **31**, 722–724 (2021).

16. Littlefair, J. E. et al. Air-quality networks collect environmental DNA with the potential to measure biodiversity at continental scales. *Current Biology* **33**, R426–R428 (2023).

17. Smith, M. et al. Geographic and temporal variations in pollen exposure across Europe. *Allergy* **69**, 913–23 (2014).

18. Ziska, L. H. et al. Temperature-related changes in airborne allergenic pollen abundance and seasonality across the northern hemisphere: a retrospective data analysis. *The Lancet Planetary Health* **3**, e124–e131 (2019).

19. Valipour Shokouhi, B., De Hoogh, K., Gehrig, R. & Eeftens, M. Estimation of historical daily airborne pollen concentrations across Switzerland using a spatio temporal random forest model. *Science of The Total Environment* **906**, 167286 (2024).

20. Sofiev, M. On impact of transport conditions on variability of the seasonal pollen index. *Aerobiologia* **33**, 167–179 (2016).

21. Bocquet, M. et al. Data assimilation in atmospheric chemistry models: Current status and future prospects for coupled chemistry meteorology models. *Atmospheric Chemistry and Physics* **15** (2015).

22. Elbern, H., Strunk, A., Schmidt, H. & Talagrand, O. Emission rate and chemical state estimation by 4-dimensional variational inversion. *Atmospheric Chemistry and Physics* **7**, 3749–3769 (2007).

23. Elbern, H., Schmidt, H., Talagrand, O. & Ebel, A. 4D-variational data assimilation with an adjoint air quality model for emission analysis. *Environmental Modelling & Software* **15**, 539–548 (2000).

24. Gaubert, B. et al. Regional scale ozone data assimilation using an ensemble Kalman filter and the CHIMERE chemical transport model. *Geoscience Model Development* **7**, 283–302 (2014).

25. Schwinger, J. & Elbern, H. Chemical state estimation for the middle atmosphere by four-dimensional variational data assimilation: A posteriori validation of error statistics in observation space. *Journal of Geophysical Research* **115** (2010).

26. Vira, J., Carboni, E., Grainger, R. G. & Sofiev, M. Variational assimilation of IASI SO2 plume height and total column retrievals in the 2010 eruption of Eyjafjallajökull using the SILAM v5.3 chemistry transport model. *Geoscientific Model Development* **10** (2017).

27. Vira, J. & Sofiev, M. On variational data assimilation for estimating the model initial conditions and emission fluxes for short-term forecasting of SOx concentrations. *Atmospheric Environment* **46**, 318–328 (2012).

28. Sofiev, M. On possibilities of assimilation of near-real-time pollen data by atmospheric composition models. *Aerobiologia* **1** (2019).

29. Gurney, K. R. et al. Towards robust regional estimates of CO2 sources and sinks using atmospheric transport models. *Nature* **415**, 626–630 (2002).

30. Nathan, B. J. et al. Source Sector Attribution of CO₂ Emissions Using an Urban CO/CO₂ Bayesian Inversion System. *JGR Atmospheres* **123** (2018).

31. Super, I., Dellaert, S. N. C., Visschedijk, A. J. H. & van der Gon, H. A. C. D. Uncertainty analysis of a European high-resolution emission inventory of CO2 and CO to support inverse modelling and network design. *Atmos. Chem. Phys.* **20**, 1795–1816 (2020).

32. Tans, P. P., Fung, I. Y. & Takahashi, T. Observational Contrains on the Global Atmospheric CO₂ Budget. *Science* **247**, 1431–1438 (1990).

33. Chandra, N. et al. Estimated regional CO₂ flux and uncertainty based on an ensemble of atmospheric CO₂ inversions. *Atmos. Chem. Phys.* **22**, 9215–9243 (2022).

34. Lian, J. et al. Can we use atmospheric CO₂ measurements to verify emission trends reported by cities? Lessons from a 6-year atmospheric inversion over Paris. *Atmos. Chem. Phys.* **23**, 8823–8835 (2023).

35. Nalini, K. et al. High-Resolution Lagrangian Inverse Modeling of CO₂ Emissions Over the Paris Region During the First 2020 Lockdown Period. *JGR Atmospheres* **127**, e2021JD036032 (2022).

36. Pauling, A., Clot, B., Menzel, A. & Jung, S. Pollen forecasts in complex topography: two case studies from the Alps using the numerical pollen forecast model COSMO-ART. *Aerobiologia* **36**, 25–30 (2020).

37. Siljamo, P. et al. A numerical model of birch pollen emission and dispersion in the atmosphere. Model evaluation and sensitivity analysis. *International journal of biometeorology* e-pub (2012).

38. Sofiev, M. et al. Multi-model ensemble simulations of olive pollen distribution in Europe in 2014. *Atmospheric Chemistry and Physics Discussions* 1–32, https://doi.org/10.5194/acp-2016-1189 (2017).

39. Sofiev, M. et al. A numerical model of birch pollen emission and dispersion in the atmosphere. Description of the emission module. *International journal of biometeorology* **57**, 54–58 (2012).

40. Zink, K. et al. EMPOL 1.0: a new parameterization of pollen emission in numerical weather prediction models. *Geoscience Model Development* **6**, 1961–1975 (2013).

41. Linkosalo, T. et al. A double-threshold temperature sum model for predicting the flowering duration and relative intensity of Betula pendula and B. pubescens. *Agricultural and Forest Meteorology* **150** (2010).

42. Sofiev, M., Siljamo, P., Ranta, H. & Rantio-Lehtimäki, A. Towards numerical forecasting of long-range air transport of birch pollen: Theoretical considerations and a feasibility study. *International Journal of Biometeorology* **50** (2006).

43. Galán, C. et al. Recommended terminology for aerobiological studies. *Aerobiologia* **33**, 293–295 (2017).

44. Ritenberga, O. et al. A statistical model for predicting the inter-annual variability of birch pollen abundance in Northern and Northern Europe. *Science of total environment* **615**, in press (2017).

45. Rojo, J. et al. Effects of future climate change on birch abundance and their pollen load. *Global Change Biology* **27**, 5934–5949 (2021).

46. Adamov, S. & Pauling, A. A real-time calibration method for the numerical pollen forecast model COSMO-ART. *Aerobiologia* **39**, 327–344 (2023).

47. Sofiev, M. et al. Construction of an Eulerian atmospheric dispersion model based on the advection algorithm of M. Galperin: dynamic cores v.4 and 5 of SILAM v.5.5. *Geoscientific Model Development* **8**, 3497–3522 (2015).

48. Hersbach, H. et al. The ERA5 global reanalysis. *Q.J.R. Meteorol. Soc.* **146**, 1999–2049 (2020).

49. Kouznetsov, R., Sofiev, M., Vira, J. & Stiller, G. Simulating age of air and the distribution of SF6 in the stratosphere with the SILAM model. *Atmos. Chem. Phys.* **20**, 5837–5859 (2020).

50. Sofiev, M., Genikhovich, E., Keronen, P. & Vesala, T. Diagnosing the Surface Layer Parameters for Dispersion Models within the Meteorological-to-Dispersion Modeling Interface. *Journal of Applied Meteorology and Climatology* **49**, 221–233 (2010).

51. de Rigo, D., Caudullo, G., Houston Durrant, T. & San-Miguel-Ayanz, J. *The European Atlas of Forest Tree Species: Modelling, Data and Information on Forest Tree Species*. e01aa69+, https://w3id.org/mtv/FISE-Comm/v01/e01aa69 (2016).

52. Champapeaux, J. L., Masson, V. & Chauvin, F. ECOCLIMAP: a global database of land surface parameters at 1 km resolution. *Meteorological Applications* **29**–32 (2005).
53. Prank, M. et al. An operational model for forecasting ragweed pollen release and dispersion in Europe. *Agricultural and Forest Meteorology* **182–183**, 43–53 (2013).

54. Hirst, J. M. An automatic volumetric spore trap. *Annals of Applied Biology* **39**, 257–265 (1952).

55. Jäger, S. et al. News. *Aerobiologia* **11**, 69–70 (1995).

56. Galán, C. et al. Pollen monitoring: minimum requirements and reproducibility of analysis. *Aerobiologia* **30**, 385–395 (2014).

57. CEN. *Ambient Air - Sampling and Analysis of Airborne Pollen Grains and Fungal Spores for Networks Related to Allergy—Volumetric Hirst Method*, 2019. (2019).

58. Adamov, S. et al. On the measurement uncertainty of Hirst-type volumetric pollen and spore samplers. *Aerobiologia* https://doi.org/10.1007/s10453-021-09724-5 (2021).

59. Oteros, J. et al. Errors in determining the flow rate of Hirst-type pollen traps. *Aerobiologia* **33**, 201–210 (2017).

60. Oteros, J. et al. Building an automatic pollen monitoring network (ePIN): Selection of optimal sites by clustering pollen stations. *Science of The Total Environment* **688**, 1263–1274 (2019).

61. Oteros, J. et al. Automatic and online pollen monitoring. *International Archives of Allergy and Immunology* **167**, 158–166 (2015).

62. Oteros, J. et al. An operational robotic pollen monitoring network based on automatic image recognition. *Environmental Research* **191**, 110031 (2020).

63. Maya-Manzano, J. M. et al. Towards European automatic bioaerosol monitoring: Comparison of 9 automatic pollen observational instruments with classic Hirst-type traps. *Science of The Total Environment* **866**, 161220 (2023).

64. Sofiev, M., Siljamo, P., Valkama, I., Ilvonen, M. & Kukkonen, J. A dispersion modelling system SILAM and its evaluation against ETEX data. *Atmospheric Environment* **40**, 674–685 (2006).

65. Meinander, O., Kontu, A., Kouznetsov, R. & Sofiev, M. Snow Samples Combined With Long-Range Transport Modeling to Reveal the Origin and Temporal Variability of Black Carbon in Seasonal Snow in Sodankylä (67°N). *Frontiers in Earth Science* **8**, 1–11 (2020).

66. Sofiev, M. et al. Bioaerosols in the atmosphere at two sites in Northern Europe in spring 2021: Outline of an experimental campaign. *Environmental Research* **214**, 113798 (2022).

67. Sofiev, M. et al. MACC regional multi-model ensemble simulations of birch pollen dispersion in Europe. *Atmospheric Chemistry and Physics* **15** (2015).

68. Kouznetsov, R. & Sofiev, M. A methodology for evaluation of vertical dispersion and dry deposition of atmospheric aerosols. *Journal of Geophysical Research* **117** (2012).

69. Sofiev, M. Extended resistance analogy for construction of the vertical diffusion scheme for dispersion models. *Journal of Geophysical Research-Atmospheres* **107**, ACH 10-1-ACH 10-8 (2002).

70. Brasseur, G. P. et al. Ensemble forecasts of air quality in eastern China – Part 1: Model description and implementation of the MarcoPolo-Panda prediction system, version 1. *Geosci. Model Dev.* **12**, 33–67 (2019).

71. Huijnen, V. et al. Comparison of OMI NO2 tropospheric columns with an ensemble of global and European regional air quality models. *Atmospheric Chemistry and Physics* **10**, 3273–3296 (2010).

72. Petersen, A. K. et al. Ensemble forecasts of air quality in eastern China – Part 2: Evaluation of the MarcoPolo-Panda prediction system, version 1. *Geosci. Model Dev.* **12**, 1241–1266 (2019).

73. Sofiev, M., Kouznetsov, R., Hänninen, R. & Sofieva, V. F. Technical note: Intermittent reduction of the stratospheric ozone over northern Europe caused by a storm in the Atlantic Ocean. *Atmos. Chem. Phys.* **20**, 1839–1847 (2020).

74. Xian, P. et al. Current state of the global operational aerosol multi-model ensemble: An update from the International Cooperative for Aerosol Prediction (ICAP). *Q.J.R. Meteorol. Soc.* **145**, 176–209 (2019).

75. Siljamo, P., Ashbrook, K., Comont, R. F. & Skjoth, C. A. Do atmospheric events explain the arrival of an invasive ladybird (Harmonia axyridis) in the UK? *PLoS ONE* **15**, e0219335 (2020).

76. Prank, M., Sofiev, M., Siljamo, P., Kauhaniemi, M. & European Aeroallergen Network. Increasing the Number of Allergenic Pollen Species in SILAM Forecasts. in *Air Pollution Modeling and its Application XXIV*, edited by D. G. Steyn and N. Chaumerliac 313–317 (Springer, 2016).

77. Tummon, F. et al. Towards standardisation of automatic pollen and fungal spore monitoring: best practises and guidelines. *Aerobiologia* https://doi.org/10.1007/s10453-022-09755-6 (2022).

78. Hansen, P. C. *Discrete Inverse Problems: Insight and Algorithms*. https://doi.org/10.1137/1.9780898718836 (Society for Industrial and Applied Mathematics, 2010).

79. Hansen, P. C. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. *SIAM Rev.* **34**, 561–580 (1992).

80. Sofiev, M. et al. European pollen reanalysis, 1980-2022, for alder, birch, and olive, v.1.1. Finnish Meteorological Institute https://doi.org/10.57707/FMI-B2SHARE.85841086F9DB46B882D750EAA9E42515 (2024).

81. Morgado, R. et al. Drivers of irrigated olive grove expansion in Mediterranean landscapes and associated biodiversity impacts. *Landscape and Urban Planning* **225**, 104429 (2022).

82. Veriankaitė, L., Siljamo, P., Sofiev, M., Sauliene, I. & Kukkonen, J. Modelling analysis of source regions of long-range transported birch pollen that influences allergenic seasons in Lithuania. *Aerobiologia* **26**, 47–62 (2010).

83. Verstraeten, W. W. et al. Reconstructing multi-decadal airborne birch pollen levels based on NDVI data and a pollen transport model. *Agricultural and Forest Meteorology* **320**, 108942 (2022).

84. Verstraeten, W. W., Kouznetsov, R., Bruffaerts, N., Sofiev, M. & Delcloo, A. W. Assessing uncertainty in airborne birch pollen modelling. https://doi.org/10.1007/s10453-024-09818-w (2024).

85. Verstraeten, W. W. et al. Attributing long-term changes in airborne birch and grass pollen concentrations to climate change and vegetation dynamics. *Atmospheric Environment* **298**, 119643 (2023).

86. Sofiev, M., Palamarchuk, J., Kouznetsov, R., Gauss, M. & CAMS modelling teams. *Annual Report on the Evaluation of the CAMS Regional Pollen Production (Daily Forecasts)*. January 2022–October 2022. 67. https://atmosphere.copernicus.eu/sites/default/files/custom-uploads/EQC-regional/Pollen/CAMS_pollen_eval_2022_v3.pdf (2024).

87. Sofiev, M. et al. Designing an automatic pollen monitoring network for direct usage of observations to reconstruct the concentration fields. *Science of The Total Environment* **16580**, https://doi.org/10.1016/j.scitotenv.2023.165800 (2023).

88. Kouznetsov, R. Silam_v5_9 used for pollen reanalysis 2023. Zenodo https://doi.org/10.5281/ZENODO.10351493 (2023).

## Acknowledgements

The authors acknowledge support of the following research funding: Horizon projects SYLVA (grant nr 101086109), CATALYSE (nbr 101057131), EO4EU (nbr 101060784), Copernicus CAMS2_40, Academy of Finland projects PS4A (nbr 318194), ALL-Impress (nbr 329215), SPORELIFE (nbr 355851) and Latvian science base funding. The following European aerobiological networks provided their data for the reanalysis: Europe: EAN, European Aeroallergen Network. Austria: Austrian aerobiological network. Belarus: Minsk aerobiological stations. Belgium: AirAllergy, Belgian aerobiological network. Croatia: Croatian aerobiological network. Czech Republic: Czech aerobiological network. Denmark: Danish Aerobiological network. Estonia: Estonian aerobiological network. Finland: Finnish Aerobiological network. France: RNA, French aerobiological
monitoring network. Georgia: Georgian aerobiological network. Germany: ePIN, Electronic Pollen Information Network, Bayerische Landesamt für Gesundheit und Lebensmittelsicherheit (LGL). Greece: Hellenic aerobiological network. Hungary: Hungarian Aerobiological Network. Iceland: Icelandic aerobiological network. Israel: Israeli aerobiological network. Italy: Italian aerobiological network. Latvia: Riga aerobiological station. Lithuania: Lithuanian aerobiological network. Luxembourg: Luxemburg aerobiological station. Netherlands: Dutch aerobiological network. Norway: Norwegian aerobiological network. Poland: Polish aerobiological network. Portugal: Portuguese aerobiological network. Romania: Romanian aerobiological network. Russia: Moscow aerobiological station. Serbia: Serbian state network for allergenic pollen measurement, City of Novi Sad; AP Vojvodina. Slovakia: Slovakian aerobiological network. Slovenia: Slovenian aerobiological network. Spain: REA Spanish Aerobiology Network, RAA Andalusian Aerobiology Network, Madrid Region Palynological Network, XAC Catalan Aerobiological Network, Extremadura network. Sweden: Swedish aerobiological network. Switzerland: Swiss aerobiological network, R. Andenmatten, R.M. Leuschner, A.G.Peeters, and L. Wick are acknowledged for pollen data prior to 1993. Türkiye: Turkish aerobiological stations. UK: MetOffice, UK. Ukraine: Ukrainian aerobiological network.

Author contributions

All authors have participated in editing and discussing the paper, have read and approved its final version. Specific contributions: Sofiev, Mikhail: study design and implementation, drafting the paper, DA module development, data processing and analysis, SILAM pollen module development. Palamarchuk, Julia: study design and implementation, drafting the paper, data processing and analysis, SILAM pollen module development. Kouznetsov, Rostislav: study implementation, DA module development, data processing and analysis. Abramidze, Tamuna: observational data preparation, processing and evaluation. Adams-Groom, Beverley: observational data preparation, processing and evaluation. Antunes, Célia M.: observational data preparation, processing and evaluation. Ariño, Arturo H.: observational data preparation, processing and evaluation. Bastl, Maximilian: observational data preparation, processing and evaluation. Belmonte, Jordina: observational data preparation, processing and evaluation. Berger, Uwe: EAN system maintenance and data analysis and quality assurance, observational data preparation, processing and evaluation. Bonini, Maira: observational data preparation, processing and evaluation. Bruffaerts, Nicolas: observational data preparation, processing and evaluation. Buters, Jeroen: observational data preparation, processing and evaluation. Carinanos, Paloma: observational data preparation, processing and evaluation. Celenk, Sevcan: observational data preparation, processing and evaluation. Ceriotti, Valentina: observational data preparation, processing and evaluation. Charalampopoulos, Athanasios: observational data preparation, processing and evaluation. Clewlow, Yolanda: observational data preparation, processing and evaluation. Clot, Bernard: observational data preparation, processing and evaluation. Dahl, Aslog: observational data preparation, processing and evaluation. Damialis, Athanasios: observational data preparation, processing and evaluation. De Linares, Concepción: observational data preparation, processing and evaluation. De Weger, Letty A.: observational data preparation, processing and evaluation. Dirr, Lukas: EAN system maintenance, data quality assurance, observational data preparation. Ekebom, Agneta: observational data preparation, processing and evaluation. Fatahi, Yalda: SILAM model evaluation, data processing and analysis. Fernández González, María: observational data preparation, processing and evaluation. Fernández González, María Delia: observational data preparation, processing and evaluation. Fernández-Rodríguez, Santiago: observational data preparation, processing and evaluation. Galán, Carmen: observational data preparation, processing and evaluation. Gedda, Björn: observational data preparation, processing and evaluation. Gehrig, Regula: observational data preparation, processing and evaluation. Geller Bernstein, Carmi: observational data preparation, processing and evaluation. Gonzalez, Roldan Nestor: observational data preparation, processing and evaluation. Grewling, Lukasz: observational data preparation, processing and evaluation. Hajkova, Lenka: observational data preparation, processing and evaluation. Hänninen, Risto: SILAM model development and evaluation. Hentges, François: observational data preparation, processing and evaluation. Jantunen, Juha: observational data preparation, processing and evaluation. Kadantsev, Evgeny: SILAM model development and evaluation. Kaspryzk, Idalia: observational data preparation, processing and evaluation. Kloster, Mathilde: observational data preparation, processing and evaluation. Kluska, Katarzyna: observational data preparation, processing and evaluation. Koenders, Mieke: observational data preparation, processing and evaluation. Lafférvová, Janka: observational data preparation, processing and evaluation. Leru, Poliana: observational data preparation, processing and evaluation. Lipiec, Agnieszka: observational data preparation, processing and evaluation. Louna-Korteniemi, Maria: observational data preparation, processing and evaluation. Magyar, Donát: observational data preparation, processing and evaluation. Majkowska-Wojciechowska, Barbara: observational data preparation, processing and evaluation. Mitrovic, Mirjana: observational data preparation, processing and evaluation. Myszkowska, Dorota: observational data preparation, processing and evaluation. Mäkelä, Mika: observational data preparation, processing and evaluation. Oliver, Gilles: observational data preparation, processing and evaluation. Östensson, Pia: observational data preparation, processing and evaluation. Pätsi, Sanna: observational data preparation, processing and evaluation. Pérez-Badia, Rosa: observational data preparation, processing and evaluation. Piotrowska-Weryszko, Krystyna: observational data preparation, processing and evaluation. Prank, Marje: SILAM pollen modules development, SILAM model development and evaluation. Przedpelska-Wasowicz, Ewa Maria: observational data preparation, processing and evaluation. Rajo, F. Javier Rodríguez: observational data preparation, processing and evaluation. Ramfjord, Hallvard: observational data preparation, processing and evaluation. Rapiejko, Joanna: observational data preparation, processing and evaluation. Rodinkova, Victoria: observational data preparation, processing and evaluation. Rojo, Jesús: observational data preparation, processing and evaluation. Ruiz-Valenzuela, Luis: observational data preparation, processing and evaluation. Rybníček, Ondrej: observational data preparation,
processing and evaluation. Saarto, Annika: observational data preparation, processing and evaluation. Sauliene, Ingrida: observational data preparation, processing and evaluation. Seliger, Andreja Kofol: observational data preparation, processing and evaluation. Severova, Elena: observational data preparation, processing and evaluation. Shalaboda, Valentina: observational data preparation, processing and evaluation. Sikoparija, Branko: observational data preparation, processing and evaluation. Siljamo, Pilvi: SILAM pollen modules development and evaluation. Soares, Joana: SILAM model application and evaluation. Sozinova, Olga: observational data preparation, processing and evaluation. Stangel, Anders: SILAM model development and operationalization. Stjepanović, Barbara: observational data preparation, processing and evaluation. Teinemaa, Erik: observational data preparation, processing and evaluation. Tyuryakov, Svyatoslav: SILAM model distribution, project management. Uppstu, Andreas: SILAM model development and evaluation. SILAM DA modules development. Vill, Mart: observational data preparation, processing and evaluation. Vira, Julius: SILAM model development and evaluation. SILAM DA modules development. Visez, Nicolas: observational data preparation, processing and evaluation. Vitikainen, Tiina: observational data preparation, processing and evaluation. Vokou, Despoina: observational data preparation, processing and evaluation. Weryszko-Chmielewska, Elżbieta: observational data preparation, processing and evaluation. Karppinen, Ari: financial, managerial, and administrative support and coordination of the SILAM team.

**Competing interests**

The authors declare no competing interests.

**Additional information**

Correspondence and requests for materials should be addressed to M.S.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_17.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre-
ative Commons licence, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons licence and your intended use is not per-
mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the
copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024

Mikhail Sofiev¹, Julia Palamarchuk¹, Rostislav Kouznetsov¹, Tamuna Abramidze²,
Beverley Adams-Groom³, Célia M. Antunes⁴, Arturo H. Ariño⁵, Maximilian Bastl⁶,
Jordina Belmonte⁷,⁸, Uwe E. Berger⁹, Maira Bonini¹⁰, Nicolas Bruffaerts¹¹,
Jeroen Buters¹², Paloma Cariñanos¹³,¹⁴, Sevcan Celenk¹⁵, Valentina Ceriotti¹⁰,
Athanasios Charalampopoulos¹⁶, Yolanda Clewlow¹⁷, Bernard Clot¹⁸, Aslog Dahl¹⁹,
Athanasios Damialis¹⁶, Concepción De Linares¹³, Letty A. De Weger²⁰, Lukas Dirr⁶,
Agneta Ekebom²¹, Yalda Fatahi¹, María Fernández González²², Delia Fernández
González²³,²⁴, Santiago Fernández-Rodríguez²⁵, Carmen Galán²⁶, Björn Gedda²¹,
Regula Gehrig²⁸, Carmi Geller Bernstein²⁷, Nestor Gonzalez Roldan²⁸, Lukasz Grewling²⁹,
Lenka Hajkova³⁰, Risto Hänninen¹, François Hentges³¹, Juha Jantunen³²,
Evgeny Kadantsev¹, Idalia Kasprzyk³³, Mathilde Kloster³⁴, Katarzyna Kluska³⁵,
Mieke Koenders³⁶, Janka Lafférsová³⁷, Poliana Mihaela Leru³⁸,³⁹, Agnieszka Lipiec⁴⁰,
Maria Louna-Korteniemi⁴¹, Donát Magyar⁴², Barbara Majkowska-Wojciechowska⁴³,⁴⁴,
Mika Mäkelä⁴⁵, Mirjana Mitrovic⁴⁶, Dorota Myszkowska⁴⁷, Gilles Oliver⁴⁸,
Pia Östensson²¹, Rosa Pérez-Badia⁴⁹, Krystyna Piotrowska-Weryszko⁵⁰, Marje Prank¹,
Ewa Maria Przedpelska-Wasowicz⁵¹, Sanna Pätsi⁴¹, F. Javier Rodríguez Rajo²²,
Hallvard Ramfjord⁵², Joanna Rapiejko⁵³, Victoria Rodinkova⁵⁴, Jesús Rojo⁵⁵,
Luis Ruiz-Valenzuela⁵⁶,⁵⁷, Ondrej Rybnicek⁵⁸,⁵⁹, Annika Saarto⁴¹, Ingrida Sauliene⁶⁰,
Andreja Kofol Seliger⁶¹, Elena Severova⁶²,⁶³, Valentina Shalaboda⁶⁴, Branko Sikoparija⁶⁵,
Pilvi Siljamo¹, Joana Soares⁶⁶, Olga Sozinova⁶⁷, Anders Stangel¹, Barbara Stjepanović⁶⁸,
Erik Teinemaa⁶⁹, Svyatoslav Tyuryakov¹, M. Mar Trigo⁷⁰, Andreas Uppstu¹,
Mart Vill⁶⁹, Julius Vira¹, Nicolas Visez⁴⁸,⁷¹, Tiina Vitikainen³², Despoina Vokou¹⁶,
Elżbieta Weryszko-Chmielewska⁵⁰ & Ari Karppinen¹
¹Finnish Meteorological Institute, Helsinki, Finland. ²Center of Allergy & Immunology, Tbilisi, Georgia. ³University of Worcester, School of Science and Environment, Worcester, UK. ⁴University of Évora, School of Health and Human Development, Department of Medical and Health Sciences & Institute of Earth Sciences - ICT, Évora, Portugal. ⁵University of Navarra, Biodiversity and Environment Institute, Pamplona, Spain. ⁶Department of Otorhinolaryngology, Medical University of Vienna, Vienna, Austria. ⁷Departament de Biologia Animal, Biologia Vegetal i Ecologia, Universitat Autònoma de Barcelona (UAB), Bellaterra, Spain. ⁸Institut de Ciència i Tecnologia Ambientals (ICTA-UAB), Universitat Autònoma de Barcelona, Bellaterra, Spain. ⁹University of Innsbruck, Department of Botany, Innsbruck, Austria. ¹⁰Department of Hygiene and Health Prevention, Agency for Health Protection of Metropolitan Area of Milan (ATS), Milan, Italy. ¹¹Mycology and Aerobiology, Sciensano, Brussels, Belgium. ¹²Center of Allergy & Environment (ZAUM), Member of the German Center for Lung Research (DZL), Technical University and Helmholtz Center Munich, Munich, Germany. ¹³Department of Botany, University of Granada, Granada, Spain. ¹⁴Andalusian Institute for Earth System Research (IISTA-CEAMA), University of Granada, Granada, Spain. ¹⁵Bursa Uludag University, Faculty of Arts and Science, Department of Biology, Aerobiology Laboratory, 16059, Görükle-Bursa, Türkiye. ¹⁶Department of Ecology, School of Biology, Faculty of Sciences, Aristotle University of Thessaloniki, Thessaloniki, Greece. ¹⁷Health, air quality, & UK pollen forecasting, UK Met Office, Exeter, UK. ¹⁸Federal Office of Meteorology and Climatology MeteoSwiss, Zurich, Switzerland. ¹⁹Department of Biological and Environmental Sciences, University of Gothenburg, Gothenburg, Sweden. ²⁰Department of Pulmonology, Leiden University Medical Center, Leiden, the Netherlands. ²¹Palynological Laboratory, Swedish Museum of Natural History, Stockholm, Sweden. ²²Sciences Faculty, University of Vigo, Ourense, 32002, Spain. ²³Biodiversity and Environmental Management, University of León, León, Spain. ²⁴Institute of Atmospheric Sciences and Climate-CNR, Bologna, Italy. ²⁵Department of Construction, School of Technology, University of Extremadura, Avda. de la Universidad s/n, Cáceres, Spain. ²⁶Inter-University Institute for Earth System Research (IISTA), International Campus of Excellence on Agri-food (ceiA3), University of Cordoba, Cordoba, Spain. ²⁷Zabludovicz Center for Autoimmune Diseases, Sheba Medical Center, Ramat Gan, Israel. ²⁸Pollen Laboratory, Department of Biological and Environmental Sciences, University of Gothenburg, Gothenburg, Sweden. ²⁹Laboratory of Aerobiology, Department of Systematic and Environmental Botany, Faculty of Biology, Adam Mickiewicz University, Poznan, Poland. ³⁰Czech Hydrometeorological Institute, Prague, Czech Republic. ³¹Unit of Immunology-Allergology, Centre Hospitalier de, Luxembourg, Luxembourg. ³²South Karelia Allergy and Environment Institute, Imatra, Finland. ³³College of Natural Sciences University of Rzeszow, Rzeszow, Poland. ³⁴The Asthma and Allergy Association, Roskilde, Denmark. ³⁵Institute of Biology, College of Natural Sciences University of Rzeszow, Rzeszow, Poland. ³⁶Elkerliek Helmond, Helmond, Netherlands. ³⁷Regional Public Health Office department of medical microbiology, bratislava, Slovakia. ³⁸Clinical Department 5, Carol Davila University of Medicine, Bucharest, Romania. ³⁹Allergology Research Laboratory, Colentina Clinical Hospital, Bucureşti, Romania. ⁴⁰Department of the Prevention of Environmental Hazard, Allergology and Immunology, Medical University of Warsaw, Warsaw, Poland. ⁴¹Biodiversity Unit, University of Turku, Turku, Finland. ⁴²National Center for Public Health and Pharmacy, Budapest, Hungary. ⁴³Aeroallergen Monitoring Centre "AMoC", Department of Immunology and Allergy, Allergy, Poland. ⁴⁴Medical University of Lodz, Lodz, Poland. ⁴⁵HUS Helsingin yliopistollinen sairaala, Jyväskylä, Finland. ⁴⁶Serbian Environmental Protection Agency, Belgrade, Serbia. ⁴⁷Jagiellonian University Medical College, Department of Clinical and Environmental Allergology, Kraków, Poland. ⁴⁸French Aerobiological Monitoring Network (RNSA), Brussieu, France. ⁴⁹University of Castilla-La Mancha, Institute of Environmental Sciences, Toledo, Spain. ⁵⁰Department of Botany and Plant Physiology, Subdepartment of Aerobiology, University of Life Sciences in Lublin, Lublin, Poland. ⁵¹Icelandic Institute of Natural History, Akureyri, Iceland. ⁵²Department of Biology, NTNU, Trondheim, Norway. ⁵³Allergen Research Center, Warsaw, Poland. ⁵⁴Department of Pharmacy, National Pirogov Memorial Medical University, Vinnytsia, Ukraine. ⁵⁵Department of Pharmacology, Pharmacognosy and Botany, Faculty of Pharmacy, Complutense University of Madrid, Madrid, Spain. ⁵⁶Department of Biology Animal, Plant Biology and Ecology, University of Jaén, Jaén, Spain. ⁵⁷University Institute of research in Olive Groves and Olive Oils, University of Jaén, Jaén, Spain. ⁵⁸University Hospital Brno, Brno, Czech Republic. ⁵⁹Masaryk University, Brno, Czech Republic. ⁶⁰Vilnius University Siauliai Academy, Siauliai, Lithuania. ⁶¹National Laboratory of Health, Environment and Food, Maribor, Slovenia. ⁶²Faculty of Biology, Moscow State University, Moscow, Russia. ⁶³Faculty of Biology, Shenzhen MSU-BIT University, Shenzhen, China. ⁶⁴Retired from Faculty of Pharmacy of the Belarusian State Medical University, Minsk, Belarus. ⁶⁵BioSense Institute Research Institute for Information Technologies in Biosystems, University of Novi Sad, Novi Sad, Serbia. ⁶⁶Stiftelsen NILU - Stiftelsen Norwegian Institute for Air Research, Kjeller, Norway. ⁶⁷University of Latvia, Riga, Latvia. ⁶⁸Laboratory of Aerobiology at Teaching Institute of Public Health dr. Andrija Štamper, Zagreb, Croatia. ⁶⁹Estonian Environmental research Institute (under Estonian Environmental Research Centre), Tartu, Estonia. ⁷⁰Department of Botany and Plant Physiology, University of Malaga, Malaga, Spain. ⁷¹Université de Lille, CNRS, UMR, 8516, LASIRE - Laboratoire de Spectroscopie pour les Interactions, la Réactivité et l'Environnement, F-59000, Lille, France. e-mail: mikhail.sofiev@fmi.fi