FLEXPART version 11: improved accuracy, efficiency, and flexibility

Numerical methods and simulation codes are essential for the advancement of our understanding of complex atmospheric processes. As technology and computer hardware continue to evolve, the development of sophisticated code is vital for accurate and efficient simulations. In this paper, we present the recent advancements made in the FLEXible PARTicle dispersion model (FLEXPART), a Lagrangian particle dispersion model, which has been used in a wide range of atmospheric transport studies over the past 3 decades, extending from tracing radionuclides from the Fukushima nuclear disaster, to inverse modelling of greenhouse gases, and to the study of atmospheric moisture cycles.

This version of FLEXPART includes notable improvements in accuracy and computational efficiency. (1) By leveraging the native vertical coordinates of European Centre for Medium Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) instead of interpolating to terrain-following coordinates, we achieved an improvement in trajectory accuracy, leading to a ∼8 %–10 % reduction in conservation errors for quasi-conservative quantities like potential vorticity. (2) The shape of aerosol particles is now accounted for in the gravitational settling and dry-deposition calculation, increasing the simulation accuracy for non-spherical aerosol particles such as microplastic fibres. (3) Wet deposition has been improved by the introduction of a new below-cloud scheme, by a new cloud identification scheme, and by improving the interpolation of precipitation. (4) Functionality from a separate version of FLEXPART, the FLEXPART CTM (chemical transport model), is implemented, which includes linear chemical reactions. Additionally, the incorporation of Open Multi-Processing parallelisation makes the model better suited for handling large input data. Furthermore, we introduced novel methods for the input and output of particle properties and distributions. Users now have the option to run FLEXPART with more flexible particle input data, providing greater adaptability for specific research scenarios (e.g. effective backward simulations corresponding to satellite retrievals). Finally, a new user manual (https://flexpart.img.univie.ac.at/docs/, last access: 11 September 2024) and restructuring of the source code into modules will serve as a basis for further development.

1 Introduction

Atmospheric transport modelling plays an important role in understanding many of the complex interactions within our atmosphere. Traditionally, Eulerian methods have been widely used for such modelling. These methods solve the atmospheric advection–diffusion equation on a fixed grid and give the concentration of tracer gases or aerosols at each grid point as a function of time. Eulerian methods offer a convenient means of accounting for all processes that occur during transport, including nonlinear atmospheric chemistry. However, they typically have a relatively high level of numerical diffusion (e.g. Reithmeier and Sausen, 2002; Cassiani et al., 2016), limiting their capability of preserving plumes over long transport distances (Rastigejev et al., 2010), and the spatial resolution of the concentration fields is limited by the grid spacing of the model. Lagrangian methods, on the other hand, follow individual “computational particles” (from now on referred to as “particle”) and simulate their transport based on large-scale winds from Eulerian input data and parameterised small-scale fluctuations. They are independent of the computational grid and can therefore produce concentration fields with potentially infinitesimal spatial resolution at relatively low computational cost, especially if several tracers are to be transported simultaneously. Another advantage of Lagrangian methods is that they can provide a direct link between regions along particle trajectories.

Atmospheric transport models using the Lagrangian method with the parameterisation of sub-grid processes are typically referred to as Lagrangian particle dispersion models (LPDMs). One prominent LPDM is the FLEXible PARTicle dispersion model (FLEXPART), which was developed nearly 3 decades ago to simulate the dispersion of radionuclides resulting from nuclear disasters such as the Chernobyl accident (Stohl et al., 1998). Since its inception, FLEXPART has proven to be a valuable tool for studying a wide range of environmental problems in both research and operational settings. Its applications now cover large parts of atmospheric research, including the simulation of the transport of heat and water in the atmosphere (e.g. Baier et al., 2022; Peng et al., 2022), volcanic and wildfire plumes (e.g. Stohl et al., 2006, 2011; Moxnes et al., 2014), transport and fall-out after nuclear accidents or explosions (e.g. Stohl et al., 2012; Arnold et al., 2015), transport of aerosols such as dust (e.g. Groot Zwaaftink et al., 2017; Ryder et al., 2019) and microplastics (e.g. Evangelou et al., 2024; Evangeliou et al., 2020), the interpretation of biogenic secondary organic aerosol compound measurements (e.g. Martinsson et al., 2017), the transport of pollutants into remote regions like the Arctic (e.g. Dada et al., 2022; Zhu et al., 2020), the interpretation of ice cores (e.g. Eckhardt et al., 2023; McConnell et al., 2024), and the modelling of emission sensitivities for greenhouse gas atmospheric concentrations (e.g. Bergamaschi et al., 2022; Vojta et al., 2022).

Besides FLEXPART, several other LPDMs exist, such as HYSPLIT (Draxler and Hess, 1998), STILT (Lin et al., 2003), TRACMASS (Döös et al., 2017), MPTRAC (Hoffmann et al., 2022), ATTILA (Brinkop and Jöckel, 2019), NAME (Jones et al., 2007), SPRAY (Tinarelli et al., 2000), and the Canadian Meteorological Centre (CMC)'s dispersion modelling suite (D'Amours et al., 2015). FLEXPART combines the following capabilities: (i) a detailed parameterisations for Gaussian (Stohl and Thomson, 1999) and non-Gaussian (Cassiani et al., 2015) turbulence in the atmospheric boundary layer (ABL), which both take into account the vertical gradient of air density; (ii) a sub-grid convection parameterisation (Forster et al., 2007); (iii) treatment of radioactive decay in the atmosphere and on the ground (Stohl et al., 1998); (iv) detailed parameterisations for dry and wet deposition and gravitational settling (Stohl et al., 2005); (v) linear chemical reaction with hydroxyl radicals (Pisso et al., 2019); (vi) the capability of running the model backward in time to create sensitivities of atmospheric concentrations and depositions to emission sources (Seibert and Frank, 2004; Eckhardt et al., 2017) that can be interfaced directly with an inverse modelling code for determining emissions (Thompson and Stohl, 2014); and (vii) a domain-filling option, where the entire atmosphere can be filled with particles (Stohl and James, 2004), which is useful for producing Lagrangian climate diagnostics (e.g. Schicker et al., 2010; Baier et al., 2022) and three-dimensional concentration fields of greenhouse gases that may serve as initial conditions for inverse modelling (Groot Zwaaftink et al., 2018; Vojta et al., 2022).

This paper describes a new model version, FLEXPART 11, that adds four important new features to FLEXPART's capabilities, as well as other improvements. First, when driven with data from the European Centre for Medium Range Weather Forecasts (ECMWF), FLEXPART 11 offers the option to keep the original model layers intact for the transport with resolved-scale winds, instead of interpolating to a coordinate system with fixed heights above terrain as previous versions did (see Sect. 3). We show that this improves the overall accuracy of FLEXPART simulations. Second, FLEXPART 11 can calculate the drag coefficient of non-spherical aerosol particles of various shapes, which improves the calculated gravitational settling of such particles substantially and thereby has a large effect on their overall simulated transport in the atmosphere (see Sect. 4). Third, a new below-cloud scavenging scheme for aerosols is introduced, as well as an improved interpolation scheme for precipitation and a new cloud layer identification scheme (see Sect. 5). Last, linear chemical reactions are implemented based on the initial work of Henne et al. (2018), who developed the FLEXPART CTM from FLEXPART 8 (where CTM stands for chemical transport model), which was first described in Groot Zwaaftink et al. (2018) (see Sect. 6). In addition to enhancing the model's core functionality, we have improved the user interface by incorporating additional options for ways of running FLEXPART (see Sect. 9).

Computational efficiency is another important consideration. Legacy models such as earlier versions of FLEXPART or HYSPLIT (Draxler and Hess, 1998) were originally designed for serial processing. An important step was made with FLEXPART 10.4, which introduced the option to use the Message Passing Interface (MPI) parallelisation (Pisso et al., 2019). However, especially with the increasing size of the meteorological input data, the memory requirement became problematic in parallel mode. For this reason, we opted for a different parallelisation strategy, using OpenMP (Open Multi-Processing), following the FLEXPART CTM (Henne et al., 2018) (see Sect. 8). Furthermore, numerous features added over a period of more than 25 years that often deviate from the original coding standards make the FLEXPART 10.4 source code difficult to understand and maintain. Models that were created recently, such as MPTRAC (Hoffmann et al., 2022), are designed from scratch to run on current hardware, for instance, by utilising the acceleration offered by graphics processing units (GPUs), while not being under the constraint of a legacy code base. However, as FLEXPART is a model that is well tested for a wide range of applications, offers many features not available in other models, and has a large user community, a complete rewrite of the code did not seem appealing. Instead, the approach that started with FLEXPART 8 and 10 to introduce features of modern FORTRAN standards into the source code was continued. In FLEXPART 11, all programme subunits have been encapsulated into modules, and more use has been made of whole-array syntax.

We validate our changes by comparing FLEXPART 11 to tracer experiments (Sect. 7). We also test the conservation of meteorological tracers (Sect. 3.2), the degree to which particles in a global-domain-filling simulation stay well mixed (Sect. 3.3), and the removal processes by reproducing the Fukushima Daiichi Nuclear Power Plant accident (Sect. 5.4).

For users unfamiliar with FLEXPART, a short overview of the FLEXPART 11 capabilities, input data, and code reorganisation can be found in the Supplement (see Appendix A). Accompanying this paper is a completely revised technical documentation, including a download and installation guide, of FLEXPART (https://flexpart.img.univie.ac.at/docs, last access: 11 September 2024); a snapshot of the code used in this work is available at https://doi.org/10.5281/zenodo.12706632 (Bakels et al., 2024). We recommend that users of FLEXPART consult this living document for updates and future code developments.

2 Input data

FLEXPART calculates particle trajectories using meteorological data interpolated in time and space from time sequences of three-dimensional fields of, e.g., wind velocities, density, temperature, specific humidity, and cloud liquid and ice water content. FLEXPART 11 supports two input formats of data from the European Centre for Medium Range Weather Forecasts (ECMWF) and from the National Centers for Environmental Prediction (NCEP) forecast systems, namely the Integrated Forecasting System (IFS) and Global Forecast System (GFS), respectively. See Appendix A1 for further details.

For the examples provided in this paper, we use the most recent re-analysis dataset of ECMWF, ERA5 (Hersbach et al., 2020), with hourly 0.5° × 0.5° data on all of the 137 model levels as input to FLEXPART.

3 Grid-scale advection based on ECMWF fields on native coordinate surfaces

While we know that FLEXPART does not quickly produce strong anomalies in particle distributions in domain-filling simulations (Stohl and James, 2004), the degree to which older versions of FLEXPART obeyed the well-mixed criterion with respect to particle positions over periods of months to years was not entirely satisfying. For this reason, we decided to switch from z to η coordinates (see Sect. 3.1), as it might reduce interpolation errors (see Sect. 3.2) and improve particle transport accuracy and particle distribution in domain-filling simulations (see Sect. 3.3).

3.1 The η-coordinate system

The ECMWF's IFS employs hybrid pressure base vertical coordinates, where the η surfaces are given by

$$\begin{array}{}\text{(1)}& {\mathit{\eta }}_k=a_k/P_{\mathrm{0}}+b_k.\end{array}$$

Here, η_k is the value of η at the kth model level; P₀= 101 325 Pa is the reference pressure; and a_k and b_k are coefficients chosen such that the levels closest to the surface follow the topography, the highest levels are pressure surfaces, and intermediate levels are hybrids between the two. The pressure in the η-coordinate system is determined by the surface pressure P_s, which varies in space and time:

$$\begin{array}{}\text{(2)}& P_k(\mathit{\lambda },\mathit{\varphi },t)=a_k+b_k{P}_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },t).\end{array}$$

Here, λ, ϕ, and t denote longitude, latitude, and time, respectively. In this system, the vertical velocity is also represented in η coordinates, where negative values indicate upward motion. Detailed information on the η-coordinate system can be found in the ECMWF – IFS documentation (2023). While ECMWF meteorological data can also be downloaded on pressure levels, pressure level data are less accurate, since they are interpolated between the η levels, and fewer levels are available.

The boundary layer turbulence scheme utilised by FLEXPART (Hanna, 1982; Stohl et al., 1998) is based on the terrain-following geometric height z as a vertical coordinate, with turbulent velocities expressed in units of metres per second (m s⁻¹). For this reason, and to avoid frequent time-consuming coordinate transformations, older versions of FLEXPART used an internal terrain-following Cartesian (z) coordinate system, and all meteorological data were interpolated to these internal z coordinates. However, this approach introduced additional interpolation errors, since the meteorological data needed to be interpolated first from η to z levels and then to particle positions.

