Regulatory Air Pollution Modelling in Poland: Is It Time for an Update?

Essentially, we have three options when selecting meteorological data necessary for modelling pollutant dispersion on a local scale:

a) Using long-term records of measurements and observations taken at surface meteorological stations, schematically classified into a set of scenarios (“meteorological situations”) with their corresponding frequencies of occurrence. This option reduces the computational cost by an order or two orders of magnitude compared to the direct use of hourly data.

b) Direct use of meteorological observations. Here, we deal with two datasets: hourly measurements and observations from surface synoptic stations (“SYNOP” messages) and upper-air data from aerological observatories (“TEMP”) – typically twice per day. The public datasets available from the Polish weather service include SYNOP data from roughly 60 stations and TEMP data from four observatories. Although radar measurements of wind speed are currently available, these data – to the best knowledge of the authors have not yet been used for computing pollutant dispersion in Poland.

The direct use of meteorological observations entered regulatory air pollution modelling at the turn of the 20th and 21st centuries, still in the context of “updated” Gaussian plume models (e.g. Weil and Brower [1984]), significantly improving the quality of the results. However, most of the limitations discussed in the previous section remain in effect. Moreover, the representativeness of data from distant stations is questionable, especially in the case of atmospheric soundings.

c) Using data from meteorological models. Here, climatological reanalysis products, which span extended time periods, play the most crucial role. By utilising historical data, they provide a continuous picture of the variability of meteorological conditions over time and space — unlike the point measurements discussed previously. However, the resolution of computational grids used in reanalyses is still insufficient in spatially diverse areas, such as mountains, coastlines, and cities. For this reason, nested high-resolution models are additionally employed in air quality studies, focusing on much smaller areas. This approach, however, requires significant computational resources and a high level of expertise. In Poland, it has been used in the development of air protection programs. Ideally, these calculations should also include data assimilation of surface observations to reduce divergences between the model-based climatology and the observational data. Unfortunately, in many studies conducted so far in the country, this step has been neglected, leading in many cases to substantial discrepancies.

The most accurate results should be obtained by choosing the last option. However, due to high computational and labour costs, and the vast volume of intermediate results generated, it is reasonable to limit high-resolution calculations to areas and periods where air quality standard exceedances have been recorded by monitoring stations. One may also note that with the progress of climate change, the consultation of long-term historical climatology becomes more and more questionable.

Emission data are often considered the primary source of uncertainty in air pollution modelling. In fact, most of this information is based on estimates rather than measurements, and different estimation methods lead to considerable discrepancies in the results. Moreover, these data are not available in real-time, and therefore calculations for current periods rely on data from previous years. In regulatory applications, declared values are used.

The lack of sufficiently accurate emission data is often cited as a justification for using oversimplified models or methods. On the other hand, using a reliable model can provide feedback to emission estimation when the model results are compared to measurements.

Terrain features play important role in the transport, dispersion and deposition of air pollutants. Topography modifies air flow, and contributes to the occurrence of spatially varying thermal stratifications in the boundary layer due to non-homogeneous exposure of air to solar radiation. The type of land cover determines its roughness, which directly affects the vertical wind profile, as well as the efficiency of turbulent mixing processes in the lowest layer of the atmosphere. Thus, the efficiency of removing air pollutants by dry subsidence depends, among other things, on the land-use type. The extent and the way in which the above-mentioned processes are taken into account in air pollution modelling depends on the class of the model.

The Polish regulatory model — i.e., the Gaussian plume model — makes simplifying assumptions (flat, homogeneous terrain, steady-state conditions), and so information about the terrain type is taken into account only during the calculation of wind profile and atmospheric diffusion (dispersion) coefficients σ_y and σ_z. The constants A and B appearing in (2) are calculated from the following formulas, developed by Nowicki [1976]: (3) A=0,088*6*m−0,3+1−lnHz0 A = 0,088*\left( {6*{m^{ - 0,3}} + 1 - ln{H \over {{z_0}}}} \right) (4) B=0,38*m1,3*8,7−lnHz0 B = 0,38*{m^{1,3}}*\left( {8,7 - ln{H \over {{z_0}}}} \right) where m is a parameter describing the atmospheric equilibrium state, H is the effective height of the emitter, and z0 is the aerodynamic surface roughness length.

