3.1 The η-coordinate system

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

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:

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 1Mean 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.

Download

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

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 2Averaged 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).

Download

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.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

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):

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:

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.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 5Global 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.

Download

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 1Statistical 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}}.$

Download Print Version | Download XLSX

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.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 7Demonstration 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.

Download

Figure 8The 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.

Download

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

Download Print Version | Download XLSX

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.

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 9Methane 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.