To mitigate these interpolation errors and improve computational accuracy, FLEXPART 11 retains the η-coordinate system whenever possible. Only in the ABL, where the Hanna turbulence scheme (Hanna, 1982) is used, are the z levels kept as the main coordinate system. However, interpolation errors are smaller in the ABL than in higher layers of the atmosphere, since the lowest η levels follow the topography, and z levels are chosen to coincide closely to η levels, with the best agreement near the surface and an almost perfect agreement for surface pressure values that are typical of the sea level. Everywhere else, the native η-coordinate system is used in FLEXPART 11 to interpolate meteorological data to the particle positions and advect the particles. The η coordinates are utilised in FLEXPART 11 by default but can be switched off at the time of the compilation by building the FLEXPART 11 executable with the addition of eta=no. In this case, meteorological data are interpolated to z coordinates as in previous FLEXPART versions. Below, we demonstrate that the switch of coordinate system clearly increases the accuracy of the model computations. We evaluate conservation errors in the quasi-conservative properties in Sect. 3.2 and the error in the particle density distribution of a global domain-filling simulation in Sect. 3.3.

3.2 Quasi-conservative properties

It is not a trivial task to prove that changes in a Lagrangian model lead to improved trajectory accuracy, since ground truth trajectories are usually not available. To show the increase in accuracy obtained as a consequence of keeping ECMWF's η-coordinate system largely intact (except for within the ABL), we followed Stohl and Seibert (1998) and evaluate the conservation of quasi-conservative properties along trajectories (Fig. 1), such as potential vorticity (PV), potential temperature (Θ), and specific humidity (q). These properties are expected to be reasonably well conserved in the absence of diabatic processes. Therefore, for a date in each season, we released 6 million particles that are globally distributed using the domain-filling option and followed their trajectories for 48 h. At heights above 5 km, diabatic processes related to surface interactions can be excluded, and furthermore, the differences between η and z levels are substantial. Therefore, we traced the trajectories of all particles between 5 and 10 km and between 10 and 20 km above ground level using FLEXPART in two simulations using η and z coordinates. We removed particles that left the specified height range, and to avoid diabatic heating by water phase changes, we removed particles that went through regions where the relative humidity was above 90 %. We only used particles outside the subtropics and tropics, excluding the zone between 40° N and 40° S, as we expect the tracer conservation in this region to be worse in general, where the geostrophic balance is weak, and deep convection is frequent. To avoid sampling biases, we selected particles from regions that keep similar density distributions over time (see Fig. 2), resulting in not including particles that reside north of 80° N or south of 80° S. We also switched off both the turbulence and convection parameterisations in FLEXPART and used short time steps (LSYNCTIME=600).

Figure 1 Mean absolute tracer conservation errors (ATCEs) in quasi-conservative particle properties. (a) Potential vorticity; (b) potential temperature; (c) specific humidity. Results are shown for altitudes between 5 and 10 km (indigo lines) and between 10 and 20 km (grey lines) above ground level as a function of time since the initiation of a trajectory. The figure illustrates these errors for both the η (solid lines) and the z vertical coordinate system (dashed lines). Each line represents the mean over around 330 000 trajectories that are selected between 40 and 80° N and 40 and 80° S and a relative humidity below 90 %. Thin lines represent one simulation in each season of 2020, starting on the first day of January, April, July, and October at 00:00 UTC. Each thick line shows the average of the four seasons. All FLEXPART parameterisations (e.g. turbulence and convection) are switched off. Transition periods between ERA5 12-hourly data assimilation windows are highlighted with lavender shades, beginning 9 h after the FLEXPART simulation's start. At the top left of each panel, the error reductions that are gained using the η coordinates as compared to the z coordinates are given in percentages, both excluding and including (in parentheses) data assimilation window transition periods.

Using these trajectories, we computed the absolute transport conservation errors (ATCEs) as a function of time along the trajectory, given by

$$\begin{array}{}\text{(3)}& \mathrm{ATCE}(\mathit{\gamma },t)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|{\mathit{\gamma }}_i\right(t)-{\mathit{\gamma }}_i(t_{\mathrm{0}}\left)\right|,\end{array}$$

where γ_i(t) and γ_i(t₀) are the quasi-conservative properties of particle i at times t and t₀ of its trajectory, and N is the total number of particles in the domain between [40° N, 80° N] and [80° S, 40° S] latitude that fulfilled our selection criteria. ATCE values indicate how strongly the considered variables, on average, deviate from the perfect conservation along the trajectories.

Figure 1 shows ATCE values for PV (left panel), Θ (middle panel), and q (right panel) as a function of time along the trajectories. The first thing to notice is that the ATCE values for PV and Θ show stepwise jumps every 12 h. These jumps occur during the 1 h transition of hourly re-analysis data from one 12 h long data assimilation window to the next in the ERA5 reanalysis production and can be explained by the dynamical inconsistencies between two different assimilation periods. The resulting step-wise increases in conservation errors are a problem inherent to the ERA5 data and not related to FLEXPART. While the inconsistencies themselves do not necessarily lead directly to trajectory position errors, they do indicate uncertainties in the analysed state of the atmosphere, which will lead to trajectory position errors. Notably, Fig. 1 also shows wiggly lines with an hourly periodicity, best visible in the left panel, indicating that the errors grow most rapidly in the middle of two consecutive hourly wind fields. These wiggles can be attributed to the fact that the temporal interpolation errors are largest when the time difference to the closest available wind field is largest (at 30 min in between two wind fields).

Most interestingly for our purpose, Fig. 1 shows that ATCE values for PV, Θ, and q for trajectories between 10 and 20 km are, respectively, 15.4 %, 60.6 %, and 16.9 % lower for the calculations using η coordinates than for those using z coordinates after 48 h, ignoring assimilation window transition periods. Corresponding differences for particles between 5 and 10 km are 8.1 %, 1.8 %, and 1.5 %, respectively. The relative differences are larger for the higher trajectory starting points, which is expected because our screening to avoid regions of turbulence, convection, and clouds is not perfect, and it is to be expected that the trajectories starting from the lower levels are more strongly affected by diabatic processes there, which can lead to large conservation errors, regardless of which coordination system is used, thereby reducing relative differences between the two trajectory calculations. In summary, our results indicate a substantial reduction in tracer conservation errors when avoiding the interpolation of meteorological data from the original η levels to a secondary z-coordinate system and thus an improved accuracy of trajectory calculations with FLEXPART 11 compared to its predecessor versions. While it is impossible to directly quantify improvements in trajectory position accuracy because of a lack of ground truth data, it is expected that better tracer conservation also indicates more accurate trajectories. Similarly, while we cannot use the dynamical tracers for quantifying conservation errors in the lower troposphere, we certainly also expect improvements in trajectory accuracy there.

3.3 Density distribution

Validation of FLEXPART based on the conservation of quasi-conservative quantities is restricted to certain regions of the atmosphere and does not directly allow us to evaluate the accuracy of trajectory positions. Another validation possibility is to test how accurately the well-mixed criterion is fulfilled. The well-mixed criterion states that if a species of passively marked particles without sources or sinks is initially mixed uniformly in the position and velocity space in a turbulent flow, it will stay that way (Thomson, 1987). Generally, Lagrangian transport models could be expected to encounter challenges in this regard, as evident from works in the broader literature discussing specific conservational issues for (stochastic) particle methods, such as Wang et al. (2015), Meyer and Jenny (2004), and Bahlali et al. (2020), making the test especially relevant. Lagrangian particle models, including FLEXPART (Stohl and Thomson, 1999; Cassiani et al., 2015), are usually tested for this criterion for simulations of turbulent dispersion in the ABL on short timescales. However, this criterion should also be fulfilled (and tested) in the atmosphere as a whole and on long timescales, similar to what was done by Cassiani et al. (2016).

To see whether the switch from z to η coordinates improved the particle distribution, we performed two 6-month global domain-filling simulations of a passive air tracer based on z and η coordinates, respectively. In total, 6 million particles were initially globally distributed, perfectly in accordance with the local air density, and ideally, the particle density distribution should closely follow the air density distribution throughout the entire simulation, even though the particles were allowed us to move through the atmosphere without any further constraint on their distribution. We used short time steps, with the basic integration time step LSYNCTIME = 600 s, and both the CTL and IFINE options set to 10 (thus allowing for reduced time steps in turbulence calculations; see Appendix A2.2) to optimise the accuracy of the simulation of turbulent transport in the ABL. For an evaluation of the results, we averaged the densities of the particles over the sixth month after the start of the simulations and compared those to the air densities given by the ECMWF ERA5 data. This is a comprehensive test, as it involves transport in the whole atmosphere and also includes sub-grid-scale parameterisations of convection and turbulence.

Figure 2 Averaged density distribution of particles as compared to the density of air given by the ERA5 data (black lines) after 6 months of running FLEXPART. The results are averaged over the sixth month (July) of the η- (solid blue lines) and z-coordinate (dashed red lines) system simulations, after initially starting with a perfectly fitting particle density distribution. The results are separated into three regions: polar regions (a), midlatitudes (b), and the tropics (c). The right-hand side of each panel shows the absolute errors in the η (blue) and z (red) simulations. Thick horizontal lavender lines represent the average ABL height (solid) and tropopause height (broken).

Figure 2 shows vertical particle density profiles averaged over polar regions, midlatitudes, and tropics, respectively. Overall, the η-coordinate system version reproduces the vertical air density profile more accurately than the z-coordinate system version. In the midlatitudes and the tropics, both model runs reproduce the vertical air density profiles relatively well; however, for the midlatitudes, the η version is clearly performing better than the z version throughout almost the entire depth of the troposphere. For the tropics, better results are seen for the lowest 2 km of the atmosphere. Especially in the ABL, the z version underestimates the air density on average by ∼2 % (midlatitudes) and ∼3 % (tropics). Only in the stratosphere, above about 15 km, does the z version perform slightly better than the η version. Both model runs overestimate the density below about 1 km in polar regions quite substantially, by up to ∼3.6 %. The z version shows somewhat better agreement there, with a weighted average of ∼1.2 % better up to ∼1 km. However, this reverses above about 1 km, and in the polar lower stratosphere, the z version shows a substantial deviation from air density, peaking at ∼20.3 %, and here the η version shows much better agreement. The underestimation of particle density in the tropical troposphere and the overestimation in the polar stratosphere by the z version show that this version tends to move too many particles from the low to the high latitudes and also to higher altitudes near the poles. This deficiency is less pronounced for the η version, which shows a much better agreement overall of particle density profiles with air density profiles in most of the atmosphere (with the exception of the lowest kilometre in the polar regions and in the stratosphere at lower latitudes). This again indicates the improved trajectory accuracy of FLEXPART 11 compared to previous FLEXPART versions. However, it is also clear that even FLEXPART 11 deviates somewhat from how well it is mixed.

We attribute these deviations to a number of possible causes. First, the ECMWF input data are not perfectly mass-consistent, with additional implicit mass flow violations likely being introduced by the interpolation procedures to generate the sub-grid-scale velocity field, following the considerations in Wang et al. (2015) and Meyer and Jenny (2004). Second, aspects related to grid-scale stochastics could also play a role, as particles may be transferred from regions with higher variability in grid-scale winds towards regions with lower variability because there is no drift correction applied to the grid-scale winds. Last, dynamical inconsistencies between data assimilation periods may lead to the production of spurious trajectories.

4 Gravitational settling

The gravitational settling of aerosol particles is called at each time step and added to the vertical motion. In FLEXPART 11, a settling module, settling_mod.f90, has been introduced and contains all the relevant procedures. The module contains several improvements for the calculation of the drag coefficient, C_d, of aerosol particles, which is needed to determine their gravitational settling velocity. While previous versions of FLEXPART only considered spherical particles, FLEXPART 11 is also able to calculate the settling velocity for non-spherical particles and takes into account different types of particle orientations (see below). One limitation is that the gravitational settling is only calculated if a particle carries a single-aerosol species with non-zero mass. This is necessary, since aerosols of different size, shape, or density exhibit different settling velocities and follow different trajectories. In contrast, a particle can represent any number of gaseous tracers, as they all follow the same trajectory.

Figure 3 Comparison of the atmospheric transport of particles of different shapes. Shown are the FLEXPART simulation results of the monthly averaged total atmospheric deposition for daily releases from a point in Vienna, Austria (white dot), in January 2018 of spheres with a diameter of 50 µm (a) and volume-equivalent straight fibres with an aspect ratio (AR) of 50 (b). In both cases, dry-deposition accounts for more than 90 % of the total deposition. Values for the monthly mean horizontal transport distance (D) and its standard deviations are also reported near the bottom in the left and middle panel. The probability density function (PDF) of horizontal travel distance of spheres (light blue) and fibres (dark blue) is displayed in panel (c).

4.1 Spherical particles

Earlier versions of FLEXPART only considered the settling of spherical particles. The drag coefficient of spheres in the laminar flow regime can be calculated with Stokes' law. However, Stokes' law is only valid for the case of

$$\begin{array}{}\text{(4)}& \mathit{Re}={\displaystyle \frac{\mathit{\rho }{d}_{\text{eq}}{w}_{\mathrm{s}}}{\mathit{\mu }}}\ll \mathrm{1},\end{array}$$