In advanced air quality models, meteorological and terrain data are entered into specialized modules known as preprocessors, whose function is to prepare in an appropriate format the full input information needed for the air quality model. For example, in the AERMOD system [Cimorelli et al. 2005], a separate terrain preprocessor (AERMAP) provides a physical relationship between terrain features and the behaviour of pollutant plumes. It generates location and height data for each receptor location while providing information that allows the dispersion model to account for the effects of stability-dependent flow patterns around or over hills.

In order to take into account the influence of surface conditions on boundary layer parameter estimates, the surface characteristics (surface albedo, Bowen ratio and effective surface roughness length) are entered into the meteorological pre-processor (AERMET). This calculates the atmospheric parameters needed by the dispersion model, such as atmospheric turbulence characteristics, mixing height, friction velocity, Monin-Obukov length, and surface heat flux.

Another example is the GRAL [Oettl 2015b], a Lagrangian particle dispersion model that is discussed later in this article. This model can be coupled with GRAMM [Oettl 2015a], an Eulerian prognostic non-hydrostatic mesoscale wind field model. It features a terrain-following grid, allowing the effects of terrain on dispersion (e.g., cold-air drainage flows) to be taken into account at all scales. GRAMM solves the conservation equations for mass, enthalpy, momentum, and humidity. A radiation model to take long- and short-wave radiation into account also exists. The surface-energy balance is calculated in a surface module, where several different land-use categories are used to define surface roughness; albedo; emissivity; soil moisture content; specific heat capacity of the soil; and the heat transfer coefficient.

As previously mentioned, in Lagrangian Particle Dispersion (LPD) models, the movement of pseudoparticles is a combination of average flow (wind) and turbulent motions, i.e., (5) dr→dt=v→¯r→,v→+v→′r→,v→=v→r→,v→ \frac{d\vec{r}}{dt}=\bar{\vec{v}}\left( \vec{r},\vec{v} \right)+\vec{v}'\left( \vec{r},\vec{v} \right)=\vec{v}\left( \vec{r},\vec{v} \right) where the first term on the right-hand side represents the velocity of the mean flow, and the second term represents the turbulent fluctuation, i.e., the random deviation from the mean value. Integrating this equation allows for the determination of a pseudoparticle’s trajectory. The term describing fluctuations is typically modelled as a random process through a stochastic differential equation of the form: (6) dv→=a→r→,v→dt+b→r→,v→dW d\vec{v}=\vec{a}\left( \vec{r},\vec{v} \right)dt+\vec{b}\left( \vec{r},\vec{v} \right)dW where the second term on the right-hand side describes random changes in velocity through the increments of the Wiener process dW, with a mean value of zero and variance dt. In homogeneous turbulence, this is the Langevin equation with (7) a→=−v→¯/t and b→=2v→¯2/τ \vec{a}=-\bar{\vec{v}}/t\ \ \ \text{and}\ \ \ \vec{b}=\sqrt{2{{{\bar{\vec{v}}}}^{2}}/\tau } where τ is the Lagrangian timescale of turbulence. With the random component set to zero, the first of the described equations describes the changes in the position of plume centres in the models discussed earlier.

The most important differences in the construction of various models of this class concern the method of modelling fluctuations (see, e.g., Oettl et al. [2001]), and the form of the kernel of the integral operator used to determine concentrations.

Among LPD models, some have been applied for regulatory purposes, and we will focus on those here. The strongest legal backing for this model is found in Germany, where a model of this class (AUSTAL2000) has served as the reference model designated for permitting processes and environmental impact assessments since 2002. ⁽³⁾ Procedures related to air protection are described in the TA Luft technical instructions published by the Ministry of the Environment, while the algorithm implemented in the model is detailed in the VDI 3945 standard [VDI Guideline… 2000]. The latest version of the model (AUSTAL v.3, the successor to the program AUSTAL2000) refers to the Annex 2 and 7 of the new version of TA Luft, which entered into force in December 2021 [TA Luft 2021].

Both the source code and executables for Windows and Linux are made available free of charge by the German Federal Environment Agency (Umweltbundesamt). ⁽⁴⁾ According to the documentation [Janicke 2024], the program supports both calculations using measurement series and the classification of meteorological situations. It allows for calculations under time-varying emission conditions, in complex terrain (through the use of the kinematic wind field model TALDIA), and around buildings; it also can assess exposure to odorous substances. In the latest version, the modelling of the wet deposition of pollutants is also available.

