The Lagrangian Atmospheric Radionuclide Transport Model (ARTM) – sensitivity studies and evaluation using airborne measurements of power plant emissions

The Atmospheric Radionuclide Transport Model (ARTM) operates at the meso-γ scale and simulates the dispersion of radionuclides originating from nuclear facilities under routine operation within the planetary boundary layer. This study presents the extension and validation of this Lagrangian particle dispersion model and consists of three parts: (i) a sensitivity study that aims to assess the impact of key input parameters on the simulation results, (ii) the evaluation of the mixing properties of five different turbulence models using the well-mixed criterion, and (iii) a comparison of model results to airborne observations of carbon dioxide (CO₂) emissions from a power plant and the evaluation of related uncertainties. In the sensitivity study, we analyse the effects of the stability class, roughness length, zero-plane displacement factor, and source height on the three-dimensional plume extent as well as the distance between the source and maximum concentration at the ground. The results show that the stability class is the most sensitive input parameter as expected. The five turbulence models are the default turbulence models of ARTM 2.8.0 and ARTM 3.0.0, one alternative built-in turbulence model of ARTM, and two further turbulence models implemented for this study. The well-mixed condition tests showed that all five turbulence models are able to preserve an initially well-mixed atmospheric boundary layer reasonably well. The models deviate only 6 % from the expected uniform concentration below 80 % of the mixing layer height, except for the default turbulence model of ARTM 3.0.0 with deviations of up to 18 %. CO₂ observations along a flight path in the vicinity of the lignite power plant Bełchatów, Poland, measured by the Deutsches Zentrum für Luft- und Raumfahrt (DLR) Cessna aircraft during the Carbon Dioxide and Methane Mission (CoMet) campaign in 2018 allowed for evaluation of model performance for the different turbulence models under unstable boundary layer conditions. All simulated mixing ratios are of the same order of magnitude as the airborne in situ data. An extensive uncertainty analysis using probability distribution functions, statistical tests, and direct spatio-temporal comparisons of measurements and model results help to quantify the model uncertainties. With the default turbulence setups of ARTM versions 2.8.0 and 3.0.0, the plume widths are underestimated by up to 50 %, resulting in a strong overestimation of the maximum plume CO₂ mixing ratios. The comparison of the three alternative turbulence models shows good agreement of the peak plume CO₂ concentrations, the CO₂ distribution within the plumes, and the plume width, with a 30 % deviation in the peak CO₂ concentration and a less than 25 % deviation in the measured CO₂ plume width. Uncertainties in the simulations may arise from the different spatial and temporal resolutions of simulations and measurements in addition to the turbulence parametrisation and boundary conditions. The results of this work may help to improve the accurate representation of real plumes in very unstable atmospheric conditions through the selection of distinct turbulence models. Further comparisons at different stability regimes are required for a final assessment of model uncertainties.

1 Introduction

Atmospheric dispersion models (ADMs) are widely used by the scientific community and authorities. They are applied to a variety of problems, such as the study of the impact of pollutant emissions on air quality (Gariazzo et al., 2007; Stohl et al., 2007; Berchet et al., 2017; Lonati et al., 2022; Shupe et al., 2022) or the dispersion of radioactive discharges to the air (Chino et al., 2011; Connan et al., 2013; Draxler et al., 2015; Arnold et al., 2015), and they can operate at the full meteorological scale, ranging from the microscale (shorter than 1 km) to the mesoscale (1 km to 1000  km) up to the synoptic or global scale (larger than 1000  km).

The Atmospheric Radionuclide Transport Model (ARTM), analysed in this study, belongs to the class of models operating at the micro-β to meso-γ scales (approx. 0.5 to 20 km). It is a Lagrangian particle dispersion model (LPDM) designed for the dispersion of radionuclides from nuclear facilities under routine operation in the planetary boundary layer (PBL).

However, any ADM has to demonstrate its applicability to the system of study. The most important method to confirm this is validation (Kleijnen, 1995; Schlesinger et al., 1979). This includes (i) sensitivity analysis, which relates the model's response to variations in the model's input parameters, and (ii) a comparison with observations, revealing whether a model is an accurate representation of the system and whether simulation results are, to a certain degree, in agreement with observations (Kleijnen, 1995; Rao, 2005).

Concerning ARTM, Hettrich (2017) performed a sensitivity study to analyse the effect of input parameters (e.g. emission rate, source geometry, stability class, and particle number) on concentrations at selected locations near the ground. Hanfland et al. (2022) provided an overview of the physical basis and mathematical formulations of the model and presented a qualitative description of the influences of different input parameters on three-dimensional plume characteristics for a general simulation setup. However, both lack in quantifying sensitivities.

Here we expand these former studies with a more systematic and quantitative sensitivity analysis. Different sensitivity coefficients are calculated, which describe the dependence of the simulation output on the input parameter stability class (SC), roughness length (z₀), zero-plane displacement factor (d), and source height (h_s) within the whole simulated PBL and rank them according to their effects on the model output.

Publications presenting comparisons of ARTM's mixing ratio simulation results with measurements are rare. Hettrich (2017) compared ARTM simulation results with measurements at a few selected locations near the surface, showing discrepancies which are related to complex orography or local thermal-induced winds that are not covered by ARTM's wind field model TALdia. Martens et al. (2012) studied the influence of a single large building close to the source on the dispersion, showing that at distances larger than 4 km the influence decreases. These comparisons covered only a few atmospheric conditions and were limited to near-surface concentration measurements.

Brunner et al. (2023) presented an intercomparison of six different atmospheric transport models, including ARTM, with airborne in situ and remote sensing carbon dioxide (CO₂) measurements, sampling the exhaust plume of the Bełchatów lignite power plant in Poland under very unstable atmospheric conditions (Fix et al., 2018). The data set comprises a considerable number of plume transects at different distances from the source and heights within the PBL, is characterised by a strong contrast between background and plume CO₂ mixing ratios, and provides a three-dimensional description of the mixing ratio field. The spatial extent of the area covered by measurements is around the maximum domain size ARTM can tackle. In the Brunner et al. (2023) study, ARTM simulations are performed using the default turbulence model of ARTM version 2.8.0, including a workaround for meandering plumes because the simulated plume appeared to be too narrow under very unstable atmospheric conditions.

ARTM's dispersion depends on the turbulence model used (Hanfland et al., 2022), and since the three-dimensional data set facilitates further analysis of the model, we investigate whether the modelled plume could become more realistic by using different turbulence models. Two new turbulence models were implemented based on the ideas of Hanna (1982) and Degrazia et al. (2000) in addition to three built-in models of ARTM. The same workaround for meandering plumes based on the default turbulence model of ARTM 2.8.0, as primarily presented in Brunner et al. (2023), is also included in this study. Furthermore, all five turbulence models are evaluated concerning turbulence model characteristics and mixing efficiency and are compared with measurements.

The structure of this study is as follows: Sect. 2 gives a short description of the atmospheric dispersion model ARTM. Section 3 introduces several different sensitivity analysis methods and presents the sensitivity of typical model simulation outputs to key input parameters. Section 4 presents the five turbulence models and assesses their performance with respect to the well-mixed condition. Section 5 shows the comparison and evaluation of ARTM simulation results for the five turbulence models with the three-dimensional airborne data. Section 6 concludes the results of this paper.

2 ARTM

The Atmospheric Radionuclide Transport Model is an LPDM developed specifically for the dispersion of radioactive emissions from nuclear facilities in an area of typically 10 km × 10 km. Its purpose is to provide annual activity concentration fields in the area around nuclear facilities under routine operation in slightly structured non-urban terrain, which are used in a follow-up step to calculate the additional exposure of the population (Hanfland et al., 2022).

The dispersion model propagates numerical particles representing radioactive tracers in space and time according to wind and turbulence fields obtained by a diagnostic approach. Meteorological data from measurements in the vicinity of the nuclear facilities are used to calculate a mass-conserving diagnostic wind field. The turbulence is obtained by a Markov process, which uses wind speed fluctuations and Lagrangian correlation times as input parameters, both depending on the Obukhov length as a turbulence parameter. This static diagnostic approach employed by ARTM is similar to that of, for example, SWIFT/micro-SWIFT or CALMET (Cox et al., 2005; Moussafir et al., 2004; Scire et al., 1998) but differs from other larger-scale LPDMs, such as FLEXPART, STILT, NAME, or HYSPLIT, which use prognostic meteorological fields from a numerical weather prediction model to drive the particle propagation (Hanfland et al., 2022; Lin et al., 2003; Stohl et al., 2005; Ryall and Maryon, 1998; Draxler and Hess, 1998). The advantage of this static diagnostic approach is its high computational efficiency with typical applications from 10 km to a few hundred kilometres for horizontal extent, depending on the number of meteorological measurement locations and the terrain (Ratto et al., 1994). Since ARTM works with one single meteorological measurement location, its simulation domain can horizontally extend to up to a few tens of kilometres depending on the terrain. The simulation domain is divided into grid cells for which the average activity concentration for the whole simulation period is calculated. A detailed model description is given in Hanfland et al. (2022).