where Re is the Reynolds number, ρ the density of air (in kg m⁻³), d_eq the volume-equivalent diameter of the particle (in m), w_s the settling velocity (in m s⁻¹), and μ the dynamic viscosity of the air (in Pa s). Substantial deviations from Stokes' law can be expected for particles with diameters larger than 10–30 µm in a laminar regime (Drakaki et al., 2022; Saxby et al., 2018; Alfano et al., 2011). To account also for larger Reynolds numbers, earlier FLEXPART versions used the drag coefficient approximation by Bird et al. (1960). However, this approximation still shows a significant settling velocity mismatch for Reynolds numbers > 400 compared to other drag coefficient models (Näslund and Thaning, 1991). Therefore, in FLEXPART 11, we replaced this scheme with the drag coefficient formulation of Clift and Gauvin (1971), which is valid for Reynolds numbers up to 10⁵ and is within 6 % of experimentally determined values (Clift et al., 2005):

$$\begin{array}{}\text{(5)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{24}}{\mathit{Re}}}\left(\mathrm{1}+\mathrm{0.15}{\mathit{Re}}^{\mathrm{0.687}}\right)+{\displaystyle \frac{\mathrm{0.42}}{\mathrm{1}+\frac{\mathrm{42}\mathrm{500}}{{\mathit{Re}}^{\mathrm{1.16}}}}}.\end{array}$$

Using this more accurate scheme changes the settling velocities for large particles substantially. For example, particles with a diameter of 100 µm have a settling velocity lower by about 16 % as compared to the old scheme.

4.2 Non-spherical particles

Experiments show that non-spherical particles experience a larger drag in the atmosphere and therefore have lower settling velocities than volume-equivalent spherical ones (Tatsii et al., 2024). To take this into account, we have extended the gravitational settling scheme and the calculation of dry-deposition velocities to allow for non-spherical particle shapes. This was done by implementing the scheme of Bagheri and Bonadonna (2016) with modifications by Tatsii et al. (2024). The scheme determines a particle's drag coefficient C_d based on its shape, providing the following three options for the orientation of a particle in the atmosphere: (i) random orientation; (ii) the particle's maximum projection area being perpendicular to the vector of gravity, which corresponds to its maximum drag or horizontal orientation and thus lowest settling velocity; or (iii) the particle's orientation corresponding to the average of the first two options. The equations defining Stokes' and Newton's drag corrections (k_S and k_N) for the three options are given in Table B2. The drag coefficient is then computed as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{C}_{\mathrm{d}}& =\mathrm{24}{\displaystyle \frac{k_{\mathrm{S}}}{\mathit{Re}}}\left(\mathrm{1}+\mathrm{0.125}{\left({\displaystyle \frac{\mathit{Re}{k}_{\mathrm{N}}}{k_{\mathrm{S}}}}\right)}^{\mathrm{2}/\mathrm{3}}\right)\\ & +\mathrm{0.46}{\displaystyle \frac{k_{\mathrm{N}}}{\mathrm{1}+\mathrm{5330}\frac{k_{\mathrm{S}}}{\mathit{Re}{k}_{\mathrm{N}}}}}.\end{array}\end{array}$$

The scheme was tested in laboratory experiments with microplastic fibres, where it was found to give best results for the average orientation option (Tatsii et al., 2024). The settling velocities of fibres can be less than one-third of those of spherical particles with the same volume, which has major implications for their atmospheric transport and deposition. Figure 3 compares the atmospheric deposition of spheres and equivalent-volume fibres from a point source release, showing the much greater distances travelled by the fibres. For example, for releases from Vienna, Austria, spherical particles with a diameter of 50 µm are deposited mainly in central Europe, whereas volume equivalent fibres with an aspect ratio of 50 are also deposited in eastern and southern Europe (Fig. 3a, b). The longer atmospheric transport of fibres results in an average horizontal travel distance that is 2.7 times greater than that of spherical particles (Fig. 3c).

In terms of technical implementation, parameters defining a particle's shape and orientation are to be provided in the SPECIES file. If the three characteristic dimensions (longest axis PLA, intermediate axis PIA, and smallest axis PSA) of a non-spherical particle (Blott and Pye, 2008) are known, the user should set PSHP=1 and provide the three-dimensional parameters in units of micrometres. If the particle dimensions are not fully known, FLEXPART can compute PLA, PIA, and PSA values for a few specific non-spherical particle shapes, assuming that the particles have the same equivalent volume as a sphere with the given diameter PDIA. For cylinders (PSHP=2), the user also needs to set the aspect ratio PASPR, i.e. ratio of the cylinder's length (PLA) to its diameter (PSA). Other options are cubes (PSHP=3), regular tetrahedrons (PSHP=4), regular octahedrons (PSHP=5), and a regular rotational ellipsoid (PSHP=6) characterised by PLA=2 ⋅ PIA=2 ⋅ PSA. PSA needs to be specified by the user. These options are particularly useful for a direct comparison of the atmospheric transport properties of non-spherical particles with those of corresponding spherical particles with the same volume.

The particle orientation in the atmosphere can be specified with the option PORIENT, with horizontal or maximum drag (PORIENT=0), random orientation (PORIENT=1), and the average of the two (PORIENT=2). If straight or bent cylindrical particles with aspect ratios equal to or greater than 20 are used, it is recommended to select the averaged orientation option, as this model best fits the experimental data on the gravitational settling of cylinders (Tatsii et al., 2024).

5 Wet scavenging

Wet scavenging is the process by which trace substances are removed from the atmosphere through precipitation. The mass of particles is reduced, and that amount is deposited to the ground. Three main changes were implemented to improve the wet-deposition scheme in FLEXPART 11. A full discussion with comprehensive tests will be provided in a forthcoming publication led by Anne Tipka.

Figure 4 Accumulated wet deposition (shading) 66 h after a point release of 1 kg of aerosols with a diameter of 0.4 µm and a density of 2×10³ kg m⁻³ for FLEXPART 10.4 (a), FLEXPART 11 with standard precipitation fields (b), and FLEXPART 11 using two additional precipitation fields (c). In total, 10 million particles were released at a constant rate between 18 January 1995 03:00 UTC and 19 January 1995 03:00 UTC from the location of the Ascó Nuclear Power Plant (41°12^′ N, 0°3^′ E; red star). Total and maximum deposition amounts are reported at the top of each panel.

5.1 New below-cloud scavenging scheme for aerosols

A major change with respect to wet scavenging was introduced in FLEXPART version 8, differentiating between in-cloud and below-cloud scavenging. For the below-cloud scavenging of aerosol particles, the scheme of Laakso et al. (2003) was used with parameters for snow, as given by Kyrö et al. (2009). This scheme was derived for small particles only, with diameters between 10 and 510 nm. For FLEXPART 11, a new parameterisation scheme has been introduced that is valid for a wider size range, namely the scheme of Wang et al. (2014), which covers particles from 1 nm (nucleation mode) to 100 µm (coarse mode). This is important especially for larger particles like dust, sea salt, or pollen.

5.2 Improved interpolation scheme for precipitation

Version 10 of FLEXPART used a nearest-neighbour method to obtain the precipitation rate and cloud information at the location of each particle at the given time; however, everywhere else in FLEXPART, linear interpolation is used. The nearest-neighbour scheme for precipitation results in unrealistic “checkerboard-like” deposition fields, which become visible when time intervals between meteorological fields are large, e.g. 3 h (Hittmeir et al., 2018), and especially in incremental deposition fields. The underlying problem, which led the authors of FLEXPART 10.4 to remove the linear interpolation of precipitation as implemented before, is related to the fact that the precipitation data are not given as an instantaneous rate but as an accumulation over a time interval (Hittmeir et al., 2018; Tipka et al., 2020). Hittmeir et al. (2018) therefore proposed a new reconstruction algorithm for Lagrangian particle dispersion models, like FLEXPART, which can be applied in conjunction with linear interpolation and which ensures the conservation of integral precipitation within each time interval, maintains continuity at interval boundaries, and conserves non-negativity.

In order to use this new interpolation scheme, the flag numpf in par_mod.f90 must be set to three. This causes FLEXPART to read in two sets of precipitation fields (a set of three large-scale and convective precipitation fields each, making it six fields in total) per input time interval instead of one. These additional fields represent two additional supporting time steps in between the original ones to represent precipitation as point values and, at the same time, preserve the integral precipitation in each time interval, guarantee continuity at interval boundaries, and maintain non-negativity (see Hittmeir et al., 2018, for more details). They can be produced with flex_extract version 7 (Tipka et al., 2020). If this is not set, precipitation will be linearly interpolated between the two precipitation fields assigned to the same times as the other quantities, which is still considered an improvement on the nearest-neighbour scheme used in FLEXPART 10.4. The example shown in Fig. 4 clearly illustrates the reduction in artefacts moving from the nearest-neighbour method to simple linear interpolation and finally the algorithm by Hittmeir et al. (2018) (numpf=3). In this example, 3-hourly input data were used, which make the artefacts more pronounced than with 1-hourly data.

5.3 New cloud layer identification scheme

The cloud identification scheme has been improved in FLEXPART 11. One problem is that the cloud water content is provided as instantaneous fields in the ECMWF meteorological input data, whereas precipitation is an accumulated forecast product. Their combination in the scavenging scheme can, at times, produce inconsistencies such as precipitation occurring when no (or only a very thin) cloud is present, and this can lead to unrealistic in-cloud scavenging rates in FLEXPART. To prevent such unrealistic combinations, the following steps are taken to identify clouds. If cloud water content is available from the meteorological data, the cloud water mass is integrated vertically and used later on for the in-cloud scavenging, as already done in version 10. Iterating from the surface upwards, the first instance of cloud water is taken as the height of the bottom of the cloud and the last instance as the top of the cloud, as in version 10. If, however, the cloud water content is not available, the presence of a cloud is assumed between the lowest level, with relative humidity greater than rhmin, and the highest level with that value. Relative humidity is calculated over water for temperatures above −20 °C and over ice for temperatures below. The value of rhmin is set in par_mod.f90 (default is 0.90). A minimum thickness of min_cloudthck (defined in par_mod.f90; default 50 m) is required for a cloud. Spurious clouds below this thickness are eliminated, as they would not be expected to produce precipitation, but if precipitation were present in the data, it would cause unrealistic in-cloud scavenging rates.

For the case of convective precipitation, a simple fix was implemented; it requires that clouds associated with convective precipitation stronger than precmin (in par_mod.f90; default is 0.002 mm h⁻¹) must have a top of at least 6 km and a bottom below 3 km (adjustable by conv_clrange in par_mod.f90). If this is not fulfilled, the cloud bottom and top will be set to a bottom and top of 0.5 and 8 km for precipitation rates below 0.1 mm h⁻¹, and to 0 and 10 km for heavier precipitation (parameters are set in par_mod.f90 to highconvp_clrange and lowconvp_clrange). The results of the preliminary tests indicate that the deposition resulting from convective clouds is not significantly influenced by the parameters used in this fix. However, further investigation is required to ascertain the full extent of the influence of these parameters and possible further optimisation.

5.4 Fukushima test

To validate that changes made to the precipitation scheme did not introduce errors, several tests have been carried out to check the behaviour of the code with respect to the removal processes. Here, we show a series of tests reproducing the Fukushima Daiichi Nuclear Power Plant accident of March 2011, where the aerosol-bound radioactive isotope caesium-137 (¹³⁷Cs) and the noble gas xenon-133 (¹³³Xe) were released in large quantities. As demonstrated by Kristiansen et al. (2016), after correcting for radioactive decay, the ratio of caesium-137 and xenon-133 can be used to evaluate the modelled aerosol lifetimes. They obtained aerosol lifetimes for 19 widely used chemical transport models (including FLEXPART 9) and compared the model results with the aerosol lifetime obtained from experimental measurements. Their results showed, for the aerosol e-folding lifetime (among all the models), a model mean of 10.7 d and a model median of 9.4 d. FLEXPART 9 resulted in a shorter aerosol lifetime of 5.8 d, using the original Hertel et al. (1995) definition of the cloud liquid water content, c_l (see Sect. A3.2). The observations, as discussed in Kristiansen et al. (2016), suggested a rather longer aerosol lifetime of 14.3 d. More recently, Pisso et al. (2019) re-evaluated the aerosol lifetime in FLEXPART 10.4 with the Grythe et al. (2017) improved scheme and obtained an aerosol e-folding lifetime of 10 d for caesium-137. Figure 5 shows the aerosol lifetime obtained for three different FLEXPART configurations using different values of ratio_incloud, which is the in-cloud replenishment rate r_icl (see Appendix A3.2 for a full description). An e-folding lifetime of 9.3 d is obtained for FLEXPART 11 when using ratio_incloud=0.005, very close to the previously obtained model median by Kristiansen et al. (2016), and therefore the default value set in par_mod.f90. A value of ratio_incloud=0.001 (results not shown here) gave a e-folding lifetime of about 20 d.