Unfortunately, the scientific documentation for the model is rather limited, and aside from the aforementioned VDI standard, it is difficult to find scientific publications on the model’s structure and its validation. Furthermore, the latter has faced sharp criticism [Schenk 2020]. Additionally, as noted by Langner and Klemm [2011], a comparison of modelling results with measurement data shows significantly better performance for the updated Gaussian model AERMOD than for AUSTAL2000. They found that AERMOD was more efficient for modelling pollutant dispersion in complex terrain, likely due to its enhanced algorithms for handling topographic influences on pollutant flow. On the other hand, AUSTAL2000 was particularly effective when comprehensive meteorological input data were unavailable, thanks to its capability to simulate dispersion under various data-limited conditions. The study highlighted the fact that AERMOD was significantly faster in terms of computation time, making it more suitable for applications requiring quick assessments. It is also worth noting that the adopted flow model (TALDIA), based on the variational method of minimising deviations at measurement points, produces results that are highly dependent on the initial approximation. It should, therefore, be used more as an assimilation method in combination with a mesoscale meteorological model or hydrodynamic model, rather than as a standalone module.

A particularly interesting example is the Austrian GRAL/GRAMM system, developed by the Graz University of Technology (Technische Universität Graz). The work on this system, motivated by the need to create modelling methods suitable for use under low wind conditions or in complex terrain, began in 1999. In 2006, a collaboration was established between Technische Universität Graz and the authorities of the state of Styria, ⁽⁵⁾ resulting in a system composed of two models: the mesoscale prognostic model, GRAMM (Graz Mesoscale Model [Oettl 2015a]) and the Lagrangian particle dispersion model, GRAL (Graz Lagrangian Model [Oettl 2015b]).

The GRAMM model is a non-hydrostatic mesoscale prognostic model based on averaged conservation equations of mass, momentum, and energy in the anelastic or Boussinesq approximation formulated in a Cartesian coordinate system. It features alternative parameterisations for turbulent transport: (i) a simple Pandolfo parameterisation based on deformation velocity, Richardson number, and mixing length; (ii) a turbulent kinetic energy (TKE) budget-based closure; and (iii) the standard k-ε model, which includes equations for the balance of turbulent kinetic energy and dissipation rate. Surface exchange processes are described using the Monin-Obukhov similarity theory and a parameterised surface heat balance. The current version of the model (23.11) does not include condensation, radiation, or atmospheric chemistry modules. The model’s validation covers 22 test scenarios described in the German standard VDI 3783-9 [VDI Guideline… 2017], which involve comparisons with wind tunnel experiment results. Although GRAMM appears fairly basic when compared to other modern mesoscale models, it provides automated integration with the GRAL model, allowing it to be used by properly-trained individuals who are not necessarily specialists in numerical atmospheric dynamics modelling.

The Lagrangian particle dispersion model GRAL possesses all the advantages of the AUSTAL2000 model and is scientifically supported by numerous publications, such as those demonstrating its effective application in complex urban environments and high-resolution simulations [Berchet et al. 2017; Petrov et al. 2024; Patiño et al. 2024]. It also features comprehensive and detailed documentation covering both the physical principles and their validation (Oettl 2024).

Validation of this model included four experiments as required by the Austrian guideline RVS 04.02.12, which establishes quality criteria for dispersion models, requiring a Normalized Mean Square Error (NMSE) of 3.0 or less and a Fraction Bias (FB) of 0.3 or less. The GRAL model was evaluated on these four datasets to assess its performance and compliance with these criteria. The first dataset, CALTRANS99, focused on the dispersion of a tracer gas (SF6) along Highway 99 in California, which is characterized by open terrain and low wind speeds. Emissions were monitored from vehicles travelling in a loop along the highway. GRAL produced an NMSE of 0.5 and an FB of 0.0, meeting the guideline criteria. Although it overestimated peak concentrations, it performed comparably to other models such as ADMS-ROADS and AERMOD.

The second dataset, A2 Biedermannsdorf, monitored NOx concentrations near the A2 highway in Austria, taking into account high traffic volumes and background concentrations. This site included obstacles, such as a six-meter noise barrier wall. GRAL showed good agreement at most observation points, though slight underestimation was noted 400 metres east of the highway. NMSE values ranged between 0.0 and 0.2, while mean deviations ranged from −0.2 to 0.2, demonstrating compliance with the guideline.

The third dataset, Ehrentalerberg Tunnel, focused on SF₆ dispersion from a tunnel portal under low wind speed conditions. It simulated complex jet stream behaviour caused by meandering wind patterns. GRAL achieved an NMSE of 0.9 and a FB of 0.1 when using the observed turbulence data, which met the guideline criteria. A second simulation using stability classes produced a higher NMSE of 2.2, but remained compliant. GRAL effectively modelled mean concentrations and concentration statistics.