3 Sensitivity study

Sensitivity analysis (SA) is a method to study the model's response to variations in input parameters in a systematic way, and it may answer the following questions. How does the uncertainty in input parameters influence the model output? Which parameters require additional research in order to reduce output uncertainty? Which parameters are the most significant or insignificant for the model's output? And does the model behave as expected when varying a certain input parameter (Hamby, 1994; Frey and Patil, 2002; Rao, 2005; Saltelli et al., 2008)? SA methods are either local or global depending on the sampled input parameter space (Saltelli et al., 2008; Morio, 2011; Zagayevskiy and Deutsch, 2015).

The results of the methods may differ depending on the shape of the input parameter space. Thus, the application of several methods is recommended (Iman and Helton, 1988; Hamby, 1995). In this work, several different local and global SA methods are therefore applied to ARTM 2.8.0 using its default turbulence model (Hanfland et al., 2022). A description of the turbulence model for unstable atmospheric conditions is given in Sect. 4.1. The application of several SA methods provides a comprehensive assessment of the response of ARTM to different input parameters.

3.1 Local sensitivity analysis methods

Local SA focuses on one single point in the input parameter space. The output of a model is represented by $Y=g(X_{\mathrm{1}},\mathrm{\dots },X_k)$, where the random variables X_i with $i=\mathrm{1},\mathrm{\dots },k$ denote the different input parameters. The representations (or values) of X_i are denoted with x_i. The input parameters X_i are varied one at a time, while all the other parameters are held constant at their reference values $x_i^{\text{ref}}$. This local SA approach is similar to estimating the partial derivative $\frac{\partial Y}{\partial X_i}$ and characterises the effect of the input parameter X_i on Y at one reference point $X^{\text{ref}}=(x_{\mathrm{1}}^{\text{ref}},\mathrm{\dots },x_k^{\text{ref}})$ (Morio, 2011).

3.1.1 Sensitivity index

The sensitivity index described by Hoffman and Gardner (1983) uses the parameters at the reference point where each parameter is varied one at a time by their full range. The sensitivity index is calculated as

$$\begin{array}{}\text{(1)}& {\mathrm{SI}}_i={\displaystyle \frac{Y_{i,\text{max}}-Y_{i,\text{min}}}{Y_{i,\text{max}}}},\end{array}$$

where Y_{i,max (min)} indicates the maximum (minimum) output value. The sensitivity index is a value between $\mathrm{0}\le {\mathrm{SI}}_i\le \mathrm{1}$ and gives the fraction of output variation caused by the varied input parameter (Hamby, 1994).

3.1.2 One-at-a-time sensitivity measure

The one-at-a-time sensitivity measure calculates the variation in the model output normalised to the largest output variation ΔY_max that had been observed for the different input parameters. Starting from the default parameter set, the parameters are varied one at a time by a percentage α. The sensitivity coefficient for the input parameter X_i is calculated as

$$\begin{array}{}\text{(2)}& {\mathrm{SM}}_i^{\mathit{\alpha }}={\displaystyle \frac{|Y_{i,+\mathit{\alpha }}-Y_{i,-\mathit{\alpha }}|}{\mathrm{\Delta }{Y}_{\text{max}}}},\end{array}$$

where $\mathrm{\Delta }{Y}_{\text{max}}=max\left(\right|Y_{l,+\mathit{\alpha }}-Y_{l,-\mathit{\alpha }}\left|\right)\forall l\in i$ (Link et al., 2018). In this work the percentages ± 25 % and ± 50 % are used for α. ${\mathrm{SM}}_i^{\mathit{\alpha }}$ is a value between $\mathrm{0}\le {\mathrm{SM}}_i^{\mathit{\alpha }}\le \mathrm{1}$, where unity identifies the input parameter with the biggest effect on the model output Y.

3.2 Global sensitivity analyses

Global SAs sample the whole input parameter space, which leads to a broader representation of the sensitivity compared to local methods but also increases computation time. A general discussion about global SAs can be found in Saltelli et al. (2008).

3.2.1 Sobol′ indices