Figure 5 Global decrease in the ratio of caesium-137 with respect to xenon-133 (decay corrected) as obtained at several measurement stations as a function of time after the Fukushima nuclear reactor accident. Panels (a), (b), and (c) show the results using values of 0.01, 0.0062, and 0.005, respectively, for the replenishment rate r_icl of cloud water during precipitation (parameter ratio_incloud in FLEXPART 11). Corresponding e-folding lifetimes of the aerosol were 6.99, 8.58, and 9.32 d, respectively. Red boxes show the ratio of the sum of caesium-137 and the sum of xenon-133 over all stations.

6 Linear Chemistry Module

The Linear Chemistry Module (LCM) is based on the initial work of Henne et al. (2018), who developed the FLEXPART CTM from FLEXPART 8, and which was first described and evaluated in Groot Zwaaftink et al. (2018). This model was an extension of the domain-filling capability of FLEXPART and added the possibility of initialising particles' mixing ratio from pre-defined fields, accounting for the influence of surface fluxes and simple linear chemistry on the particles' mass, and sampling the particle-mixing ratios at user-defined receptor locations.

The essence of the FLEXPART CTM code has been implemented into FLEXPART 11, which can be run in all its usual configurations but now includes the CTM configuration. Note that the CTM configuration is renamed the Linear Chemistry Module, reflecting the fact that it is not a separate model to FLEXPART and that only linear chemical reactions are implemented.

LCM uses the possibility in FLEXPART to fill a global domain with particles that are constantly advected throughout the atmosphere. This scheme was introduced by Stohl and James (2004) for air tracers and has been used for studies of heat and water transport in the atmosphere. However, each particle can carry multiple chemical species. Emissions of a species are not simulated as releases of new particles representing the emitted mass, as it would be done in a standard FLEXPART run, but are taken up by particles passing by close to the emission source (note that emission in this sense can be positive or negative fluxes; e.g. carbon dioxide could be removed from the atmosphere), thereby changing the mass of the species carried by these particles. The mass of the particles for each species can also change due to first-order chemical reactions. Deposition of particles is at present not considered.