The fourth dataset, Kaisermuehlen Tunnel, examined NOx dispersion from a 2,150-m tunnel in Vienna that was characterized by complex topography and traffic patterns. The evaluation covered two distinct wind direction scenarios over a 10-month period. GRAL demonstrated variable performance across sites, with NMSE values ranging from 0.4 to 2.4 and FB values between −0.3 and 0.3. While average concentrations were captured well, peak values were sometimes underestimated. Results complied with guidelines, although higher uncertainty was observed due to rough emission estimates.

Overall, GRAL satisfied the Austrian quality criteria across all datasets, confirming its suitability for pollutant dispersion modelling in Austria. Although localized deviations and uncertainties were noted, especially in complex terrain, the model demonstrated reliable performance under varied conditions.

Validation of the GRAL model also includes 15 test cases, according to the German VDI 3783-9, and 26 diffusion experiments. This latter set encompasses all experiments included in the Model Validation Kit (MVK) package [Olesen 1995]; diffusion experiments under low wind conditions with low emission sources (such as the Idaho and Raaba experiments, where GRAL performed well, while AUSTAL2000 and ADMS completely failed); the Gratkorn experiment, in a mountain valley; the Idaho Falls experiment with a linear emission source; concentration measurements around highways (Elimäki); in a street canyon (Göttinger Strasse, Hanover; Frankfurter Allee, Berlin; Hornsgatan, Stockholm); in a parking lot (Vienna); near a pig farm (Uttenweiler, Germany; Roager, Denmark); near an experimental reactor (EOCR, Idaho); around gas compressors (Texas, Kansas); and at the exits of road tunnels (Ninomiya, Hitachi, Enrei, Japan). In all but one case, the GRAL model provided the best results among the compared models, or results that were very close to the best. Additionally, the model was tested for odour nuisance forecasting around in four experiments in livestock farms, and yielded positive results.

As an example of GRAL’s applications, it is worth mentioning the recent study by Petrov et al. [2024], who conducted high-resolution (2-m and 30-m) sensitivity analysis of modelled NOx emission distribution in Sofia under various meteorological conditions, represented by predefined wind and atmospheric stability. This research highlighted GRAL’s ability to account for wind flow dynamics and produce spatially-detailed pollutant concentration maps. Although challenges were encountered in replicating dispersion effects in complex environments, the study underscored GRAL’s computational efficiency for long-term simulations and its suitability for analysing emission scenarios to improve air quality or assess population exposure.

Similarly, Berchet et al. [2017] evaluated high-resolution (10 m) GRAL simulations of NOx concentrations over Zürich, Switzerland. Their study demonstrated the compliance of GRAL’s overall performance with the objective criteria of the European Commission expert panel FAIRMODE (Forum for Air-quality Modelling in Europe) [Thunis et al. 2012, Pernigotti et al. 2013]. Despite a persistent tendency to overestimate concentrations, the model was able to successfully capture averaged spatial distributions in different temporal scales, with an emphasis on representing properly diurnal cycles.

Another example is provided by Patiño et al. [2024], who examined the suitability of dispersion models with different levels of complexity for use in urban planning. They found that GRAL/GRAMM, along with a more advanced large eddy simulation (LES) model, PALM, could effectively reflect wind-flow patterns and concentration distributions. Despite limitations related to emission representation and sensitivity to meteorological inaccuracies (present in both simulations), GRAL/GRAMM was recognized as a valuable component of the urban air quality assessment framework. Although it has a simplified structure compared to the LES model PALM, it demonstrated satisfactory performance while requiring much fewer computational resources.

Still another important comparison is presented by Oettl et al. [2001], who evaluate the performance of a Gaussian plume model, CAR-FMI, and GRAL in predicting NO_x concentrations under low wind speed conditions. Model results were compared with the field measurements from a major road in Elimäki, Finland. The performance of the CAR-FMI model declined when the wind direction became nearly parallel to the road, as well as under the lowest wind speed conditions. In contrast, the GRAL model exhibited greater stability across different wind speeds and directions, providing more accurate results.

GRAL/GRAMM was also compared to CALPUFF in a complex urban environment [Ward, Rollings 2021; Rollings 2022] and showed superior results due to its ability to account for terrain obstacles such as trees, buildings and fences. On the other hand, modeling the impact of microscale terrain features poses problems, primarily with flow modeling. The practical application of full hydrodynamic models at this scale is limited by the immense computational power required.