The variance-based Sobol′ indices use variance decomposition to calculate indices of different orders (Sobol', 1993). Usually, only two key Sobol′ indices are determined: the first-order index S_i and the total-effect index S_Ti.

For the first one, the conditional expected value of the model output $E_{X_{\sim i}}\left(Y\right|X_i)$ with a constant value of X_i and varying values for every other input parameter X_∼i is computed. For different realisations of X_i, $V_{X_i}\left[E_{X_{\sim i}}\right(Y\left|X_i\right)]$ reflects the variance of the model output Y originating from a variation in the input parameter X_i. The first Sobol′ index is then given by

$$\begin{array}{}\text{(3)}& S_i={\displaystyle \frac{V_{X_i}\left[E_{X_{\sim i}}\right(Y\left|X_i\right)]}{V\left(Y\right)}},\end{array}$$

where V(Y) is the unconditional variance of the output where all X_i values are varied. $V_{X_i}\left[E_{X_{\sim i}}\right(Y\left|X_i\right)]$ cannot be larger than V(Y), and thus for the sensitivity coefficient $\mathrm{0}\le S_i\le \mathrm{1}$ is valid. This index is called the first-order sensitivity index as it does not take higher-order effects (i.e. interactions between different input parameters) into account (Saltelli et al., 2008).

The second index considered here is the total-effect index. It takes higher-order terms into account, which might be important depending on the model. The total effect is calculated as

$$\begin{array}{}\text{(4)}& S_{\mathrm{T}i}=\mathrm{1}-{\displaystyle \frac{V_{X_{\sim i}}\left[E_{X_i}\right(Y\left|X_{\sim i}\right)]}{V\left(Y\right)}},\end{array}$$

where $V_{X_{\sim i}}\left[E_{X_i}\right(Y\left|X_{\sim i}\right)]=V_{X_{\sim i}}[E_{X_i}\left(Y\right|X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{\dots },X_{i-\mathrm{1}}$, X_i+1, …, X_k)]≤V(Y) is the total variance of all input parameters except X_i. As the first-order index, this total-effect index is a value between 0 and unity, where a value of 0 indicates no influence of X_i on the output Y, while unity indicates a strong influence (Saltelli et al., 2008). A comprehensive description of the method is given by Saltelli et al. (2008).

For the analysis presented here, the Python library SALib (Herman and Usher, 2017) is used for the quasi-random sampling with low discrepancy in the input parameter space, after Joe and Kuo (2008), as well as for the calculation of the Sobol′ indices. It furthermore allows for the estimation of the 95 % confidence intervals (Herman and Usher, 2021).

3.2.2 δ method

In comparison with the Sobol′ indices, the δ method takes the complete density distribution of the model output into account, which ensures the conservation of all the information of the output density distribution (Borgonovo, 2007). The probability density function of X_i is denoted $f_{X_i}\left(x_i\right)$. The sensitivity coefficient δ_i for the input parameter X_i is calculated using the marginal density distribution of the input parameter $f_{X_i}\left(x_i\right)$ and the difference between the unconditional and the conditional density functions f_Y(y) and $f_{Y|X_i}\left(y\right)$ of the model output with fixed representation X_i=x_i as (Borgonovo, 2007)

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}_i={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\int f_{X_i}\left(x_i\right)[\int |f_Y\left(y\right)-f_{Y|X_i}\left(y\right)\left|\mathrm{d}y\right]\mathrm{d}{x}_i.\end{array}$$

δ_i represents the total effect of an input parameter X_i on Y. It can take a value between 0 and unity ($\mathrm{0}\le {\mathit{\delta }}_i\le \mathrm{1}$), where 0 means that the output is independent of X_i (Plischke et al., 2013). The same library (SALib; Herman and Usher, 2017) was used to apply the δ method and estimate the 95 % confidence interval.

3.3 Model setup for sensitivity analyses

For the SA a simple model setup with a single source and constant wind was chosen. Figure 1 illustrates the simulation domain with the plume sampled at a height of 20.5 m above ground level (a.g.l.). The simulation domain has an extent of 10 km × 1.5 km × 1.5 km in the x, y, and z directions. The x direction is defined in a west–east orientation and the y direction in a north–south orientation. The domain is divided into grid cells with a horizontal resolution of 50 m. Vertically, the 1.5 km high simulation domain is divided into 19 levels of varying thickness gradually increasing from the lowest layer (3 m thick) to the top simulation layer (300 m thick). The level thicknesses are shown in Table S1 in the Supplement. The point source is located at the coordinates (x = 25 m; y = 25 m). The vertical position is varied during the SAs. A constant westerly wind (270°) was used for the entire simulation period of 24 h with a velocity of 1 m s⁻¹ at 10 m height. In order to focus on the evolving dispersion pattern, the topography is assumed to be a flat surface. The gaseous krypton isotope ⁸⁵Kr with a decay constant of λ_decay = 2.05 × 10⁻⁹ s⁻¹ was used as a tracer. This results in a decay of less than 0.02 % within the simulation period and can therefore be neglected. The emission source is represented as a source with a constant activity rate of 1 Bq s⁻¹.

Figure 1 The x–y plane of the simulation domain for the sensitivity analyses with the activity concentration at a height of 20.5 m. The generic plume is simulated for SC (neutral), z₀ (0.5 m), d (6), and h_s (20 m) with a wind speed of 1 m s⁻¹ from the west at 10 m height and an emission rate of 1 Bq s⁻¹. The activity concentration distribution in Bq m⁻³ is given in a logarithmic scale.

Table 1 Input parameters and their values and ranges. The default parameters for local SAs are marked with ^∗.

¹ For the global SAs the given values are sampled. ² For the global SAs the values are sampled continuously within the range. ³ For the global SAs the values are sampled with a resolution of 1 m within the range.

The input parameters stability class SC, roughness length z₀, zero-plane displacement factor d, and source height h_s are assumed to be the key parameters of ARTM. The zero-plane displacement d₀ depends on d as $d_{\mathrm{0}}=d\cdot z_{\mathrm{0}}$. When using d instead of d₀ the input parameters for the SA are independent. This allows for the analysis of the unbiased effects of input parameter variations. These parameters and their ranges are summarised in Table 1. Only a discrete set of six SCs is allowed in ARTM. The range of z₀ values allowed by ARTM is limited to the roughness lengths that correspond to typical land covers in the vicinity of German nuclear facilities. German authorities recommend a value of d=6 for the zero-plane displacement factor (TA Luft, 2002; VDI 3783 part 8, 2017). The range for the variation in d is centred on this value and limited to forest canopy heights typical of mixed forest (Lang et al., 2022). In Table 1, the parameter values for the reference point for the local SAs are marked with a * symbol. For global sensitivity analyses, the whole parameter range is sampled.

The target quantities of the SAs are two important characteristics of the plume: (i) the plume volume, which is a measure of the tracer dispersion and is closely linked to the maximum mixing ratio, and (ii) the distance between the location of the maximum activity concentration and the source position at ground level, which is of special interest for radiation exposure assessment.

Table 2 Sensitivity coefficients of local and global sensitivity analyses for the plume volume and for the distance between the source and the maximum concentration at the ground level. For the Sobol′ indices (S_i and S_Ti) and the δ method (δ_i) 95 % confidence intervals are given as well. Coefficients with very large relative confidence intervals are marked with ^∗, and coefficients of one method, which cannot be distinguished within their confidence intervals, are marked with ^†.

Table 3 Ranking of the influence of the input parameters on the plume volume and on the distance between the source and the maximum concentration at ground level for local and global sensitivity analysis methods.

3.4 Results of the sensitivity analyses

The results of the calculations of the local and global sensitivity coefficients are summarised in Table 2. Concerning the plume volume, all SA methods compute the highest SA coefficients for the stability class. Although less prominent, this is also observable for the distance between the source and maximum concentration at the ground level, except for sensitivity index SI_i.

For ${\mathrm{SM}}_{z_{\mathrm{0}}}^{\mathrm{25}}$ no value can be calculated because a variation of ± 25 % from the reference roughness length value does not lead to a change in the categorial z₀ value. For the other parameters, the two different ranges of variation (α=25 and α=50) for ${\mathrm{SM}}_i^{\mathit{\alpha }}$ provide valuable additional information. For example, it can be seen from Table 2 that the deviations between the coefficients of ${\mathrm{SM}}_i^{\mathrm{25}}$ and ${\mathrm{SM}}_i^{\mathrm{50}}$ are small for the plume volume, while they are large for the distance between the source and the maximum concentration. The influence of the input parameters seems to be rather linear for the plume volume but highly non-linear for the distance between the source and maximum concentration at the ground.

For the global SA methods, both target quantities show a distinct importance not only of first-order (direct influence of one single input parameter) but also of higher-order (including interactions of two or more input parameters) effects. A small difference between S_i and S_Ti shows a large first-order effect, as can be seen for the plume volume. On the contrary, a large difference reveals a small first-order effect rather than a higher-order effect, as can be seen for the distance between the source and maximum concentration at ground level. This agrees with the conclusions that can be drawn from the δ_i coefficients. The sum of the sensitivity coefficients for the plume volume ${\sum }_i{\mathit{\delta }}_i=\mathrm{1.05}\pm \mathrm{0.01}\cong \mathrm{1}$ indicates that the effects of variation in the input parameters on variation in the plume volume are separable; i.e. interactions between input parameters play a minor role. For the distance between the source and maximum concentration, the sum of the sensitivity coefficients ${\sum }_i{\mathit{\delta }}_i=\mathrm{0.690}\pm \mathrm{0.008}\textcolor{red}{}\mathrm{1}$ indicates the important role of cross-interactions between the input parameters (Borgonovo, 2007). Contrary to the findings of the Sobol′ indices that some input parameters have negligible influences, the δ method suggests that the output characteristics are sensitive to all parameters. This difference could be due to the different amounts of information processed by the two methods. While the Sobol′ indices compare conditional and unconditional variances of the output distribution, the δ method takes the entire output distribution into account.

Some of the global SA coefficients have very large relative confidence intervals and cannot be distinguished from 0 (marked with ^∗). Others cannot be distinguished from each other within their confidence intervals (marked with ^†). Increasing the sample size of 24 576 further would be necessary to get smaller confidence intervals, but this would also increase the computation time (Herman and Usher, 2021).

Based on the coefficients from Table 2, the input parameters were ranked according to their importance as summarised in Table 3. The rankings obtained for the individual SA methods differ not only for the two target quantities but also between different methods. The overall ranking, which is simply computed as the sum (Σ) of the different methods, is provided in the second-last column.

The most unambiguous result is that all SA methods show the plume volume to be the most sensitive to the SC. The ranks for the other input parameters, in contrast, are not the same. At this point we want to mention that the ranking of ${\mathrm{SM}}_i^{\mathrm{25}}$ given in Table 3 is the average of all possible rankings for this method when taking into account that there is no coefficient for ${\mathrm{SM}}_{z_{\mathrm{0}}}^{\mathrm{25}}$. The rankings of the remaining local SA methods, SI_i and ${\mathrm{SM}}_i^{\mathrm{50}}$, are in agreement with each other for the plume volume, while the rankings of the global SA methods disagree. Compared to the rankings for the plume volume, those for the distance between the source and maximum concentration at the ground level are less uniform.

The overall rankings for both target quantities differ from each other, which emphasises that different target quantities are not necessarily sensitive to the same input parameters. Both target quantities are the most sensitive to the SC, which is thus a potential source of high uncertainty. The source height h_s has little influence on the plume volume, but it is the second most important parameter for the distance between the source and maximum concentration at the ground level.

4 Turbulence models for unstable conditions and their performance in the well-mixed test

In LPDMs the turbulent motion is described via a Markov chain approach in the form of a Langevin equation (Lin and Gerbig, 2013). ARTM uses the wind speed fluctuation σ and Lagrangian correlation time T_L as input parameters for this Markov chain approach (Hanfland et al., 2022). The turbulence model specifies these two parameters and thus influences the tracer dispersion and hence the simulated mixing ratio field (Katharopoulos et al., 2022).

Figure 2 Vertical profiles of the model characteristics of the five turbulence models, ARTM2, ARTM3, PRFMOD, MODHANNA, and DEGRAZIA, under very unstable atmospheric conditions. Wind speed fluctuation σ: (a) along-wind direction σ_x, (c) crosswind direction σ_y, and (e) vertical direction σ_z against the normalised height (normalised to the boundary layer height). The corresponding Lagrangian correlation time T_L is shown in (b), (d), and (f). The turbulent kinetic energy per unit mass TKE is shown in (g).

4.1 Description of turbulence models

The turbulence model implemented in ARTM 2.8.0 as the default model is not widely used in the scientific community. Besides this, it has been reported by Janicke and Janicke (2011) that it sometimes underestimates plume dispersion. Therefore, they introduced a modified turbulence model, leading to stronger dispersion, which can optionally be activated in ARTM. In 2022, the new version (3.0.0) of ARTM was released. It implements a new turbulence model according to the Association of German Engineers (VDI) guideline VDI 3783 part 8 (2017). All three models deviate from the model suggested by Hanna (1982), which is quite widely used and thoroughly tested against tracer release experiments. However, in this model the turbulence may abruptly change between SCs. To overcome this issue of discontinuity, Degrazia et al. (2000) proposed a continuous description of the turbulence throughout all atmospheric conditions. Since the measurement data set used for the comparison of simulations and observations was collected under unstable atmospheric conditions, the following analyses focus on unstable stratification. The wind speed fluctuation σ and the Lagrangian correlation timescale T_L of the five turbulence models are presented in Eqs. (6)–(26) for unstable stratification, and their profiles are displayed in Fig. 2. For the following quantities we define the x components along the average horizontal wind direction, the y components perpendicular to it in the horizontal plane, and the z components in the vertical direction. Although the zero-plane displacement is used in ARTM (GRS, 2015) to displace the wind profile vertically to account for the influence of obstacles, for the sake of simplicity it is not included in the following equations.

The first model, the default boundary layer model (BLM) of ARTM 2.8.0, was initially suggested by Kerschgens et al. (2000) and is based on the work of Lenschow et al. (1980), Panofsky et al. (1977), Hicks (1985), and Gryning et al. (1987). It describes profiles for the wind speed fluctuations as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\sigma }}_x=\mathrm{2.4}\cdot u_{\ast }{\left(\mathrm{1}+\mathrm{0.01486}{\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\cdot \exp \left({\displaystyle \frac{-z}{h_{\mathrm{m}}}}\right),\text{(7)}& {\displaystyle }{\mathit{\sigma }}_y=\mathrm{1.8}\cdot u_{\ast }{\left(\mathrm{1}+\mathrm{0.03522}{\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\cdot \exp \left({\displaystyle \frac{-z}{h_{\mathrm{m}}}}\right),\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_z=\mathrm{1.3}\cdot u_{\ast }{\left[{\left(\mathrm{1}-\mathrm{0.8}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\mathrm{3}}\cdot {\displaystyle \frac{-z}{\mathit{\kappa }L}}+\exp {\left({\displaystyle \frac{-z}{h_{\mathrm{m}}}}\right)}^{\mathrm{3}}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

where u_∗ is the friction velocity, h_m is the mixing layer height, κ=0.4 is the von Kármán constant, L is the Obukhov length, and z is the height above ground level (VDI 3783 part 8, 2002; Hanfland et al., 2022). This model is called ARTM2 in the following.

The second turbulence model available in ARTM is based on ARTM2 with a modification in the exponents and in the prefactor of the crosswind component as follows (Janicke and Janicke, 2011):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{\sigma }}_x=\mathrm{2.4}\cdot u_{\ast }{\left(\mathrm{1}+\mathrm{0.01486}{\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\cdot \exp \left({\displaystyle \frac{-\mathrm{0.3}\cdot z}{h_{\mathrm{m}}}}\right),\text{(10)}& {\displaystyle }{\mathit{\sigma }}_y=\mathrm{2.0}\cdot u_{\ast }{\left(\mathrm{1}+\mathrm{0.03522}{\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\cdot \exp \left({\displaystyle \frac{-\mathrm{0.3}\cdot z}{h_{\mathrm{m}}}}\right),\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_z=\mathrm{1.3}\cdot u_{\ast }{\left[{\left(\mathrm{1}-\mathrm{0.8}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\mathrm{3}}\cdot {\displaystyle \frac{-z}{\mathit{\kappa }L}}+\exp {\left({\displaystyle \frac{-\mathrm{0.3}\cdot z}{h_{\mathrm{m}}}}\right)}^{\mathrm{3}}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

This model leads to wider plumes and is called PRFMOD in the following.

In addition to the two previous models, we added a turbulence model to ARTM based on ARTM2 modified with formulations used in other ADMs (Stohl et al., 2005). This model uses σ_z from ARTM2 given in Eq. (8), but the horizontal wind speed fluctuations

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_x={\mathit{\sigma }}_y=u_{\ast }{\left(\mathrm{12}+{\displaystyle \frac{h_{\mathrm{m}}}{\mathrm{2}\left|L\right|}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\end{array}$$

are equal to the equations suggested by Hanna (1982). The aim of this modification is to increase the turbulent kinetic energy and to analyse the effect of the horizontal wind speed fluctuations on the dispersion. In the following, this model is called MODHANNA.

The Lagrangian correlation times of the three models are given according to Kolmogorov's theory as (Luhar and Britter, 1989)

$$\begin{array}{}\text{(13)}& T_{\mathrm{L}i}={\displaystyle \frac{\mathrm{2}\cdot {\mathit{\sigma }}_i^{\mathrm{2}}}{{\mathrm{C}}_{\mathrm{0}}\cdot \mathit{\eta }}},\end{array}$$

with the Kolmogorov constant C₀=5.7 and the dissipation rate of the turbulent kinetic energy

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathit{\eta }=& max\mathit{\{}{\displaystyle \frac{u_{\ast }^{\mathrm{3}}}{\mathit{\kappa }z}}\left[{\left(\mathrm{1}-{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\mathrm{2}}+{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right]\\ & +{\displaystyle \frac{-u_{\ast }^{\mathrm{3}}}{\mathit{\kappa }L}}\left[\mathrm{1.5}-\mathrm{1.3}{\left({\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\right],{\displaystyle \frac{u_{\ast }^{\mathrm{3}}}{\mathit{\kappa }z}}\mathit{\}}.\end{array}\end{array}$$

The fourth model is the default model of the new version of ARTM (3.0.0) with the wind speed fluctuations given as (VDI 3783 part 8, 2017)

$$\begin{array}{}\text{(15)}& {\displaystyle }{\mathit{\sigma }}_x=\mathrm{2.4}\cdot u_{\ast }{\left[\mathrm{1}+\mathrm{0.01486}{\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\cdot \exp \left(-\mathrm{0.9}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{3}}},\text{(16)}& {\displaystyle }{\mathit{\sigma }}_x=\mathrm{2.0}\cdot u_{\ast }{\left[\mathrm{1}+\mathrm{0.02568}{\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\cdot \exp \left(-\mathrm{0.9}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

and

$$\begin{array}{}\text{(17)}& {\mathit{\sigma }}_z=\mathrm{1.3}\cdot u_{\ast }{\left[{\left(\mathrm{1}-\mathrm{0.8}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\mathrm{3}}\cdot {\displaystyle \frac{-z}{\mathit{\kappa }L}}+\exp {\left(-\mathrm{0.9}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\mathrm{3}}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

The Lagrangian correlation timescales are calculated via the turbulent diffusion coefficient K_i as

$$\begin{array}{}\text{(18)}& T_{\mathrm{L}i}={\displaystyle \frac{K_i}{{\mathit{\sigma }}_i^{\mathrm{2}}}},\end{array}$$

with

$$\begin{array}{}\text{(19)}& K_j=\mathrm{0.9}{\displaystyle \frac{u\left(z\right)\cdot h_{\mathrm{m}}}{\mathrm{100}\cdot u_{\ast }}}{\mathit{\sigma }}_j\end{array}$$

for the horizontal component j and

$$\begin{array}{}\text{(20)}& K_z=\mathit{\kappa }{u}_{\ast }z{\left[{\left(\mathrm{1}-\mathrm{0.8}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)}^{\mathrm{4}}{\displaystyle \frac{\mathrm{9}z}{-L}}+\exp \left(-\mathrm{3.6}{\displaystyle \frac{z}{h_{\mathrm{m}}}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}\end{array}$$

for the vertical component (VDI 3783 part 8, 2017). This model is called ARTM3 in the following.

We implemented a fifth model, which, in contrast to the previous four turbulence models that are based on similarity theory, is based on the spectral distribution of the turbulent kinetic energy of the boundary layer and was presented by Degrazia et al. (2000). For very unstable boundary conditions the wind speed fluctuations are given as

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathit{\sigma }}_x=\mathrm{0.53}\cdot u_{\ast }{\left({\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\text{(22)}& {\displaystyle }{\mathit{\sigma }}_y=\mathrm{0.61}\cdot u_{\ast }{\left({\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

and

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\mathit{\sigma }}_z=& \mathrm{0.54}\cdot u_{\ast }{\left({\displaystyle \frac{-h_{\mathrm{m}}}{\mathit{\kappa }L}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ & \times {\left\{\mathrm{1.8}\left[\mathrm{1}-\exp \left({\displaystyle \frac{-\mathrm{4}z}{h_{\mathrm{m}}}}\right)-\mathrm{0.0003}\cdot \exp \left({\displaystyle \frac{\mathrm{8}z}{h_{\mathrm{m}}}}\right)\right]\right\}}^{\frac{\mathrm{1}}{\mathrm{3}}}\end{array}\end{array}$$

with the Lagrangian correlation time

$$\begin{array}{}\text{(24)}& T_{\mathrm{L}i}={\displaystyle \frac{l_i}{{\mathit{\sigma }}_i}},\end{array}$$

where l_i is the Lagrangian correlation length given as

$$\begin{array}{}\text{(25)}& l_{\mathrm{x}}=l_{\mathrm{y}}=\mathrm{0.21}\cdot h_{\mathrm{m}}{\left(\mathrm{0.01}{\displaystyle \frac{h_{\mathrm{m}}}{-L}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\end{array}$$

and

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{l}_{\mathrm{z}}=& \mathrm{0.14}\cdot h_{\mathrm{m}}{\left(\mathrm{0.01}{\displaystyle \frac{h_{\mathrm{m}}}{-L}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\\ & \times \left\{\mathrm{1.8}\left[\mathrm{1}-\exp \left({\displaystyle \frac{-\mathrm{4}z}{h_{\mathrm{m}}}}\right)-\mathrm{0.0003}\cdot \exp \left({\displaystyle \frac{\mathrm{8}z}{h_{\mathrm{m}}}}\right)\right]\right\}.\end{array}\end{array}$$

In this work this turbulence model is denoted as DEGRAZIA.

The turbulent kinetic energy per unit mass is determined as (Stull, 1988)

$$\begin{array}{}\text{(27)}& \mathrm{TKE}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\sigma }}_x^{\mathrm{2}}+{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_z^{\mathrm{2}}\right).\end{array}$$

4.2 Evaluation of turbulent mixing

The degree of preservation of well-mixed conditions is a key quality indicator for any LPDM, similar to the preservation of mass in an Eulerian model. It tests whether an initially uniform distribution of a tracer in an incompressible flow remains uniform as postulated by the second law of thermodynamics (Sawford, 1986; Thomson, 1987; Lin and Gerbig, 2013; Bahlali et al., 2020). Exactly fulfilling this criterion is challenging, but it is important to quantify the degree of deviation from this ideal behaviour to judge the magnitude of systematic model biases and whether these biases are acceptable. In this work the tests of the turbulence models are limited to the case of very unstable atmospheric conditions because the observations that are used for the comparison were collected under very unstable atmospheric conditions, too.

The well-mixed condition test can characterise the vertical mixing homogeneity of a model. For these tests, simulation domains with periodic horizontal boundaries and reflecting vertical boundaries are used. This virtually expands the simulation domain to infinite extent and prevents the simulation from losing tracer mass. The whole simulation domain serves as a volume source where 115 200 numerical particles are inserted uniformly within the first simulation hour because ARTM does not take into account the changing density with height. The domain size is 2000 m × 2000 m × 1100 m in x, y, and z directions, with a horizontal (vertical) resolution of 200 m (25 m). The vertical extent of the domain is equal to the assumed mixing depth. A temporally constant wind profile for unstable atmospheric conditions, as described in Hanfland et al. (2022), with a wind speed of 2.3 m s⁻¹ at 10 m height and a direction of 270 ^∘ (westerly) is used. The selected wind speed originates from measurement sites under very unstable stratification conditions in Germany. For the evaluation, the hourly mean concentration and its standard deviation were derived for each vertical level.

Figure 3 Profiles of the concentration normalised to the mean concentration $c{\stackrel{\mathrm{‾}}{c}}^{-\mathrm{1}}$ (a–e) of the different turbulence models, ARTM2, ARTM3, PRFMOD, MODHANNA, and DEGRAZIA, after 1 (red lines) and 2 h (dashed–dotted blue lines) for periodic lateral simulation domain boundaries and reflecting bottom and top boundaries. In (b) the x-axis scale changes at $c{\stackrel{\mathrm{‾}}{c}}^{-\mathrm{1}}=\mathrm{0.9}$. (f) Time series of the normalised concentration at normalised height $z{h}_{\mathrm{m}}^{-\mathrm{1}}\approx \mathrm{0.3}$ for the ARTM2 model, which is indicated by the dashed horizontal line in (a). The x-axis scale changes at 10 h.

The concentration profiles of the different turbulence models for very unstable PBL conditions are shown in Fig. 3. The concentrations of the state of mixing after 1 (red line) and 2 h (dashed blue line) are shown. Concentration values are normalised to the mean concentration ($c{\stackrel{\mathrm{‾}}{c}}^{-\mathrm{1}}$), and the height is normalised to the mixing depth ($z{h}_{\mathrm{m}}^{-\mathrm{1}}$). We used the same initial numerical particle distribution for all turbulence models to eliminate possible differences arising from different initial distributions.

The concentration profiles after 1 h differ from the uniform distribution $c{\stackrel{\mathrm{‾}}{c}}^{-\mathrm{1}}=\mathrm{1}$. This indicates a certain degree of segregation of the numerical particles, but most deviations are less than 5 % (vertical dashed lines). The largest deviations can be found at the top of the PBL for the ARTM3 model (> 35 %) and the DEGRAZIA model (> 15 %). The profiles of the ARTM2 and MODHANNA turbulence models are very similar since they both contain the same vertical turbulence parametrisation. The PRFMOD turbulence model differs slightly from the ARTM2 model due to modifications described in Eq. (11). The profile of ARTM3 shows trends of dilution and accumulation similar to ARTM2, PRFMOD, and MODHANNA but magnified in its extent. The profile of the DEGRAZIA turbulence model shows a different shape because of its different formulation of the turbulence parameters (see Eqs. 23, 24, and 26).

By t = 2 h, the dilution of the concentration has further increased at the bottom and top of the PBL, and the accumulation at $z{h}_{\mathrm{m}}^{-\mathrm{1}}\approx \mathrm{0.3}$ (horizontal dashed black line) has further increased partly beyond 5 % but well below 10 % deviation for ARTM2, PRFMOD, MODHANNA, and DEGRAZIA. For ARTM3, the dilution at the ground almost vanishes, while the dilution above $z{h}_{\mathrm{m}}^{-\mathrm{1}}=\mathrm{0.8}$ increases to 40 %, and the accumulation in the middle of the PBL increases to 18 %. After the second hour, no further changes are observed, as can be seen in Fig. 3f for the ARTM2 turbulence model at $z{h}_{\mathrm{m}}^{-\mathrm{1}}\approx \mathrm{0.3}$. Time series for other heights and other turbulence models show similar behaviour and are given in Sect. S2 in the Supplement.

This well-mixed condition test shows that the simulation result systematically overestimates the concentration values at $z{h}_{\mathrm{m}}^{-\mathrm{1}}\approx \mathrm{0.3}$ for the ARTM2, PRFMOD, and MODHANNA models after the second hour. Near the surface, which is important for estimation of exposure to the population, the concentration values are underestimated. In both cases, the errors are only 5 %–6 %. At the top of the PBL, the models underestimate the expected concentration significantly. The ARTM3 turbulence model shows the smallest deviation from the mean domain concentration near the ground, but it overestimates the concentration in the middle of the PBL before substantially underestimating the mixing layer top. Below $z{h}_{\mathrm{m}}^{-\mathrm{1}}=\mathrm{0.8}$ the DEGRAZIA turbulence model performs the best. At the top of the PBL the model decreases well below the expected concentration. All the tested turbulence models can be assumed to be acceptable for simulations in very unstable atmospheric conditions, but the partly large deviations of the concentration from the expected values at certain heights have to be taken into account when interpreting model results. Under low-wind conditions (1 m s⁻¹ at 10 m height), the deviation from the uniform concentration is similar (see Sect. S2 in the Supplement).

5 Comparison of ARTM simulation with airborne observations

The comparison of atmospheric dispersion simulation results with measurements near the ground is not sufficient to derive any conclusions about the three-dimensional structure of simulated emission plumes. To study the agreement of simulated and observed plume dispersion it is inevitable to use observations that resolve the structure of the real plume. Since ARTM simulates the emissions of nuclear facilities with source heights of mainly 100 to 200 m, it is useful to choose observational data originating from similar height levels. In this work, we present a comparison of ARTM simulations with airborne CO₂ observations within the PBL. In such a case, the comparison is challenging because of the turbulent character of the PBL. As pointed out by Brunner et al. (2023), observations only provide snapshots of the real world, while simulations provide one realisation of a multitude of stochastic representations of the real world. Simulations with slightly perturbed initial conditions could result in different dispersion patterns of the plume. Furthermore, simulation results and observations may have different spatial and temporal resolutions and uncertainties, which complicates the comparison of simulations with observations (Farchi et al., 2016). Thus, in this work, the comparison of simulation results with observations for the five turbulence models is given using rather general plume characteristics, such as the plume width per transect and maximum mixing ratios.

5.1 Observational data

The aircraft observations used for this investigation originate from the Carbon Dioxide and Methane Mission (CoMet 1.0) (Fix et al., 2018; Luther et al., 2019; Fiehn et al., 2020; Gałkowski et al., 2021; Wolff et al., 2021; Krautwurst et al., 2021; Andersen et al., 2023; Brunner et al., 2023). The campaign took place in May and June 2018 and involved three aircraft performing in situ and remote sensing measurements. The objective was to study CO₂ and methane (CH₄) emissions from different sources in Europe, including power plants, and to compare the different observational methods.

For the evaluation of ARTM, airborne in situ CO₂ measurements in the vicinity of the Bełchatów lignite power plant in Poland were used (Fiehn et al., 2020; Kostinek et al., 2021). An overview map with the CO₂ mixing ratios along the flight path is shown in Fig. 4. The in situ measurements were performed on 7 June 2018 between 13:00 and 15:00 UTC aboard the DLR Cessna Grand Caravan 208B. One transect on the upwind side of the emitter was performed at the beginning in order to derive the mean background CO₂ mixing ratio $c_{{\mathrm{CO}}_{\mathrm{2}}}$ = 401.2 ppmv. The exhaust plume of the power plant was probed during several transects on the downwind side at heights between 500 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ and 1.7 $\mathrm{km}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ They form two wall patterns at meridional distances of approx. 13 km (wall 1) and 23 km (wall 2) and a single transect at approx. 6 km from the source.

Figure 4 Map showing the flight path of the DLR Cessna aircraft in the vicinity of the Bełchatów lignite power plant (blue star), colour-coded by the in-situ-measured CO₂ values. Transects were performed both east (upwind side) and west (downwind side) of the emitting power plant. The red box indicates the simulation domain.

(KTA 1508, 2017)(TA Luft, 2002)(TA Luft, 2002)(SkyscraperPage, 2023)(Brunner et al., 2023)(Brunner et al., 2023)(Farr et al., 2007)

Table 4 Input parameters needed by ARTM that are constant during the simulation run.

CO₂ was measured with a cavity ring-down spectroscopy analyser (G1301-m, Picarro) specifically modified for the airborne deployment operating at 0.5 Hz. The CO₂ measurement uncertainty is ± 0.15 ppmv, and the temporal resolution was increased to 1 s by interpolation to make the data comparable with other data collected during the campaign. Details of the measurement equipment and related uncertainties are described by Klausner et al. (2020). The sampling repetition and the velocity of the aircraft result in a spatial distance of about 140 m between the 0.5 Hz data points. The Picarro instrument measures CO₂, methane, and water vapour sequentially, and thus the values are representative of the last third of the measurement interval. Observational data for wind direction, wind speed, and flight height are shown in Sect. S3 in the Supplement.

5.2 Model setup

We chose a simulation domain of 33.3 km × 33.3 km × 1.9 km that covers the horizontal extent of the flight trajectory and vertically extends beyond the mixing layer depth by four simulation levels. The horizontal resolution was 150 m. The extent of the simulation domain with the location of the emission source (two stacks over a distance of 300 m) is shown in Fig. 4. Vertically, the grid spacing gradually increases from 3 to 35 m until 100 m height is reached. Above, 50 m level thickness was used. All level thicknesses are given in Table S2 in the Supplement.

ARTM requires several input parameters: SC, z₀, d₀, orography, several source-specific parameters, and wind speed and direction at one location in the simulation domain. Since there were no stationary ground-based wind measurements available, wind direction, wind velocity, SC, and mixing layer height were derived from the airborne measurements. The actual emission rates are unknown. However, Brunner et al. (2023) estimated the overall CO₂ emission rate according to the generated electrical power of the power plant, resulting in 1503.0 kg s⁻¹ during the measurement flight. This corresponds to 123 % of the annual mean emission rate of 38.4 Mt CO₂ reported by the power plant to the European Pollutant Release and Transfer Register (E-PRTR) for the year 2018. The description of the derivation of SC, z₀, d₀, and the plume rise as well as the orography data are given in Sect. S4 in the Supplement. The parameter values and the origin of the orography data are summarised in Table 4.

ARTM requires radionuclide emission rates in Bq s⁻¹; thus we chose CO₂ consisting of the radioactive isotope ¹⁴C for the simulations. Its decay constant λ=5730 ± 40 years leads to a decay of 5.5 × 10⁻⁶ %, which is negligible within the simulation period. Thus, ARTM's internal emission rates in Bq s⁻¹ can be used as an equivalent for a mass rate in kg s⁻¹ and to convert activity concentration into mixing ratio.

The wind speed (4.4 m s⁻¹) and directions driving the simulation were derived from one flight transect (13:28:03 to 13:33:14 UTC) at a distance of ≈ 13 km to the west of the power plant at a height of ≈ 600 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ This transect is located close to the middle of the simulation domain and is therefore assumed to be representative. The histogram of the wind directions of the transect is shown in Fig. 5. Based on this histogram, two different setups of the model were selected:

• i.

  A single wind direction of 120° (mean of the distribution) was selected, assuming that the wind fluctuations are part of the turbulence spectrum and should therefore be represented by the turbulence parametrisation of ARTM.

• ii.

  Two different wind directions were used alternatingly to drive ARTM, a direction of 106° (mean of all directions < 120°) and a direction of 134° (mean of all directions > 120°). This assumes that part of the wind variation is due to mesoscale variability that cannot be resolved by ARTM's turbulence scheme.

Figure 5 Histogram of the wind directions of the transect chosen for the determination of the wind direction and wind velocity. The transect covers a duration from 13:28:03 to 13:33:14 UTC with the mean position 53.31° N, 19.15° E. The mean measurement height is 599 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ The mean value of the wind direction is 120° (red line). The mean value of the wind directions below 120° is 106° and above 120° is 134° (dashed blue lines).

The first setup was applied to all turbulence models, while the second setup was only tested for ARTM2. The hourly sequence of wind inputs for the model is summarised in Sect. S5 in the Supplement.

Figure 6 Modelled CO₂ mixing ratio for the case of one wind direction in (a) ARTM2, (b) ARTM3, (c) PRFMOD, (d) MODHANNA, and (e) DEGRAZIA and two wind directions in (f) ARTM2. The wind directions and speeds are given in Table S3 in the Supplement; the input parameters are given in Table 4. Mixing ratios of the simulated plume (averaged over the simulation time) at heights between 750 and 800 m and the in situ data along the flight path between 700 and 800 m are shown in logarithmic scale in ppmv. The case of two wind directions in (f) shows the mean CO₂ mixing ratio of 2 subsequent hours for the duration of the measurement flight from 13:00 to 15:00 UTC. The background CO₂ mixing ratio of 401 ppmv is subtracted from the observation.

5.3 Horizontal dispersion

The mixing ratio maps simulated with the five turbulence models at a height of 750 to 800 m are shown in Fig. 6 together with the observations between 700 and 800 m. We subtracted the background CO₂ mixing ratio of 401 ppmv from the observations to make them comparable to the simulation results.

The simulated and observed mixing ratios of the plumes are of the same order of magnitude. The simulated plumes show the mean wind direction to be in agreement with the observed one; however, the meandering behaviour of the real plume can be observed at transects 1, 2, and 3 in Fig. 6, revealing that this behaviour is not covered by all turbulence models. The mixing ratio profile in the lateral (crosswind) direction simulated by ARTM resembles a Gaussian distribution. This is expected for a constant wind direction and wind speed (Thykier-Nielsen et al., 1999).

The different turbulence models clearly affect the simulated plume widths. The ARTM2 turbulence model simulates the narrowest plume. The ARTM3 model results in a slightly wider plume, but compared to the observations both are too narrow. The PRFMOD and DEGRAZIA turbulence models show much broader plumes that cover the observed one to a large extent. The widest plume is simulated by the MODHANNA turbulence model and is in good agreement with the observed plume width. The width of the plumes of the turbulence models is mainly attributed to the horizontal wind speed fluctuations and Lagrangian correlation times displayed in Fig. 2. The highest values for σ_y and T_Ly are simulated by MODHANNA, PRFMOD, and DEGRAZIA followed by ARTM3 and ARTM2 in the upper half of the PBL, in agreement with the simulated plume width.

Figure 7 Simulated (red line) and measured (black line) CO₂ data along the flight path together with the flight height (dotted blue line) within the simulation domain. The (a) ARTM2, (b) ARTM3, (c) PRFMOD, (d) MODHANNA, and (e) DEGRAZIA turbulence models and the (f) ARTM2 turbulence model with two alternating wind directions. The transects shown in Fig. 6 are shaded in grey.

Figure 7 shows the simulated and observed plumes of the different turbulence models together with the flight height above ground level along the flight path. Transects 1, 2, and 3 are shaded in grey. Data above the simulated boundary layer top are excluded from the figures. In Fig. 7, the simulated maximum CO₂ mixing ratios of ARTM2 are larger at all transects compared to the observations. Within the simulated boundary layer, this deviation reaches 300 % at 14:08 UTC and is attributed to the too-narrow simulated plume. With increasing plume width of the different turbulence models, the maximum mixing ratios decrease (see Fig. 7a–e). The ARTM3, PRFMOD, MODHANNA, and DEGRAZIA turbulence models simulate mixing ratio peaks similar to or below the observed values. Due to dispersion the mixing ratio maximums decrease with increasing distance from the source for all models, in agreement with the observation. It is important to point out that simulated mixing ratio values are highly dependent on the emission rates.

The simulation gives 1 h averages of the exhaust plume, which is expected as the mean of several realisations of meandering plumes. It is not expected that simulated values are much larger than the observed ones but can occur if the width of the simulated plume or the mixing layer depth are underestimated or the emission rate is overestimated.

An alternative to modelling the meandering behaviour via the turbulence is the usage of alternating wind directions for subsequent simulation hours for the ARTM2 turbulence model to explicitly simulate the meandering plume (Figs. 6f and 7f). Simulation results from subsequent hours are combined by calculating the average concentrations. The wind direction derivation is explained in Sect. 5.2. This method generates the widest plume covering the observations and mimics the structure of two maxima at transect 1. However, these two observed maxima originate from snapshots of the meandering plume and are not expected to be reproduced by the time-averaged simulation. Moreover, physically unrealistic plateaus of mixing ratios are simulated in wall 1 and a single narrow mixing ratio peak in wall 2, which is a result of the alternating wind directions. Mixing ratio maps of simulations and observations at other selected heights are given in Sect. S6 in the Supplement.

Figure 8 Cross section of the simulated plumes at wall 1 of the observations for the different turbulence models: (a) ARTM2, (b) ARTM3, (c) PRFMOD, (d) MODHANNA, and (e) DEGRAZIA. (f) The case of two wind directions for ARTM2. The x-axis y distance is in a south–north orientation. The dashed line at 1650 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ marks the simulated mixing layer top.

5.4 Vertical dispersion

For the analysis of the vertical plume behaviour, the cross sections of the simulated plumes at wall 1 are presented in Fig. 8. The narrowest simulated plume is obtained by the ARTM2 model and underestimates the width of the observed plume at heights from 600 m to 1400 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ The plume of the ARTM3 model is slightly wider throughout the PBL. In both the ARTM3 and the PRFMOD model, the values of σ_y (T_Ly) decrease (increase) with height (see Fig. 2). While these opposing trends cancel each other out for the ARTM3 model, they lead to a slight increase in lateral dispersion with height for the PRFMOD model. The vertical profiles of σ_y and T_Ly of the MODHANNA model shown in Fig. 2c and d appear to lead to a slightly increasing dispersion, too. The DEGRAZIA model, in contrast, shows a constant behaviour for both σ_y and T_Ly. Below 200 m the widths of all five simulated plumes decrease towards the surface.

All turbulence models show a slight decrease in the mixing ratio with increasing height at a constant distance from the source (see Fig. 7), which agrees with observations. From the cross sections at wall 1 (Fig. 8), the average horizontal mixing ratio profiles are derived and shown in Fig 9. Except for the DEGRAZIA model, the decreasing mixing ratio with increasing height above 600 m can be recognised here as well. However, in Fig. 7 the highest maximum mixing ratios at wall 1 (13:25 to 14:10 UTC) occur at transect 2 for the measurements. The simulations, having the highest concentration values, instead show very similar peaks at the two lowest transects in wall 1.

Figure 9 Profile of the average horizontal mixing ratio of the six simulation cases at wall 1 (see Fig. 8). The dashed line at 1650 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ marks the simulated mixing layer top.

In contrast to the Gaussian lateral mixing ratio distribution of the plume in Fig. 8a, the ARTM2 turbulence model with two alternating wind directions (Fig. 8f) shows the uniform mixing ratio distribution in the plume's core region (mixing ratio > 1 ppmv), as shown in Fig. 7f.

The cross sections of wall 2 in Fig. S12 in the Supplement show a similar behaviour of the plumes. The measured data show a large variation in the plume width on the different transects, emphasising the meandering and turbulent character of the real plume. Furthermore, it can be recognised that the real plume is not entirely recorded; the transects are too short at this wall.

5.5 Validation and uncertainty evaluation

In order to quantify the simulations' uncertainty, we investigate the deviations of the simulated and the observed CO₂ mixing ratios in the plume through probability distributions (PDs), comparisons of integrated plume mass, and point-to-point mixing ratio comparisons.

Figure 10 Probability distribution (bars) and cumulative probability distribution (lines) of the simulated and measured mixing ratios of the five turbulence models. The PDs are normalised according to the maximum mixing ratio of the measurements, and the integral of the simulated and measured PDs is equal. Mixing ratio values below 1 ppmv are not considered in the PDs and CPDs.

The deviation of model results and measurements in a plume can be accessed by the comparison of the PDs and the cumulative probability distributions (CPDs) of simulated and observed CO₂ mixing ratios in the plume. The PDs of the simulation and measurement are normalised to the maximum mixing ratio of the measurements with the integrals of the simulated and measured distributions being equal. To get rid of the mixing ratio fluctuation in the excess mixing ratio of the measurement, mixing ratio values below 1 ppmv are not taken into account. The PDs and CPDs of the five different turbulence models and the observations for all transects below the simulated boundary layer top are given in Fig. 10. The PDs of the simulated and measured plume show the occurrence of mixing ratio values relative to the maximum mixing ratio of the measurements. There is an overestimation of the simulated maximum mixing ratios for the ARTM2 and ARTM3 turbulence models. The high number of data points at approx. 20 % of the maximum mixing ratio of the measurements is due to the fine structure, shoulders beside peaks, and broad indistinct peaks of the plume not represented in the simulations. MODHANNA can be identified as the turbulence model that shows the best agreement with the observations concerning the PDs and CPDs; i.e. the occurrence of the mixing ratio values is the most similar to that of the measurement. To quantify the similarity and to decide whether simulations and measurements are significantly different, three statistical test were applied: the Z test, the Kolmogorov–Smirnov (KS) test, and the Cramér–von Mises (CvM) test (Conover, 1980; Wilks, 2006; University of Oregon, 2020). The Z statistic represents the distance between the means of two PDs normalised to the standard error. Statistics below 2 indicate no significant difference between the distributions. Additional interpretation limits are given in Sect. S7 in the Supplement. The KS statistic represents the supremum of the distance between two CPDs, while the CvM statistic is proportional to the integral of the distances between two CPDs. For both we assumed a significance level of 0.05. The statistics and their p values (in brackets) are summarised in Table 5. The three statistical tests show that all simulated mixing ratio distributions differ significantly from the observed one. Nevertheless, the statistics can be used to rank the models. MODHANNA shows the best agreement with the observations; i.e. the distribution of mixing ratio values in the transects is the most similar to that of the observations compared to the other turbulence models. The statistical tests rank ARTM3 second, but this may be biased by the statistical tests being very sensitive to deviations in the regions of the PDs with high numbers of low mixing ratio values. We want to point out that the results do not mean that the MODHANNA model produces mixing ratio peaks that are structured like the observed ones, but the relative occurrence of mixing ratio values is the most similar.

Table 5 Z statistics, Kolmogorov–Smirnov (KS) statistics, and Cramér–von Mises (CvM) statistics of the mixing ratio distributions of the five turbulence models. The p values are given in brackets. The significance level is 0.05.

To compare the simulation results, the integral of the mixing ratio values along the flight path below the simulated boundary layer top (see Fig. 7) within the plume is shown in Table 6. We used the method above to get rid of the baseline fluctuations in the excess mixing ratios to calculate the integrals. This procedure is also applied to the simulations. Except for ARTM3, there is good agreement between the modelled and measured data: the deviation is less than 13 %. Concerning ARTM3, there is a strong vertical gradient in the simulated mixing ratios of the plume above 700 m, as is illustrated in Figs. 8b and 9. Tracers are more strongly diluted (accumulated) in the upper (lower) half of the PBL than in the other turbulence models. This corresponds to the findings presented in Sect. 4. Since the flight path is mainly located in the upper half of the PBL, the integral along the flight path results in a lower value for ARTM3. The higher mixing ratios in the lower half of the PBL might become important when simulations are used for radiation exposure assessment. The results suggest that the original assumption of the emission rate may not deviate much from the actual value. However, observations below 600 m are necessary to get a more complete comparison of the simulated and actual plumes.

Table 6 Integrals of mixing ratio values (values below 1 ppmv are not considered) along the flight path for simulation A_sim and observation A_obs within the simulated PBL (see Fig. 7) given in ppmv km and their ratio.

Figure 11 Density scatterplots of the simulated mixing ratios of the five turbulence models, ARTM2 (a), ARTM3 (b), PRFMOD (c), MODHANNA (d), and DEGRAZIA (e), against the observations. Single data points in a bin are indicated with +; more data points in a bin are colour-coded. The regression with slope m=1 (dotted black line) represents the identification of the simulation with the measurement, the dotted–dashed red line represents a deviation from the regression with m=1 by a factor of 2, and the solid black line represents the orthogonal linear fit to the data points. Fac2 gives the percentage of data points with deviations of not more than a factor of 2 from the regression m=1. The residual variance of the orthogonal fit is given by ${\mathit{\sigma }}_{\text{res}}^{\mathrm{2}}$.

The deviation between the simulations of the five turbulence models and the observation at a specific position can be assessed using density scatterplots as given in Fig. 11. All mixing ratio values larger than 1 ppm along the flight path below the simulated boundary layer top are considered. The regression with slope m=1 is shown with a dashed line and represents the equality of the simulated and observed mixing ratios. The deviation from this equality by a factor of 2 is confined by the dashed–dotted red lines. It is not expected to find a lot of data points at the regression m=1 due to the fundamental differences in the data set properties of the simulation and observation. However, a large number of data points within a deviation of a factor of 2 decreases the uncertainty. The percentage of data points within these borders is represented as Fac2 (Fig. 11). Low Fac2 values can also be explained by the large number of measurement data points outside the simulated plume because they are too narrow. The smallest Fac2 is derived for the ARTM2 model. This coincides with the unbalanced distribution of the data points around the regression m=1, with the simulation overestimating the observed mixing ratios and simultaneously simulating too-narrow plumes. This is represented by the orthogonal regression of data points (black line) given in the figure with a slope above 3. The residual variance ${\mathit{\sigma }}_{\text{res}}^{\mathrm{2}}$ quantifies the scattering of data points. The large value for ARTM2 indicates a less compact data point distribution. ARTM3 shows more balanced and compact data but still distinctly overestimates mixing ratios and underestimates plume widths. PRFMOD, MODHANNA, and DEGRAZIA show similar properties with Fac2 > 50 %; more compact, well-balanced data; and less overestimated mixing ratios and underestimated plume widths, with the MODHANNA model performing slightly better for the given measurement and turbulence conditions.

The uncertainty in the CO₂ measurement device of ± 0.15 ppmv is at least 1 order of magnitude smaller than the measured enhanced CO₂ concentrations. Thus, the measurement uncertainty has only a minor impact on the evaluation results.

6 Conclusions

In this work we present an extensive evaluation of ARTM with three different elements: a sensitivity analysis, an analysis of turbulence models, and a comparison with aircraft observations. Based on the sensitivity analysis, we identified the stability class to be the most important input parameter, followed by the roughness length, the source height, and the displacement height factor. Therefore, special care has to be taken to determine the stability class for a simulation because uncertainties in this parameter cause large uncertainties in model results. This emphasises the general disadvantage of the rather coarse stability class concept being used to describe atmospheric turbulence. A finer classification or a continuous parameter such as the Obukhov length could be a better option but would generally require detailed measurements of turbulence parameters such as friction velocity and sensible heat flux.

In addition to the three turbulence models already implemented in ARTM 3.0.0, two other turbulence models, the MODHANNA model and the DEGRAZIA model, were implemented and tested. Evaluation of the models by applying the well-mixed condition test showed that the ARTM2, PRFMOD, and MODHANNA models produced a moderate deviation of an initially uniform concentration profile under unstable atmospheric boundary layer conditions. Underestimations of the uniform concentration occur primarily at the ground and at the top of the boundary layer with up to 10 %, while overestimations occur in between with up to 7 % at one-third of the PBL. Under the same conditions the ARTM3 model produces the strongest deviations: up to 4 times higher than the other models. However, near the ground ARTM3 performs the best. The DEGRAZIA model showed a less inhomogeneous profile with deviations from the uniform concentration of 5 % or less below z=0.8h_m. The discrepancies under 6 % below 80 % of the boundary layer height show good mixing properties of the unstable planetary boundary layer for the ARTM2, PRFMOD, MODHANNA, and DEGRAZIA models.

Three-dimensional airborne in situ observational data measuring a power plant emission plume were compared to ARTM simulation results. The time resolution of ARTM results is 1 h, which is larger than the expected timescale of the observed meandering plume, and, therefore, ARTM is expected to capture the time-integrated real plume and not the fine structures on small scales. ARTM simulated the mean wind in agreement with the observations throughout the simulation domain. The different turbulence models simulate plume mixing ratios of the same order of magnitude as the measurements, although the exact mixing ratio values depend on the emission rate. ARTM2 underestimates the plume spread under very unstable conditions and overestimates the maximum mixing ratio by a factor of 2 or more. The ARTM3 model produces only slightly wider plumes in the lateral direction but lower maximum mixing ratio values at the upper half of the PBL. This is attributed to the inhomogeneous vertical mixing and the horizontal turbulence parametrisation of the ARTM3 turbulence model under unstable conditions. The other turbulence models, PRFMOD, MODHANNA, and DEGRAZIA, simulate a wider plume spread in the range of the measurements. Maximum mixing ratios are close to the measurements, and the integrals of the mixing ratios of the simulations and observations along the flight path are comparable. The models were evaluated with measurements at heights larger than 600 m and hence do not cover the heights below 600 m. The differences in the temporal and local resolutions of the simulations and measurements lead to differences in the distributions of mixing ratio values. According to Fig. 10, under unstable conditions, all turbulence models underestimate the occurrence of mixing ratio values by around 20 % of the maximum mixing ratio of the measurements, which is a result of the fine structure of the plumes. The smallest deviations in PDs and CPDs are found for MODHANNA. Using point-to-point comparisons, the ARTM2 model showed the largest deviations from the measured plume; ARTM3 showed better agreement; and PRFMOD, MODHANNA, and DEGRAZIA showed comparable good performances, with MODHANNA matching the measurements slightly better than the others.

With the results of this study we showed that ARTM is able to simulate the extension and mixing ratios of a plume under unstable stratification conditions when the proper turbulence model is used. The ARTM3 model was shown to be a suitable turbulence model for radiation exposure assessment when conservative long-term simulations are requested. However, the PRFMOD, MODHANNA, and DEGRAZIA models simulate the exhaust plume closer to real exhaust plumes under the given conditions and under the limitations of temporal and spatial uncertainties. Within this validation using the in situ data from the Bełchatów power plant, the MODHANNA turbulence model performed slightly better. However, this ranking cannot be generalised to other stability conditions without further comparison studies. Further analyses with known emission terms with different atmospheric turbulence properties, stable and neutral conditions, are necessary and can lead to better validation of ARTM. We encourage the comparison of other turbulence formulations, such as the one suggested by Hanna (1982). This would also reveal differences between the model of Hanna (1982) and the DEGRAZIA model and extends the analyses given by Carvalho et al. (2002) to larger heights. The collection of measurement data in the upper and lower half of the PBL and transects sampling the entire extent of a plume are beneficial. Moreover, the use of the Obukhov length as a measure of atmospheric stability is encouraged.