In the LCM, the mixing ratios are output by sampling the particles on a grid (specified by OUTGRID) at receptor locations (specified by RECEPTORS and, new in FLEXPART 11, for satellite receptors (specified by NetCDF files generated from a satellite pre-processor, prep_satellite, and can be created using software obtainable from https://git.nilu.no/flexpart/flexinvertplus, last access: 24 October 2024). The mixing ratios for each species are calculated using the ratio of the mass of the species over the mass of air, where the mass of air is always carried by species number 1. This method was chosen over using the mass of the species divided by the air density and volume, which is the standard method because it leads to less noisy results and avoids problems with the spurious accumulation of particles (see Sect. 3.3).

The LCM mode is switched on in FLEXPART in COMMAND with the switch LCMMODE=1. In addition, the following options in COMMAND should be set: MDOMAINFILL=1, IND_SOURCE=1, and IND_RECEPTOR=1. When these options are set, mixing ratios will be calculated using the ratio of masses method. The domain-filling mode uses the RELEASES file to define the species, the domain, and the total number of particles to use. For LCM, the first species in RELEASES must always be AIRTRACER.

A full discussion with tests and examples will be provided in a forthcoming publication led by Rona Thompson. Nudging to observations, as implemented in FLEXPART CTM (Groot Zwaaftink et al., 2018), is not yet available in the LCM module but will be added in the near future.

7 Comparison with tracer experiments

In Sect. 3, we have shown how quasi-conservative properties of particles are better conserved, better mixed, and better fulfilled when using η instead of z vertical coordinates. However, we also made large modifications in other parts of the source code, and therefore it was important to validate FLEXPART 11 for a large range of different applications, including output of gridded properties. One excellent validation possibility is to compare the model against tracer experiments. We chose to follow Stohl and Koffi (1998) and applied FLEXPART 11 to the first European Tracer Experiment (ETEX) (Nodop et al., 1998) and the Cross APpalachian Tracer EXperiment (CAPTEX) (Raynor et al., 1984; Ferber et al., 1986). Following Stohl and Koffi (1998), we compute the following:

• the normalised mean square error, $\text{NMSE}=(\mathrm{1}/N){\sum }_{i=\mathrm{1}}^N(P_i-M_i{)}^{\mathrm{2}}/(\stackrel{\mathrm{‾}}{P}\cdot \stackrel{\mathrm{‾}}{M})$, where P and M represent the model predictions and measurements, respectively, and N is the number of measurements;

• the root mean square error, $\text{RMSE}=\sqrt{(\mathrm{1}/N){\sum }_{i=\mathrm{1}}^N(P_i-M_i{)}^{\mathrm{2}}}$, where P and M represent the model predictions and measurements, respectively, and N is the number of measurements;

• the Spearman rank-ordered correlation coefficient (SCC);

• the Pearson correlation coefficient (PCC);

• the figure of merit in space, $\text{FMS}=\mathrm{100}{A}_{\text{p}}\cap A_{\text{m}}/A_{\text{p}}\cup A_{\text{m}}$, where A_p and A_m represent the model-predicted and model-measured concentrations above 0.1 ng m⁻³, respectively;

• the fraction of model predictions that lie within 2 times (FA2) and 5 times (FA5) the measured values; and

• the frequency of over- and underpredictions, $\text{FOEX}=\mathrm{100}(N_{P_i>M_i}/N-\mathrm{0.5})$, with $N_{P_i>M_i}$ as the number of overpredictions.

We computed the statistical measures for two different settings of FLEXPART 11, using the η coordinates and z coordinates. To ensure consistency with the previous FLEXPART version, we also calculated the statistics for FLEXPART 10.4, using the same input options as for FLEXPART 11. We used ERA5 meteorological input data, simulated a passive tracer and used time steps shorter than 10 % of the Lagrangian turbulent timescale (settings CTL=10 and IFINE=10). We use an output grid with 0.1° × 0.1° resolution.

Figure 6 Comparison of the concentration fields obtained by FLEXPART 11 (blue contours), using η coordinates on a 0.1°×0.1° grid, and the ETEX measurements (markers) 48 h after the ETEX-1 release. Small grey circles represent each zero measurement that corresponds to a zero model result. Large grey circles represent the stations with non-zero measurements at which the model also found a signal. Indigo triangles represent stations with non-zero measurements at which the model found no significant signal, and indigo crosses show where the opposite is the case. Shading within the circles and triangles shows the concentrations measured at the corresponding stations. The initial release location near Rennes is marked by a red star. Note that results using z coordinates and FLEXPART 10.4 do not show significant differences.

7.1 ETEX

On 23 October 1994, 340 kg perfluoromethylcyclohexane (PMCH) was emitted near Rennes, France. The release lasted 12 h, and measurements were taken for the next 90 h across 168 stations distributed all over Europe (see Fig. 6). We recreated this experiment by releasing 1 million particles equally spread over the 12 h time frame with altitudes between 8 and 20 m above ground level.

Table 1 Statistical measures of how well FLEXPART 10.4 and FLEXPART 11, using η- and z-coordinate systems, respectively, perform compared to the ETEX observations (top) and compared to the CAPTEX aircraft measurements (bottom). CAPTEX results show the mean statistics over six releases. n gives the total number of model–observation pairs used for the analysis. The RMSE is given in $\mathrm{ng}{\mathrm{m}}^{-\mathrm{3}}.$

The results of the statistical measures are shown in Table 1. The versions using the z- and η-coordinate systems give almost identical results. This is not unexpected, as most of the gas stays in the turbulent ABL, which is where the z-coordinate system is always used. Importantly, FLEXPART 11, using the z-coordinate system version, produced similar results as FLEXPART 10.4, which demonstrates that the many code changes did not result in undesired performance losses.

7.2 CAPTEX

During the Cross APpalachian Tracer EXperiment (CAPTEX), seven separate PMCH releases were made in North America (Ferber et al., 1986), and both ground-based and airborne PMCH concentration measurements were made. Even though all aircraft measurements were made below 3053 m a.s.l., and many within the ABL, the values at these altitudes are more likely to be affected by the changes in the coordinate system than those at the ground level; therefore, we present FLEXPART results only for the aircraft measurements. To compare the difference in densities between the aircraft measurements and the FLEXPART output, we use the mass of particles present in a volume of a 0.5° by 0.5° by 100 m box around the measurement locations.

In Table 1, we list the means of the statistical measures for all CAPTEX releases, excluding release 6, for which the number of valid observations was very small. We show the mean of the statistical indicators, instead of highlighting the spread in accuracy between the different releases. Table 1 mainly shows that using the η-coordinate system as compared to the z-coordinate version leads to no substantial differences. We also see no large systematic differences between FLEXPART 10.4 and FLEXPART 11.

8 OpenMP parallelisation

Pisso et al. (2019) implemented a pure Message Passing Interface (MPI) method for the parallelisation of FLEXPART. MPI often performs equally well or better compared to other multi-processing options without excessive communication between processes and potentially low overhead. However, FLEXPART runs on input meteorological data that is ever growing in size, and without domain decomposition, each MPI process needs to keep a copy of the same meteorological data in memory. This leads to a large memory footprint that can easily exhaust the memory of a compute node. Consequently, we have followed a different parallelisation strategy for FLEXPART 11 that is based on OpenMP, which shares memory between processes. While the OpenMP option limits the number of parallel processes to the number of cores available on a compute node, this is not a significant limitation for typical FLEXPART applications. Furthermore, because of the linear nature of FLEXPART calculations, a trivial parameterisation by splitting a simulation into several separate run instances, each running on a different node, is always possible for particularly large simulations. One drawback of OpenMP compared to serial code is that it is more difficult to modify. To avoid errors related to parallel execution, users unfamiliar with OpenMP are strongly advised to make changes in the form of subroutines and functions and avoid the use of global variables.

8.1 Parallelisation strategy

The way in which the computational time used in a FLEXPART run is partitioned between different parts of the code depends on the set-up and thus on the peculiarities of the application considered. For example, if many particles are released from a single point and only a small computational domain is used, then initially most of the time is spent on particle trajectory computations. However, when a global high-resolution domain for the meteorological input data is used in a global domain-filling simulation, a significant share of the CPU time is spent on the convection parameterisation (see Appendix A2.1). On the other hand, if relatively few particles are used, then computations related to the gridded meteorological input data (e.g. coordinate transformations) are taking a larger share. For this reason, the OpenMP parallelisation was implemented throughout the various components of the model, trying to avoid bottlenecks for all these cases, especially for the most common set-ups.

We parallelised all particle-based computations, apart from their initial release in the releaseparticles subroutine. On top of that, we parallelised the convection, wet and dry deposition, and the vertical coordinate transformation of the fields. Last, we parallelised the computations needed for the output for both the gridded output and the particle dump.

The total number of threads for gridded-output-related computations can be set by the user (MAXTHREADGRID$=<$number of threads> in the COMMAND file) with a maximum of the general requested number of threads. The reason for this is that parallelisation is applied to the particles, and therefore each thread needs its own set of output grid variables. Depending on the resolution of the output grid, this can significantly increase the memory required. On top of that, depending on the number of particles in the simulation, the related computational overhead might outweigh the speed-up. The default of MAXTHREADGRID is set to a maximum of 1 task, thus switching off parallel computations related to the grid. A number of threads >1 should only be set after having tested with the specific set-up to be used. For example, for a simulation with 10 million particles and a 0.5° by 0.5° global output grid, using more than 16 threads will degrade the performance (with the optimum number being 6), while for a simulation with only 10 000 particles and a 0.1° by 0.1° output grid, parallelisation is not at all efficient, and MAXTHREADGRID should be left at its default value of 1.

8.2 Code performance and scalability

To evaluate the scalability of the OpenMP parallelisation across different FLEXPART running options, we select a few test cases across the spectrum of typical FLEXPART applications.

• Case Tracer. This is a domain-filling simulation with particles representing a passive air tracer distributed across the globe following air density. Every hour, all particle information but no gridded output is written to NetCDF files (see Sect. 9.1). We run FLEXPART for 5 h, using 10 min time steps, and we set the turbulence options CTL=10 and IFINE=10. This option covers a large range of different computations within FLEXPART (i.e. convection and turbulence in all possible conditions and particle output) and is therefore a difficult condition to obtain perfect scalability for.

• Case Aerosol. This involves aerosol particles that are initially homogeneously distributed in the bottom 100 m across the globe. Every hour, gridded properties are written to NetCDF files on the same horizontal resolution as the input data (0.5° by 0.5° global grid) and cover four vertical levels. FLEXPART is run for 5 h, using 15 min time steps and turbulence options CTL=10 and IFINE=4. Case Aerosol simulations generally take much longer than Case Tracer simulations. On the one hand, this is because of the extra computations in the wet and dry deposition and gravitational settling routines. On the other hand, this is because all particles start within the ABL, where solving the Langevin equations of the turbulence parameterisation requires very short time steps. We set MAXTHREADGRID to 16, meaning that the gridded computations are using a maximum of 16 threads.

• Case Nuclear. A single xenon-133 point release, using the CBL option for the skewed turbulence (see Appendix A2.2) and CTL=40 and IFINE=5, with a nested input domain and nested gridded output over Europe with 0.25° × 0.25° resolution. Unlike with MPI parallelisation, OpenMP parallelisation is only active within specified regions of the code. With this third case, we want to complete our demonstration of the parallelisation of all possible parts of FLEXPART, including radioactive decay, skewed turbulence, and computations on the nested grid. Gridded computations are conducted on a single thread.

We use two different computers to demonstrate the performance of FLEXPART 11 across platforms. The first one is the local department server (called “Jet”), which has 9 Intel Skylake compute nodes, each with 2×20 cores and 80 threads in total and 768 GB memory. The second one is the Vienna Scientific Cluster 5 (VSC-5), which has 564 AMD EPYC Milan compute nodes, each with 2×64 cores and 256 threads in total, 512 GB RAM, and 1.92 TB solid-state drive (SSD) storage. For both computers, we use the GNU FORTRAN compiler with the optimisation flag -O3, and since we have Intel central processing units (CPUs) on Jet, we use -march=skylake-avx512 for the compiler to correctly recognise the hardware architecture. The compiler versions are the GNU compiler collection (GCC) 8.5.0 on Jet and GCC 12.2.0 on the VSC-5. We set export OMP_PLACES=cores and export OMP_PROC_BIND=true for the best use of non-uniform memory access (NUMA) and therefore the fastest performance. To make the comparison between the two versions clearer, we do not make use of the hyperthreading capabilities of the computers. This means that the number of OpenMP threads is equal to the number of tasks used by MPI, and we use threads and cores interchangeably when referring to OpenMP processes.

Figure 7 Demonstration of the strong scalability of the OpenMP parallelisation of FLEXPART 11. Shown is the computation time as a function of the number of cores used (N_core). Perfect strong scalability is represented by the slope of the straight and thin dotted black lines. The left column of the panels, shaded in blue, shows the results for globally distributed passive tracer particles (Case Tracer), using 10⁶ (a1), 10⁷ (a2), and 10⁸ (a3) particles. The middle column, shaded in magenta, shows results for globally distributed aerosols (Case Aerosol), using 10⁵ (b1), 10⁶ (b2), and 10⁷ (b3) particles. The right column, shaded in maroon, shows results for xenon-133 emitted from a single point source (Case Nuclear), using 10⁵ (c1), 10⁶ (c2), and 10⁷ (c3) particles. In black, the total of 5 h running time is plotted for FLEXPART 11 (solid black lines). The MPI-parallelised FLEXPART 10.4 equivalent is plotted for reference (solid grey lines). The dashed coloured lines represent the break-down of the computational time spent in different components of FLEXPART 11; the blue lines represent the total computational time per hourly time step minus the input/output (I/O) operations (green lines) and initialisation of the run (orange lines). Note that for some cases, green lines cannot be seen, as the I/O operations fall below the limits of the y axis but are of a similar order to the I/O operations of the same case with fewer particles. The runs represented by the thick lines and the FLEXPART 10.4 run were executed on the VSC-5 machine, and the thin dashed–dotted black lines show the total FLEXPART 11 running time on our local server, Jet. Note that the black lines are not the sum of all other lines, since the black lines show the total time that the simulation took, while the I/O operations and computations are shown per hour of simulation time. This was done for clarity of presentation, i.e. to better separate the lines from each other.

Figure 8 The memory footprint of FLEXPART 11 (solid lines for z; dotted lines for η-coordinate system) and FLEXPART 10.4 (dashed lines) simulations as a function of the number of cores. Shown are the results for Case Tracer, a domain-filling simulation using 1 million (thick purple), 10 million (orchid), and 100 million (thin grey) particles. Note that the FLEXPART 11 cases using 1 and 10 million particles largely overlap.

Table 2 Weak scaling and serial comparison of different applications described in Sect. 8.2. All values are given in percentages.

We investigate both the strong and weak scalability of our OpenMP parallelisation for each case. Strong scaling refers to the speed-up by increasing the number of threads. Ideally, the simulation time should scale close to the inverse number of cores, which is $\mathrm{1}/N_{\text{cores}}$ for all cases (e.g. ideally giving a speed-up of 10 when using 10 times the number of cores). We separated the initialisation, which only happens once and would thus become relatively less important for longer model runs than our rather short 5 h test cases, and input/output (I/O) operations, which are very dependent on the output needs of the user (e.g. output frequency of the particle dump), from the computations that are done every time step, which probably best represent the overall model performance in realistic cases. Weak scaling refers to the relative speed when increasing the size of the problem in proportion with increasing the number of cores used for the computation (e.g. running 10 times the problem size on 10 times the computational resources should ideally take the same amount or less time). In Table 2, this is defined by the time it takes to run A particles on 1 core compared to running B particles on 10 cores, i.e. $\left[t\right(B_{\mathrm{10}})-t(A_{\mathrm{1}}\left)\right]/t\left(A_{\mathrm{1}}\right)$. In addition, we provide a serial comparison, which is done by comparing the time it takes to run A particles and run a single B particle 10 times, both on one core, i.e. $\left[t\right(B_{\mathrm{1}})-\mathrm{10}t(A_{\mathrm{1}}\left)\right]/\mathrm{10}t\left(A_{\mathrm{1}}\right)$. In our case, we only increase the number of particles but keep the size of the meteorological input data constant for all parallel runs, since most users do not switch between data types. That means that the problem size does not entirely increase 10-fold when running 10 times the number of particles.

Scalability of computation times is not the only issue to consider for evaluating the FLEXPART performance. Equally important is the memory requirement of the model. Therefore, we also compare the memory footprint of FLEXPART 11 with the MPI-parallelised version of FLEXPART 10.4.

8.2.1 Case Tracer

The left column of Fig. 7 shows how the strong scalability for Case Tracer depends on the number of particles used in a simulation by benchmarking FLEXPART 11 using 1 million (top), 10 million (middle), and 100 million particles (bottom). Every panel also shows a comparison to the FLEXPART 10.4 MPI version. For 100 million particles, FLEXPART 10.4 with MPI and FLEXPART 11 are performing similarly, while for smaller particle numbers, FLEXPART 11 clearly scales better than FLEXPART 10.4. The scalability of the main computations (blue lines) is dependent on the number of particles that are used; for example, on VSC-5, when using 32 cores, we find a ∼1.8-, ∼2.5-, and ∼2.8-fold increase in computational time compared to the perfect scaling for the three problem sizes of 1 million, 10 million, and 100 million particles, respectively. On Jet, the performance is worse, with ∼2.3, ∼3.0, and ∼4.1, for the three problem sizes, respectively. The large difference between VSC-5 and Jet is likely the result of using a more recent version of the GCC compiler on VSC-5 compared to Jet (GCC 12.2.0 instead of GCC 8.5.0), which optimises the non-parallelised regions more effectively. The time it takes to write the NetCDF particle dump files is less than 10 % of all the computations per step when using only a few cores. This increases to up to 43.1 % when using 128 cores because the writing of NetCDF files cannot be parallelised when using OpenMP. Since the output size for the Case Tracer increases linearly with the number of particles used, the total scalability suffers from the relatively inefficient writing of output files when using many cores.

The weak scalability depends heavily on the number of particles of the simulation. We find a decrease of more than 50 % in the computation time when running 10 million particles on 10 cores instead of running 1 million particles on 1 core (see the top row of Table 2), which can be explained by the parallelisation of the convection computations, since these take an (almost) equal amount of time, regardless of the number of particles in a grid column. Increasing the number of particles from 10 million to 100 million, however, increases the computation time by more than a third. Profiling indicates the likely cause to be inefficient communication between NUMA regions. Therefore, for large problem sizes, it could be worthwhile to parallelise by splitting the problem into multiple smaller ones that are simulated separately. However, fewer CPU hours in total will be consumed when simulating more particles at once on a single core; the serial comparison in Table 2 shows that in the range of 1 million to 100 million particles, it always takes less time to simulate more particles at once than splitting them up into multiple smaller runs.

8.2.2 Case Aerosol

For Case Aerosol (middle column of Fig. 7), we chose to use 100 000, 1 million, and 10 million particles for our benchmark testing, separating the computation in the same way as described in Case Tracer. Generally, the performance of FLEXPART 11 is similar to or better than that of FLEXPART 10.4, with the exception of the serial version, where FLEXPART 10.4 is faster. Large improvements in performance for version 11 are especially visible when using smaller numbers of particles on more cores, where the MPI version suffers from more overhead.

The main computations (without I/O and initialisation) of Case Aerosol generally have better scalability than Case Tracer, coming close to perfect scaling for all particle number cases, especially when using up to eight cores. For the example of 32 cores, scaling is only 1.7 (2.5), 1.4 (2.0), and 1.4 (1.9) times slower for 100 000, 1 million, and 10 million particles on VSC-5 (Jet), respectively, than perfectly scaled code would be. The I/O per step, reading meteorological data, and the organisation and writing of information to the grid in NetCDF format are only improving slightly with the number of cores, but there is a significant bottleneck only for the smallest problem size (100 000 particles). Here, it takes up ∼4 % of the time of all computations in serial mode and ∼57 % when using all 128 cores as a consequence of the excellent scaling of the other computations. The largest problem size reports a minimum of ∼0.2 % (serial mode) and a maximum of ∼8 % (128 cores) time being spent on I/O operations.

The weak scalability of the main computations is adequate, as shown in the middle row of Table 2, taking about a fifth less time to compute a time step of 1 million particles on 10 cores as compared to 100 000 particles on 1 core but taking slightly more time for the larger problem size. When running FLEXPART in a serial mode, fewer CPU hours are consumed by simulating many particles at once compared to splitting the problem into multiple smaller simulations, although this gain decreases when moving towards larger problem sizes.

8.2.3 Case Nuclear

For Case Nuclear (right column of Fig. 7), we chose a single release of 100 000, 1 million, and 10 million particles. Similar to Case Aerosol, the strong scaling of the main computations for 1 million and 10 million particles is reasonably good and only 1.3 (1.4) and 1.3 (1.2) times slower, respectively, than perfect scaling when using 32 cores on VSC-5 (Jet). For the 100 000 particle run, this becomes 2.1 (2.9) times slower. For all problem sizes and number of cores, FLEXPART 11 outperforms FLEXPART 10.4. The I/O computations take no more than ∼3.4 % of all computations per step in serial mode and a maximum of ∼36 % when using the full 128 cores for the smallest problem size. For the largest problem size, this ranges between ∼0.08 % (serial) and ∼8.2 % (128 cores).

The weak scaling relation is within the range of the other two cases. Table 2 shows that a run using 1 million particles divided on 10 cores takes approximately the same time as simulating 100 000 particles on 1 core. And a run using 10 million particles divided on 10 cores takes between ∼10 %–15 % more time than a run of 1 million particles on only 1 core. However, when running FLEXPART on a single core, fewer CPU hours are consumed by simulating many particles at once compared to splitting the problem into multiple smaller simulations, although this gain is marginal when requiring a total number of more than 10 million particles.

8.2.4 Memory requirements

In the previous section, we showed that in most cases, OpenMP parallelisation results in better speed-up than the MPI parallelisation of FLEXPART 10.4. However, the main advantage of using OpenMP over MPI is that data (especially the meteorological input fields) are shared between the processes instead of requiring separate copies on each core. This reduces the memory footprint substantially for runs on shared-memory systems. Typically, work group servers and most high-performance computers today provide many cores per node but only a limited amount (some GB) of RAM per core. For example, as shown in Fig. 8, when running FLEXPART 11 on Case Tracer with 10⁶ particles, approximately 11 GB of memory are used, regardless of the number of cores used. An equivalent simulation with FLEXPART 10.4 on 32 cores, requires about 132 GB – more than 10 times the memory footprint of FLEXPART 11. Running FLEXPART 10.4 on one full node with 128 cores on the VSC-5 exceeds the node's available memory of 500 GB. We note that without parallelisation, FLEXPART 11 needs somewhat more memory than FLEXPART 10.4 (e.g. 11 GB instead of 7 GB for the example above), which is mainly due to the extra memory needed for the η-coordinate system, moving to a necessary double-precision default for some variables in FLEXPART 11, and added functionality.

9 Running FLEXPART: new input and output options

In this section, we give an overview of new user functionalities of FLEXPART 11. For general information about how to run FLEXPART, we strongly recommend users to check the instruction manual (https://flexpart.img.univie.ac.at/docs, last access: 11 September 2024), since this page will be updated with any future additions and changes.

9.1 Particle output

In addition to the gridded output, FLEXPART can also write out information related to each particle, its position, and meteorological information associated with it. In FLEXPART 10.4, this was possible already, and for reasons of efficiency, unformatted binary files per output time step were produced. Since users found it cumbersome to work with these, for FLEXPART 11, particle information is now written to NetCDF files instead. With HDF5-based compression, NetCDF files are about ∼50 %–60 % smaller than their binary equivalents. Another advantage of the NetCDF output is that it is easier to modify the content or structure of the information to be written, especially with respect to safe post-processing. Taking advantage of this, it is now possible to specify which particle properties are to be written out. The desired parameters can now be specified in a new option file, PARTOPTIONS. It is now required if particle properties are to be written out (i.e. when IPOUT=1 is set in the COMMAND file). Currently, the available variables are particle position (longitude, latitude, and height in metres above ground), PV, specific humidity, air density, temperature, pressure, particle mass, cumulative loss of mass by dry and wet deposition (separately), settling velocity, 3D grid-scale particle velocities, the height of the ABL, the tropopause, and the topography. Each property can be written out as an average over the preceding output time interval and/or as instantaneous values.

In previous FLEXPART versions, the memory location of terminated particles was reused for newly released particles in order to minimise the memory requirement. This, however, complicates the tracking of individual particles. Hence, in FLEXPART 11, if particle output is switched on (IPOUT=1), terminated particles are kept in the simulation, but values associated with them are set to“NaN” (not a number) instead of being overwritten by newly released particles in the NetCDF output. This comes with no additional computational cost, but it may require more memory than running without the particle output option switched on. The previous behaviour of overwriting terminated particles when using particle output can be restored, as explained in the instruction manual (https://flexpart.img.univie.ac.at/docs, last access: 11 September 2024), and will be needed for domain-filling simulations with a limited domain, where particles are continuously produced at the inflowing boundary and destroyed at the outflowing boundary.

9.2 Starting a simulation from user-defined particle data

To give users complete control over the initial conditions of a simulation, in FLEXPART 11 we implemented an option to start a simulation from particle information provided in an additional input file in the NetCDF format. This file must hold information about the initial particle positions, their time of release, and their initial mass (see the manual for a full description). This option is used if IPIN=3 is set (IPIN=4 in the case of a simulation to be restarted). The user then needs to provide a file part_ic.nc in the options directory, as specified in the pathnames file. The RELEASES file in the input directory is ignored in this case. This option is particularly useful if release geometries are complicated and latitude–longitude–altitude boxes are not suitable. In Bucci et al. (2024), this option was used to produce a global release at the ocean surface. Other possible applications are FLEXPART runs for satellite observations, with non-trivial geometry of the pixel at the Earth's surface and height-varying kernel sensitivities. In such a case, for a backward simulation to determine the sensitivities of column-averaged molar fractions to emissions and/or initial conditions, one would fill the volume relevant for a specific satellite measurement with particles, horizontally homogeneous, and with a vertical profile proportional to the pressure-interval-weighted satellite retrieval kernel. Figure 9 shows such an example of a satellite observation over India.

Figure 9 Methane emission sensitivity produced with a 7 d FLEXPART backward simulation from a TROPOspheric Monitoring Instrument (TROPOMI) methane remote-sensing observation (Schneising et al., 2023). (a) Emission sensitivity obtained for the location of the TROPOMI observation. The inset (b) shows all TROPOMI methane observations that were made in the region inside the blue square, with one single observation near the centre highlighted by a red outline. The FLEXPART backward simulation was done for this single observation. (c) Vertical profile of the FLEXPART particle density per pressure interval (given as number of particles per hPa). The number of particles, P_n, released in layer n of the 20 TROPOMI retrieval layers was calculated as $P_{\text{n}}=P\cdot A_{\text{n}}\cdot W_{\text{n}}$, where P is the total number of particles used for the simulation (1 million), A_n is the averaging kernel value for retrieval layer n, and W_n is the pressure weight. The particles are uniformly distributed horizontally, based on the geometry of the satellite data pixel (i.e. the pixel marked with a red outline in panel (b)). Also shown in panel (c) is the averaging kernel and the air density. FLEXPART particle releases were made using the user-defined initial particle conditions option (file part_ic.nc).

9.3 Restarting a simulation

As a FLEXPART run may be interrupted (e.g. because of server maintenance or reaching the maximum allowed wall-clock time), already in FLEXPART 10.4 there was an option to restart a simulation from a saved simulation state. However, this required the user to request a particle dump to the binary files, which was done at every output time step; this could substantially slow down the whole simulation. Furthermore, the gridded deposition data lacked the information from the previous run. In FLEXPART 11, we created a separate option in the COMMAND file which sets the time interval for the output of the binary restart files (LOUTRESTART). The content of these files has been amended so that the restarted simulation will receive all the relevant information. For restarting, set IPIN=1 in the COMMAND file and rename or link the desired start point file from restart_XXX file to restart.bin in the output directory. It is recommended to set LOUTRESTART to a relatively large value in order to minimise the computational and storage overhead.

9.4 Dynamical allocation of arrays

In FLEXPART 11 it is no longer necessary to specify dimensions of (nested) input fields, receptors, or particle- or species-related arrays at the compilation time by manually setting the dimensions in the par_mod.f90 file. These arrays are now dynamically allocated as required by the size of the input data, enhancing the user-friendliness. Another advantage is that it allows us to containerise FLEXPART.

9.5 Updated SPECIES files

Together with the model code, a set of input files describing the physical–chemical behaviour for various species is provided as templates (SPECIES_xxx, with “xxx” being a number in the model input directory); these properties can and should be set by the users and are at their own discretion. Each file is dedicated to one “species” which represents a real-world species or group of species with a specific behaviour in the model. The information contained in that file is used to calculate wet and dry deposition, gravitational settling (radioactive) decay, or destruction due to the reaction with the hydroxyl radical, OH (OH fields have to be prescribed). The parameter values in the template files provided contain suggested parameter values; short references (name of the first author, year of publication, or a web link) are provided to the sources from which the parameter values were taken. A more detailed description is given in Appendix B.

10 Platform for interaction of users and developers

FLEXPART is a community model, allowing everyone to modify the code according to their specific needs, to fix bugs, to add functionality, and to communicate with the current developers. Since its inception, it has been available under the General Public License (GPL) V3 as free and open software. A web-based project management and bug-tracking system is available for use by developers and for interaction with users. With FLEXPART 11, this has been switched from a trac hosted by the Central Institute for Meteorology and Geodynamics (now GeoSphere Austria) to a GitLab instance hosted at University of Vienna. It includes a continuous integration (CI) functionality. The University of Vienna provides this open GitLab instance within the Phaidra (Permanent Hosting, Archiving and Indexing of Digital Resources and Assets) service for long-term research data storage (PHAIDRA, 2008). Currently, a few simple tests are bundled with FLEXPART, which are used in the CI. They should be seen as a starting point for future enhancements. The continuous code testing supports community code development and provides additional information on the deployment of FLEXPART under various scenarios (e.g. compilers and dependencies). It further allows us to provide a FLEXPART container that is ready for deployment and easy porting of FLEXPART to various high-performance computing environments or cloud instances.

11 Conclusions

We have presented the new version 11 of FLEXPART. A large number of development steps have been carried out which improve the accuracy, enhance computational properties, and add functionality.

The following new features improve the accuracy of FLEXPART simulations:

1. The possibility of performing transport calculations in the native η vertical coordinate system for ECMWF data, instead of internally converting all fields to vertical coordinates in metres, results in reduced absolute transport conservation errors in PV and other quasi-conservative variables and a better, but still not perfect, fulfilment of the well-mixed criterion in global domain-filling simulations.

2. An improved parameterisation scheme for the wet scavenging of aerosols below clouds, an improved precipitation interpolation method, and a more transparent cloud layer identification method were introduced, resulting in the removal of artefacts that were present in previous model versions.

3. More accurate drag coefficients for aerosol particles of different sizes and shapes, which can be characterised in all three dimensions, in some cases dramatically improve the simulation of the atmospheric transport of coarse-mode aerosols.

4. The implementation of the optional linear chemistry module turns FLEXPART into a simple linear chemistry transport model. When activated, it is possible to initialise particle-mixing ratios, incorporate the influence of surface fluxes and linear chemistry, and sample particle-mixing ratios at receptor locations.

In terms of computational enhancements, the switch from MPI to OpenMP parallelisation brings better memory utilisation and improved scalability across a wide range of applications. However, we note that significant improvements could still be achieved using co-arrays or a hybrid OpenMP–MPI implementation, especially in combination with a domain decomposition of the meteorological input data.

Finally, new use options were introduced in FLEXPART 11, and the SPECIES template files have now been documented. The option to replace the conventional particle release information (fixed release rates for a list of given spatiotemporal domains) with arbitrary initial particle conditions adds flexibility and, for certain cases, efficiency. The option to restart a FLEXPART simulation after an unforeseen or planned termination has been made more complete and more computationally efficient. Particle information can now be dumped in a configurable manner in a compressed NetCDF format.

Finally, to aid the user community, an instruction manual (https://flexpart.img.univie.ac.at/docs, last access: 11 September 2024) has been written as a living document. A continuous integration environment based on GitLab is available for current developers and those who wish to contribute to the further development of FLEXPART.

Appendix A: FLEXPART overview and code reorganisation

In FLEXPART, the modelling of computational particle movement involves a combination of interpolating properties from meteorological input data and computing properties that are intrinsic to each particle. This hybrid approach allows for a more comprehensive and realistic representation of particle behaviour in the atmosphere. Computational particles can represent one or more species, which may be gases or aerosols with certain properties. For brevity, computational particles are hereafter referred to as particles.

This overview of the code is divided into three parts: meteorological input data used by FLEXPART (see Appendix A1), computations related to the transport of particles (see Appendix A2), and computations affecting the properties of particles (see Appendix A3). We will only provide short summaries, since these three parts have all been documented in publications accompanying previous releases (e.g. Stohl et al., 2005; Pisso et al., 2019) and are described in more detail in the manual at https://flexpart.img.univie.ac.at/docs (last access: 11 September 2024). The main purpose here is to document how these computations are now organised within the source code, following its restructuring.

A1 Input data

FLEXPART advances particles based on interpolated meteorological fields, namely grid-scale three-dimensional fields of wind velocities, temperature, specific humidity, cloud liquid water, and ice content, as well as precipitation and various surface fields. In principle, any gridded dataset could be used. Data formats and coordinate systems used, as well as differences in the meteorological variables provided, would, however, make a generic input interface rather complex. The main FLEXPART code described here supports two input formats, namely data from the ECMWF and NCEP forecast systems (IFS and GFS, respectively). For ECMWF data, the flex_extract software package (Tipka et al., 2020) is provided to extract, process, and store the required fields for use as FLEXPART input, including support for ECMWF's reanalyses. Note that in the case of ECMWF-based input, data on all model layers are used, whereas NCEP-based input comprises only pressure-level fields. Other formats could be read in by writing the appropriate gridcheck and readwind subroutines and adding them to the windfields_mod.f90 module. A full list of necessary input fields can be found in Table B1. With the exception of the option of reading three precipitation fields per input time interval (see Sect. 5), nothing has changed compared to FLEXPART 10.4.

Other variants of FLEXPART were developed to accommodate fields produced by specific meteorological models., e.g. FLEXPART-WRF (Brioude et al., 2013), FLEXPART-COSMO (Henne et al., 2016), FLEXPART-AROME (Verreyken et al., 2019), FLEXPART-HIRLAM (Foreback, 2023), and FLEXPART–NorESM/CAM (Cassiani et al., 2016). There also exists a FLEXPART version for very high-resolution simulations (1 km), where turbulence parameterisations have been adapted to account for the fact that turbulence is partly already resolved by the meteorological input data (Katharopoulos et al., 2022).

A2 Particle transport

The transport of particles is central to a LPDM and straightforward. FLEXPART advances particles using motion vectors. At each time step, these motion vectors combine the grid-scale wind velocity ($\overline{\mathit{v}}$) from linearly interpolated meteorological input data, the parameterised turbulent wind velocity (v_t; see Appendix A2.2), and, for aerosols, the settling velocity (v_s; see Appendix A2.3). The resultant motion vector for that time step is $\mathit{v}=\overline{\mathit{v}}+{\mathit{v}}_{\mathrm{t}}+{\mathit{v}}_{\mathrm{s}}$. The new particle position is then given by

$$\begin{array}{}\text{(A1)}& \mathit{X}(t+\mathrm{\Delta }t)=\mathit{X}\left(t\right)+\mathit{v}(\mathit{X},t)\mathrm{\Delta }t,\end{array}$$

where Δt is the internal time step of the model. In addition, particles may be displaced vertically due to convection (see Appendix A2.1).

As shown in Fig. A1, the time manager module (timemanager_mod.f90) initiates, at each time step, the various processes. Particles first undergo convection (if required), which only affects their vertical position, and then the advance module (advance_mod.f90) is called, separately for each particle. Depending on the vertical position of the particle, different turbulence parameterisations are selected (see Appendix A2.2). After having derived the displacement vector and obtained a first guess for the new location of the particle, one Petterssen correction step (Petterssen, 1940) is executed for greater numerical accuracy to move the particle to its final new position. This is a second-order numerical scheme, which is accurate enough for the purpose of trajectory calculations, given sufficiently short time steps (Seibert, 1993). Particle positions are updated within the particle module (particle_mod.f90), which is directly or indirectly called from the advance module.

Figure A1 Simplified call tree for the transport of particles; green boxes represent modules containing subroutines (blue boxes) responsible for updating the location of particles through time.

A2.1 Convection

Vertical transport by moist, sub-grid-scale convection is parameterised using the scheme of Emanuel and Živković-Rothman (1999). In FLEXPART 11, all convection-related computations are located in the convection module (conv_mod.f90). The steps required are handled by convmix and include the computation of the matrix describing the convective redistribution of mass in a grid column (calcmatrix), making use of the convection scheme of Emanuel and Živković-Rothman (1999) (convect and tlift), followed by the actual convective displacement of particles in redist. The sorting algorithm used in convmix has been changed to use the FORTRAN Standard Library (https://fortran-lang.org/, last access: 11 September 2024), which codes the algorithm of Musser (1997). This new sorting algorithm is 4 orders of magnitudes faster; it is located in a separate module, namely sort_mod.f90.

A2.2 Turbulence

In FLEXPART 11, all turbulence-related computations have been organised in subroutines inside the turbulence module. This includes the turbulence computations inside of the atmospheric boundary layer (ABL), based on the Hanna and Chaughy 1982 scheme (Hanna, 1982), with added features from Ryall and Maryon (1998) and a skewed turbulence option by Cassiani et al. (2015); the turbulence above the ABL, a parameterisation based on pseudo-random number sampling and the length of the computational time step (Legras et al., 2003); and the mesoscale turbulence, parameterised using the method described by Maryon (1998), which uses the standard deviation of the wind velocities in the wind fields at the surrounding grid points. These turbulence schemes are called from the advance module and can now be completely turned off by an option in the COMMAND file by setting LTURBULENCE=0. Mesoscale turbulence is switched off by default and can be switched on in addition by setting LTURBULENCE_MESO=1. Below, we give a short description of the different turbulence schemes and options that are currently implemented in FLEXPART.

Turbulence in the atmospheric boundary layer

Turbulence in the ABL is parameterised, assuming a Markov chain based on the Langevin equation (Thomson, 1987). The Lagrangian timescales and turbulent velocity statistics for Gaussian turbulence are computed following either the scheme of Hanna (1982) or that of Ryall and Maryon (1998). Both distinguish neutral, stable, and unstable conditions. Turbulent motions can either be calculated using FLEXPART's fixed synchronisation time step, Δt_s, set by LSYNCTIME in the COMMAND file or use shorter time steps, depending on the Lagrangian timescale Δt_L. The former option is computationally very efficient but is inaccurate, as velocity auto-correlation cannot be taken into account if Δt_s>Δt_L. It is suitable only for long-range transport processes, where details of turbulence in the ABL are less important. Adaptive time steps are enabled by CTL>1 in the COMMAND; CTL determines by which factor the computational time step is kept shorter than Δt_L. For Gaussian turbulence, CTL values of at least 5 are recommended, and for the skewed turbulence scheme, even larger values are needed. With CTL>1, the time steps for the vertical component of turbulent motion can be further refined by setting IFINE>1, with the value of IFINE determining the factor by which the time step is further reduced.

FLEXPART offers the option to choose between Gaussian (Stohl and Thomson, 1999) and skewed turbulence (Cassiani et al., 2015). Gaussian turbulence is a suitable approximation for stable and neutral conditions. The error when using Gaussian turbulence under convective conditions is relatively small for sufficiently long transport distances, where particles become well mixed throughout the ABL. Only for short-range (i.e. ≲ 1 h) dispersion in the convective ABL, however, is it recommended to use the skewed turbulence scheme (Cassiani et al., 2015) by setting the switch CBL to 1 in the COMMAND file, as this scheme considerably increases computation time.

Turbulence in the free troposphere and stratosphere

Turbulent velocities above the ABL, both in the stratosphere and in the troposphere, are computed following Legras et al. (2003), using a constant vertical diffusivity (D_z=0.1 m² s⁻¹) to compute vertical turbulent velocities in the stratosphere and a constant horizontal diffusivity (D_h=50 m² s⁻¹) to compute horizontal turbulent velocities in the free troposphere. These default values can be modified in the input files from version 10.4 onwards. A linear transition layer of 1 km is used between the troposphere and stratosphere. The tropopause is defined as the first stable layer to fulfil the thermal tropopause criterion (i.e. the vertical temperature gradient is smaller than 0.002 K km⁻¹). Standard deviations of each of the velocity components ($i=\mathrm{1},\mathrm{\dots },\mathrm{3}$) are then obtained as

$$\begin{array}{}\text{(A2)}& {\mathit{\sigma }}_{v_i}=\sqrt{{\displaystyle \frac{\mathrm{2}{D}_i}{\mathrm{\Delta }t}}}.\end{array}$$

Horizontal diffusion is neglected in the stratosphere and vertical diffusion in the free troposphere. Already 20 years ago, an attempt was made to include a clear-air turbulence (CAT) parameterisation in FLEXPART, but due to the difficulties inherent to this problem (CAT can only be diagnosed probabilistically, and it is hard to establish Lagrangian timescales for it), it was not pursued further. Recently work done for the NAME model (Mirza et al., 2024) may show a future way forward – also for FLEXPART.

A2.3 Gravitational settling

The gravitational settling of particles is called at each time step. The settling velocity of each particle is computed following Näslund and Thaning (1991), with a temperature-dependent dynamic viscosity, according to Sutherland (1893), and then added to the vertical motion. In FLEXPART 11, a settling module, settling_mod.f90, has been introduced containing all the relevant procedures. While previous versions of FLEXPART only considered spherical particles, FLEXPART 11 is also able to calculate the settling velocity for non-spherical particles and takes into account different types of particle orientations (see Sect. 4). Gravitational settling is only calculated if a particle carries a single aerosol species with non-zero mass. This is necessary since aerosols of different size, shape, or density exhibit different settling velocities; however, a particle can only follow a single trajectory.

A3 Evolution of particle properties

In the previous subsection, we discussed how particles are advanced, using parameterisation schemes and meteorological data interpolated to the particle position. In addition, FLEXPART tracks certain properties of one or more so-called species for each individual particle. The main property considered is the “mass” of the particles. Upon release, particles are assigned a certain mass, usually representing a fraction of the total mass to be released as specified in file RELEASES. In the most simple case, this mass is physical mass; if applied to radioactivity, it can be understood to represent activity. In the case of backward simulations, it represents mass (or activity) mixing ratios (Seibert and Frank, 2004). Nevertheless, for simplicity, we often refer just to the mass (roman font) of a particle.

During the simulation, FLEXPART takes into account various processes that may alter a particle's mass (or what mass actually represents) over time. For instance, particles can undergo dry and wet deposition, by which their mass is reduced, as described in Appendix A3.1 and A3.2. Additionally, radioactive decay and chemical reaction with the hydroxyl radical are other important processes that lead to a reduction in particle mass over time (Appendix A3.3). It should be noted that particles are hypothetical entities used for discretisation in the Lagrangian framework. They should not be considered actual aerosol particles. By considering the evolving mass of the particles, FLEXPART can simulate the transport, dispersion, and fate of gases and aerosols in a realistic manner. It allows the model to account for the varying lifetimes of different substances based on their specific properties and the environmental conditions they encounter.

A3.1 Dry deposition

Gases and aerosols can both undergo dry deposition, with the result of losing mass to the Earth's surface. The (dry-) deposition velocity v_d is calculated for a reference height (h_ref) of 15 m above ground. The value of h_ref can be changed by the user, keeping in mind that it should not be too low, so that a statistically sufficient number of particles are considered, and not too high, as the calculation of dry deposition becomes inaccurate if 2 h_ref falls outside the constant-flux layer. Dry deposition is applied for all particles below 2 h_ref as

$$\begin{array}{}\text{(A3)}& \mathrm{\Delta }m=m\left(t\right)\left[\mathrm{1}-\exp \left({\displaystyle \frac{-v_{\mathrm{d}}\left(h_{\text{ref}}\right)\mathrm{\Delta }t}{\mathrm{2}{h}_{\text{ref}}}}\right)\right],\end{array}$$

where Δm is the amount of mass transferred from the particle to the dry-deposition field. The deposition velocity is a function of height above ground and defined as

$$\begin{array}{}\text{(A4)}& v_{\mathrm{d}}\left(z\right)=-{\displaystyle \frac{F_C}{C\left(z\right)}},\end{array}$$

where F_C is the deposition flux to the surface, and C(z) is the concentration of the species at height z. This deposition velocity can be prescribed as a constant value (in the associated SPECIES file) or can be computed within FLEXPART, using the resistance method as described in Wesely and Hicks (1977) for gases and Slinn (1982) for aerosols. For using the resistance method, certain chemical and/or physical properties of the species considered must be provided by the user in the respective SPECIES file. For a full description of the dry-deposition scheme, see Stohl et al. (2005).

In FLEXPART 11, all routines related to dry deposition are collected in drydepo_mod.f90. For aerosols, the module now also considers new formulations of gravitational settling velocity for spherical and non-spherical particles as described in Sect. 4.

A3.2 Wet deposition

The wet-deposition scheme in FLEXPART distinguishes between aerosols and gases and between in-cloud and below-cloud scavenging. In FLEXPART 11, cloud parameters and precipitation are interpolated to the particle position, and the scheme for below-cloud scavenging of aerosols has been replaced. Routines related to wet deposition are now collected in the wetdepo_mod.f90 module. The revised wet-deposition scheme is explained in more detail in Sect. 5.

The in-cloud scavenging is based on the approach of Hertel et al. (1995), with advancements introduced by Grythe et al. (2017). The in-cloud scavenging coefficient, λ_i (s⁻¹), is defined by (Hertel et al., 1995)

$$\begin{array}{}\text{(A5)}& {\mathit{\lambda }}_{\mathrm{i}}={\displaystyle \frac{S_{i}I}{H_i}},\end{array}$$

where S_i is the dimensionless scavenging ratio between the concentration of a substance in precipitation and the concentration in air, I is the precipitation rate (m s⁻¹), and H_i is the cloud depth at which the precipitation occurs (in m). For gases, S_i is defined as

$$\begin{array}{}\text{(A6)}& S_i=\left({\displaystyle \frac{\mathrm{1}-c_{\mathrm{l}}}{H_{\text{eff}}RT}}+c_{\mathrm{l}}\right).\end{array}$$

Here, H_eff is Henry's law constant, R is the perfect gas constant, and T is the temperature. c_l is the cloud liquid water content (expressed in m³ of water per m³ of cloud air), a dimensionless quantity.

For particles, S_i is defined as

$$\begin{array}{}\text{(A7)}& S_i={\displaystyle \frac{f_{\text{nuc}}}{c_{\mathrm{l}}}},\end{array}$$

where f_nuc is the fraction of aerosol that is activated. Grythe et al. (2017) introduced a new parameterisation for f_nuc, which takes into account the different efficiencies of aerosols in acting as cloud condensation nuclei or ice nuclei.

In the original Hertel et al. (1995) work, c_l was obtained from the purely empirical relationship as follows:

$$\begin{array}{}\text{(A8)}& c_{\mathrm{l}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}\cdot I^{\mathrm{0.36}}.\end{array}$$

This formulation is still used in FLEXPART when the cloud liquid water content c_l is not available. Aside from the new parameterisation of f_nuc, Grythe et al. (2017) also introduced an improved expression for λ_i that can be expressed in terms of S_i,

$$\begin{array}{}\text{(A9)}& S_i=f_{\text{nuc}}{\displaystyle \frac{H_{i}I{\mathit{\rho }}_{\text{water}}}{c_{\text{TWC}}I}}{r}_{\text{icl}},\end{array}$$

where c_TWC is the cloud total water content (kg m⁻²), and I ρ_water is the mass of water precipitating per unit time and area (kg s⁻¹ m⁻²). The dimensionless empirical constant r_icl is introduced to account for the replenishment rate of cloud water during precipitation and can be set in FLEXPART by incloud_ratio in par_mod.f90. This accounts for the replenishment of cloud water from condensing water vapour (Grythe et al., 2017). Noting that the cloud water washout ratio (R_w in s⁻¹) can be defined as

$$\begin{array}{}\text{(A10)}& R_{\mathrm{w}}={\displaystyle \frac{I{\mathit{\rho }}_{\text{water}}}{c_{\text{TWC}}}},\end{array}$$

introducing the in-cloud replenishment correction, r_icl, gives

$$\begin{array}{}\text{(A11)}& R_{\mathrm{w}}={\displaystyle \frac{I{\mathit{\rho }}_{\text{water}}}{c_{\text{TWC}}}}{r}_{\text{icl}},\end{array}$$

and results in a decreased washout ratio (i.e. increased washout timescale) by r_icl<1. We note that the current definition of r_icl (a dimensionless empirical parameter) is different from the original one used in Grythe et al. (2017), since they included the density of water in the empirical constant that was therefore dimensional, i.e. r_icl,Grythe = r_iclρ_water. This has now been changed in the FLEXPART 11 code for clarity and readability. Removing the density of water from the empirical constant means that the value of r_icl,Grythe=6.2 kg m⁻³, as reported in Grythe et al. (2017), now corresponds to r_icl=0.0062 (dimensionless). For completeness, as in Grythe et al. (2017), we rewrite the in-cloud scavenging coefficient in terms of the new expression for S_i directly as

$$\begin{array}{}\text{(A12)}& {\mathit{\lambda }}_{\mathrm{i}}=f_{\text{nuc}}{\displaystyle \frac{I{\mathit{\rho }}_{\text{water}}}{c_{\text{TWC}}}}{r}_{\text{icl}}.\end{array}$$

As a fallback in the case of lacking cloud water data, FLEXPART 11 uses the simple parameterisation for the total scavenging, which was used since the early versions of the model and which is common for regulatory nuclear applications (BASE, 2012), namely Λ=A I^B. Parameters A and B have to be provided in the SPECIES file as weta and wetb. Note that we use Λ here instead of λ to be in agreement with previous literature referring to bulk values.

In previous FLEXPART versions, the precipitation rate was augmented on the base of a sub-grid-scale parameterisation, considering that only a fraction (F<1) of a grid cell would actually experience precipitation (Stohl et al., 2005). It was based on horizontal-resolution (150 km × 150 km) data (Hertel et al., 1995). As it is unclear which values would be appropriate for various finer grids of current or future meteorological models, we considered it justified to remove this parameterisation in FLEXPART 11. The wet scavenging in convective clouds will need to receive further attention in the future, ideally coupling it to the convection scheme.

A3.3 Radioactive and chemical decay

FLEXPART is able to account for radioactive and/or chemical decay of particles by defining a half-life $T_{\mathrm{1}/\mathrm{2}}$ (parameter pdecay >0 in the corresponding SPECIES file; there is no decay if pdecay <0). The decay affects the mass on particles travelling through the atmosphere and deposited mass as follows:

$$\begin{array}{}\text{(A13)}& m(t+\mathrm{\Delta }t)=m\left(t\right)e^{\mathrm{-}\mathit{\alpha }\mathrm{\Delta }t},\end{array}$$

with α being the decay constant,

$$\begin{array}{}\text{(A14)}& \mathit{\alpha }={\displaystyle \frac{\ln \mathrm{2}}{T_{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

The treatment of radioactive and/or chemical decay remains unchanged compared to previous versions. Decay is computed alongside dry and wet deposition and can be found in the corresponding modules. Decayed mass is subtracted from the mass of each particle in the time-manager module (timemanager_mod.f90).

A3.4 Chemical reactions

With the introduction of the linear chemistry module into FLEXPART, simple linear reactions can generally be computed (e.g. OH and Cl reactions). Loss processes related to reactions with radicals are represented as a first-order linear approximation in FLEXPART 11 — identical to FLEXPART 10.4 (Pisso et al., 2019) for OH but now expanded to any linear reaction. The chemistry-related mass loss, m, is calculated as

$$\begin{array}{}\text{(A15)}& m(t+\mathrm{\Delta }t)=m\left(t\right)e^{\mathrm{-}\mathit{\kappa }\mathrm{\Delta }t},\end{array}$$

with Δt being the reaction time step given by lsynctime. The temperature-dependent reaction rate κ (s⁻¹) is then calculated as

$$\begin{array}{}\text{(A16)}& \mathit{\kappa }=C{T}^N{e}^{\frac{-D}{T}}{c}_{\text{reagent}},\end{array}$$

where C, N, and D are the species-specific constants defined in the SPECIES file (parameters pcconst, pnconst, and pdconst, respectively; turned off with negative values). T is the absolute temperature, and c_reagent is the hourly concentration of the reagent (Atkinson, 1997).

The OH radical is the most important oxidant in the troposphere, and although atmospheric chemistry can be highly non-linear, a first-order linear loss approximation using prescribed OH fields is often still adequate, e.g. to simulate methane in the atmosphere (for FLEXPART 11 there are monthly averaged 3°×5° OH fields with 17 vertical layers, following the GEOS-Chem model by Bey et al. (2001) and read in from NetCDF files). Hourly OH variations are accounted for by modifying the monthly fields with the hourly photolysis rate of ozone j, based on a simple parameterisation for cloud-free conditions and depending on the solar zenith angle,

$$\begin{array}{}\text{(A17)}& \text{OH}={\displaystyle \frac{j}{j^{*}}}{\text{OH}}^{*},\end{array}$$

where j^* and OH^* are the monthly mean photolysis rate and OH concentration taken from the OH fields, respectively.

In order to be able to use chemistry fields with higher spatial and temporal resolution (as in, e.g., Fang et al., 2016, for OH), the user has to implement the reading and utilising of such data in the chemistry-related FLEXPART subroutines (readreagents, getchemfield, readchemfield, getchemhourly, and chemreaction located in the chemistry_mod.f90 module).

Appendix B: SPECIES files

Species in FLEXPART are either gases or particles (the case of air tracer can be considered an inert gas); each category requires different parameters to be specified in the respective SPECIES file. Values of parameters which are not required are ignored. The parameter values contained in the template SPECIES files bundled with FLEXPART are listed in Table B3 for gaseous substances and in Table B4 for particulate substances.

B1 Gaseous species

The parameters required for the dry deposition of gases are the inverse of the diffusivity relative to that of H₂O (D; PRELDIFF), the reactivity relative to that of O₃ (f₀; PF0, originally taken from Wesely (1989)), and Henry's constant (PHENRY; describing solubility). The relative diffusivity of gases (CO, SO₂, CH₄, C₂H₆, PCB−28, γ−HCH (lindane), and N₂O) could be calculated approximately using the square root of the ratio of molar weights between water and the respective species. The source for the relative reactivity f₀ is Table 2 in Wesely (1989), with values for CO and CH₄ taken from Clifton et al. (2022).

A collection of values of Henry's constant for many atmospheric compounds was recently compiled by Sander (2023); the list of values can also be found at https://www.henrys-law.org/henry/ (last access: 11 September 2024). For most of the species, they are close to or identical to the previously used values from Wesely (1989). More significant differences were found for HNO₃, SO₂, NH₃, and HNO₂. Now, values from Sander (2023) haven been used, except for the persistent organic pollutants, where well-established values exist from other sources, namely for PCB28 from Mackay et al. (2006) (as in FLEXPART 10.4) and for γ-HCH from Sahsuvar et al. (2003) (modified).

Concerning the parameters relevant for the wet deposition of gases, in-cloud scavenging depends on Henry's constant and the diffusivity, which were already discussed above. The values of weta for the coefficient of the simplified, fallback scavenging parameterisation were originally taken from Asman (1995) for gases. Now, some values have been modified, as described in Table B3.

Radioactive noble gases are inert and do not undergo relevant wet or dry deposition. Depending on their half-life and the timescales under consideration, their decay may be considered. Otherwise they are to be treated like air tracers. Radioactive decay (for gaseous and particulate species) may either be simulated directly in the FLEXPART simulation or be considered during post-processing. The latter approach may be useful if several nuclides with different half-lives but otherwise identical properties are simulated, so that not only CPU time but also memory and output file size can be reduced. In some cases, activities are desired for a given reference time – then it would also not make sense to include the decay in the simulation. For applying a decay correction to the mean concentration or deposition during an output time interval for species simulated without decay, see Seibert et al. (2013). Half-life data of different nuclides can be obtained from the International Atomic Energy Agency's LiveChart of Nuclides at https://www-nds.iaea.org/relnsd/vcharthtml/VChartHTML.html (last access: 11 September 2024), or the U.S. National Nuclear Data Center at Brookhaven National Laboratory at https://www.nndc.bnl.gov/ (last access: 11 September 2024); half-life values given in days need to be multiplied by 86 400 to obtain the value in seconds as needed for FLEXPART.

B2 Particulate species

The most important parameter for both the dry and wet deposition of particles is their diameter, or, if maxndia is set to >1 in par_mod.f90, their size distribution, which in FLEXPART is assumed to be log-normal. Thus, a mean geometric diameter and a logarithmic standard deviation have to be specified.

For the mean geometric particle diameter (PDIA) for accumulation-mode aerosol, such as ammonium nitrate or sulfate, as well as radionuclides bound to such aerosol particles, a typical value of 0.4 µm is proposed, as in previous versions. Note that this is a quantity that is variable according to the environmental conditions. We have also introduced a template for fine and coarse dust, with typical diameters of 0.2 µm and, respectively, 10 µm.

Radioactive iodine as released from nuclear reactors in the case of accidents mainly consists of gaseous elemental iodine (I₂). This gas tends to deposit on accumulation-mode aerosol particles. Some radioiodine may also be present as methyl iodide (CH₃I) (Nair et al., 2000). A new SPECIES file was created for methyl iodide, with wet-deposition parameters of elemental and methyl iodide taken from Asman (1995) and Päsler-Sauer (2000). Concerning dry deposition, Henry's constant and reactivity data could not be found in the literature, so a mean value for the deposition velocity (v_dep) of $\mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ was introduced for methyl iodide (Müller and Pröhl, 1993). Particle-bound iodine behaves like other particle-bound nuclides (see below).

Concerning the wet deposition of particle-bound caesium-137, Van Leuven et al. (2023) have tried to adjust the FLEXPART deposition parameters based on inverse modelling; however, as the wet-deposition scheme has been changed from version 10.4 (used by these authors) to version 11, and the inverse method is also subject to significant uncertainties, we generally refrain from recommending their values.

Bolton (1980)(Berkowicz and Prahm, 1982)Berkowicz and Prahm (1982)

Table B1 List of input variables.

Table B2 Set of equations to calculate the drag coefficient C_d of aerosol particles for three types of orientation based on the shape correction scheme of Bagheri and Bonadonna (2016). L, I, and S are the longest, intermediate, and smallest axes of a particle, respectively. ρ_p and ρ_f are the density of the particle and the density of the fluid, respectively. k_S is Stokes' and k_N Newton's drag correction, F_S is Stokes' form factor, F_N is Newton's form factor, d_eq is the diameter of a sphere of equivalent volume, and α and β are sigmoidal functions of the particle-to-fluid density ratio of ${\mathit{\rho }}^{\prime }\equiv \frac{{\mathit{\rho }}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{f}}}$. If non-spherical particles are assumed to be ellipsoids or fibres with $L/I>\mathrm{20}$, the term $\left(\frac{d_{\mathrm{eq}}^{\mathrm{3}}}{LIS}\right)$ can be neglected for the calculation of F_S and F_N (i.e. F_S=fe^1.3; F_N=f²e). f and e are the flatness (PSA/PIA) and elongation (PIA/PLA) of the particles, respectively. For more details, see Bagheri and Bonadonna (2016).

Table B3 List of gaseous FLEXPART species, including their parameter values and references, if available. The asterisk (^*) denotes values modified in FLEXPART 11, and the double asterisk (${}^{**}$) denotes newly added species. $T_{\mathrm{1}/\mathrm{2}}$ refers to the half-life in seconds; weta and wetb relate to dimensionless below-cloud scavenging coefficients; reldiff, f0, and Henry's constant (in M atm⁻¹) relate to parameters for dry deposition; M relates to the molecular weight in mol; and cconst, dconst, and nconst relate to C (in ${\mathrm{cm}}^{\mathrm{3}}\mathrm{molec}{.}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$), D (in K), and dimensionless N, respectively, used to calculate reaction rates (see Eq. A16).

n/a: not applicable. If no other values are known, $\text{weta}=\mathrm{1}E-\mathrm{5}$ and wetb=0.8 are recommended for aerosol particles (BASE, 2012). ^a Asman (1995); ^b Wesely (1989); ^c Clifton et al. (2022); ^d Atkinson (1997); ^e Mackay et al. (2006); ^f Sahsuvar et al. (2003); ^g Grythe et al. (2017); ^h Tunved et al. (2013); ⁱ Groot Zwaaftink et al. (2022); ^j Päsler-Sauer (2000); ^k https://www.henrys-law.org/henry/ (last access: 11 September 2024).

Table B4 List of particulate FLEXPART species, including their parameter values and references, if available. The single asterisk (^*) denotes values modified in FLEXPART 11, and the double asterisk (${}^{**}$) denotes newly added species. $T_{\mathrm{1}/\mathrm{2}}$ refers to the half-life in second; crain_aero and csnow_aero are below-cloud scavenging efficiencies; ccn_aero and in_aero are in-cloud scavenging parameters; density gives the density of an aerosol particle (in kg m⁻³); dia is the diameter of an aerosol particle in metres; dsigma is the geometric standard deviation. A template for particulate radionuclides and an irregularly shaped particle is provided in the repository.

n/a: not applicable. If no other values are known, $\text{weta}=\mathrm{1}E-\mathrm{5}$ and wetb=0.8 are recommended for aerosol particles (BASE, 2012). ^a Asman (1995); ^b Wesely (1989); ^c Clifton et al. (2022); ^d Atkinson (1997); ^e Mackay et al. (2006); ^f Sahsuvar et al. (2003); ^g Grythe et al. (2017); ^h Tunved et al. (2013); ⁱ Groot Zwaaftink et al. (2022); ^j Päsler-Sauer (2000); ^k https://www.henrys-law.org/henry/ (last access: 11 September 2024).