From the perspective of using the model as a regulatory or reference tool, the terms of its availability under the GNU Public License ver. 3 are of fundamental importance. This also applies to its graphical interface, which, in the case of other models, is most often developed and provided as closed binary code for a fee by third-party companies. The graphical interface, as seen in Fig. 1(a), enables for convenient entry and modification of parameters related to the model settings, such as resolution; number and placement of modelling levels; emission source characteristics; meteorological data; and topography. Within the meteorological data section, users can produce diagnostic graphics, such as wind roses and wind speed histograms, as seen in Fig. 1(b). The results of pollution dispersion simulations can be displayed as concentration isolines superimposed on a topographic map or satellite image. Contaminant concentration maps can also be represented at different predefined heights [Fig. 2].

Figure 1.

(a) The graphical interface of the GRAL model. (b) A wind rose automatically generated in the GRAL model interface window

Figure 2.

Example of modelled concentration maps from the GRAL model of a pollutant emitted from a point source, superimposed on a satellite map: (a) at a height of 3 m and (b) at a height of 6 m.

Another model worth noting is LAPMOD [Bellasio et al. 2017; Bellasio et al. 2018], developed by the Italian company Environware. ⁽⁶⁾ LAPMOD is a complete modelling system, since it is fully coupled with the diagnostic meteorological model CALMET; it works in conjunction with the LAPEMI emission modelling subsystem and the LAPOST result processing and analysis subsystem. It also has another meteorological preprocessor, LAPMET, that can use US-EPA AERMOD meteorological input files (surface and profile). The model has been validated using the standard data set included in the MVK package. The company provides the source code and documentation for these programs free of charge; however, the licence terms restrict their redistribution.

The model offers several alternative parameterizations for buoyancy effects in plumes of hot gases, as well as multiple alternative kernels for the integral operators that can be used to calculate concentrations based on the set of pseudoparticle coordinates; the latter can be highly useful in development work. Another useful feature of the LAPMOD model is the availability of a module designed for solving inverse problems, i.e., source attribution. Analysis demonstrated by ul Haq et al. [2019] highlights the ability of the LAPMOD model to properly simulate the trajectory of released particles within a complex terrain. This feature can be useful for simulating significant extreme events, as shown in Yang et al. [2024], where authors analysed local-scale dispersion of atmospheric Caesium-137 following the Fukushima Daiichi nuclear accident using the CALMET-LAPMOD coupled modelling system.

The final model worth mentioning here is FLEXPART ⁽⁷⁾ [Stohl et al. 1998; Stohl et al. 2005; Pisso et al. 2019; Bakels et al. 2024]. This model has primarily been used for simulating the dispersion of pollutants on a large spatial scale — i.e., continental. Two noteworthy aspects are the extensive list of over 200 scientific publications related to the model that are included in the Web of Science database, and the availability of the program’s source code under the GNU Public License. Another valuable feature is its ability to calculate source-receptor matrices, which are essential for source apportionment. Since it was developed primarily for research rather than regulatory applications, it lacks a graphical user interface, which is often expected by beginners or occasional users. Nevertheless, it is a robust, state-of-the-art model with significant potential for adaptation to practical applications, such as validation of national-scale emission inventories, as presented in Henne et al. [2016].

FLEXPART was also used to investigate the transport pathways and source contributions to black carbon particles in the Arctic region [Zhu et al. 2020], where it proved to be a good emission source detector, attributing near-surface concentrations of black carbon over the Arctic to local anthropogenic emissions in the northern regions of Russia, as well as source regions in nearby Nordic countries. The model has been also successfully used in an integrated forecasting system in the recent study by Foreback et al. [2024], where it simulated back trajectories during heavy pollution events (under stagnant, low-wind conditions and abrupt synoptic weather changes) in Beijing.

Before concluding, it is worth adding a historical note. An LPD model was also developed in Poland during the late 1980s to provide meteorological support for the Żarnowiec Nuclear Power Plant [Uliasz 1990a; Uliasz 1990b]. That work was carried out using an IBM-PC/AT class personal computer, equipped with a 16-bit Intel 286 processor and 1 MB of RAM—four to five orders of magnitude below the capabilities of today’s laptops. Later on, it was also used to analyse air pollution levels in the Izera Mountains region (the so-called “Black Triangle”) [Uliasz et al. 1994], and continued to be developed for several more years [Uliasz 1993; Uliasz et al. 1996; Uliasz, Olendrzyński 1997]. Unfortunately, work on the model was discontinued before the end of the 20th century.