FALL3D-8.0: a computational model for atmospheric transport and deposition of particles, aerosols and radionuclides – Part 1: Model physics and numerics

This paper presents FALL3D-8.0, the last version release of an open-source code with 15+ years of track record and a growing number of users in the volcanological and atmospheric communities. The code has been redesigned and rewritten from scratch in the framework of the EU Centre of Excellence for Exascale in Solid Earth (ChEESE) in order to overcome legacy issues and allow for successive optimisations in the preparation of the code towards extreme-scale computing. However, this baseline version already contains substantial improvements in terms of model physics, solving algorithms, and code accuracy and performance. The code, originally conceived for atmospheric dispersal and deposition of tephra particles, has been extended to model other types of particles, aerosols and radionuclides. The solving strategy has also been changed, replacing the former central-difference scheme for a high-resolution central-upwind scheme derived from finite volumes, which minimises numerical diffusion even in the presence of sharp concentration gradients and discontinuities. The parallelisation strategy, input/output (I/O), model pre-process workflows and memory management have also been reconsidered, leading to substantial improvements on code scalability, efficiency and overall capability to handle much larger problems. All these new features and improvements have implications on operational model performance and allow, among others, adding data assimilation and ensemble forecast in future releases. This paper details the FALL3D-8.0 model physics and the numerical implementation of the code.

1 Introduction

FALL3D is an open-source offline Eulerian model for atmospheric passive transport and deposition based on the so-called advection–diffusion–sedimentation (ADS) equation. The model, originally developed for inert volcanic particles (tephra), has a track record of 50+ publications on different model applications and code validation, as well as an ever-growing community of users worldwide, including academic, private and research institutions and several institutions tasked with the operational forecast of volcanic ash clouds. The first versions of FALL3D (v1.x series) appeared back in 2003 (Costa and Macedonio, 2004); at that time, the code was serial and written in Fortran 77. Code improvements at different levels have been continuously incorporated since then, including relevant milestones leading to code version upgrades, e.g. the coupling with 1-D buoyant plume theory (BPT) models to define eruption column sources (v2.x, 2004), the introduction of the Lax–Wendroff (LW) central-difference scheme for solving the ADS equation (v3.x, 2005) and other algorithmic improvements (Costa et al., 2006), full code rewriting in Fortran 90 and distributed memory parallelisation by means of Message Passing Interface (MPI) (v5.x, 2007), first implementations of operational workflows to forecast ash cloud dispersal and fallout (Folch et al., 2008, 2009), and several other improvements (e.g. de la Cruz et al., 2016) until the v7.3.4 release in 2018.

Along these 15+ years, FALL3D has been used in multiple applications (e.g. Folch, 2012) including, among others, assessment of hazards from tephra fallout at different volcanoes (e.g. Scaini et al., 2012; Selva et al., 2014; Sandri et al., 2016), impacts of ash cloud dispersal on civil aviation (e.g. Sulpizio et al., 2012; Bonasia et al., 2013; Biass et al., 2014; Scaini et al., 2014), obtaining relevant eruption source parameters (e.g. Folch et al., 2012; Parra et al., 2016; Poret et al., 2017), impacts of past super-eruptions on climate, environment and humans (e.g. Costa et al., 2012, 2014; Martí et al., 2016), operational forecast of ash clouds and tephra fallout (e.g. Bear-Crozier et al., 2012; Collini et al., 2013; Poulidis et al., 2019), modelling of ash resuspension events (e.g. Folch et al., 2014; Mingari et al., 2017) or model validation (e.g. Scollo et al., 2010; Corradini et al., 2011; Osores et al., 2013). However, as in other long-lived community codes, code legacy limitations have appeared with time on, e.g. lack of code performance and poor scalability on hundreds/thousands of cores, constraints on portability and adaption to emerging hardware architectures and difficulties for code refactoring that is needed in order to widen the spectrum of model applications. With time, properly addressing these aspects required of substantial code refactoring or even code rewriting, a substantially time-consuming task in terms of human effort. This has recently been possible in the framework of the European Centre of Excellence for Exascale in Solid Earth (ChEESE), which includes FALL3D as one of its flagship codes.

This paper describes FALL3D-8.0, a new model version upgrade in which the code has been completely written from scratch, mostly in Fortran 90 but mixed with some Fortran 2003 functionalities. In addition to dramatic improvements on different levels (extended model physics and applications, numerical algorithmic and code performance), v8.0 also provides with a baseline that will allow the incorporation of developments and optimisations. However, the heterogeneity of model users has been considered when rewriting the code, which can still run on platforms spanning from a laptop to a large supercomputer.

This paper starts first outlining what is new in FALL3D-8.0 (Sect. 2) with respect to the previous code release (v7.3.4) and then describes the physical model and related governing equations and parameterisations (Sect. 3). The next section focuses on the numerical implementation (Sect. 4), including coordinate mappings and scaling, spatial discretisation and a new solving strategy based on the Kurganov–Tadmor scheme (Kurganov and Tadmor, 2000) that can be combined either with a fourth-order Runge–Kutta or with a first-order Euler to integrate explicitly in time. These two solver options allow users to choose, respectively, between better solver accuracy or higher computational efficiency. In any case, it will be shown how the Kurganov–Tadmor finite-volume-like formulation is much less diffusive than the previous scheme implemented in v7.x, an important feature when one aims at modelling substances encompassing sharp gradients of concentration. Sections 5 and 6 focus, respectively, on the (pre-process) model workflow and on the new code parallelisation strategy and related memory optimisations, including a comparison of model performance and scalability with respect to v7.x. This paper shows only one illustrative example of FALL3D-8.0 model results for ash dispersal from the 2011 Cordón Caulle eruption. A companion paper (Prata et al., 2020) contains the second part, including a detailed FALL3D-8.0 model validation for several simulations that are part of the new benchmark suite of the code. Finally, Sect. 7 wraps the conclusions of the paper and outlines which will be the next steps for further model optimisation.

2 New features in v8.0

FALL3D-8.0 introduces substantial improvements at different levels. From the point of view of model physics, the following points apply:

• The code has been generalised to deal with species different from tephra, including other kinds of particles, aerosols and radionuclides. Different types of species have been defined and can be simulated using independent sets of bins.

• Weibull and bi-Weibull particle total grain size distributions (TGSDs) can now be generated in addition to lognormal distributions. On the other hand, FALL3D-8.0 can estimate TGSDs of tephra particles directly from magma viscosity (magma composition) and eruption column height following the fit proposed by Costa et al. (2016a).

• An “effective bin” flag has been introduced. For a given species, only those bins having this flag “on” are actually simulated, whereas bins tagged as “off” are frozen in terms of transport and used only for source term characterisation. This option has been added because volcanic plume source parameterisations depend on the whole granulometric spectrum of particles but, at the same time, model users are often interested only on a subset of the particle spectra (e.g. fine volcanic ash for long-range tephra dispersal simulations).

• For the species “tephra”, several classes of particle aggregates can now be defined in certain aggregation options, differently than the single-class aggregation option available in v7.x.

• Model parameterisations for physics have been revised to include more recent studies. Meteorological drivers have also been updated to add recent datasets (e.g. ERA5) and to remove deprecated options.

• Periodic boundary conditions are now possible, permitting simulations on domains spanning from local to semi-global (pole singularities still remain).

• For some species, the initial model condition can be furnished from satellite retrievals. This “data insertion” option is a preliminary step towards a full model data assimilation using ensembles.

From the point of view of numerics and code performance, the following points apply:

• The solving strategy has been changed to a high-resolution central-upwind scheme (Kurganov and Tadmor, 2000), optionally combined with a fourth-order Runge–Kutta explicit scheme. This replaces the former LW central-difference scheme, which was known to be over-diffusive in the case of sharp concentration gradients or discontinuities.

• New coordinate mappings have been added. These include a new vertical σ-coordinate system with linear decay that smooths numerical oscillations over complex terrains (Schar et al., 2002).

• A new parallelisation strategy exists based on a full 3-D domain decomposition. The former trivial parallelisation on particle bins in v7.x has been removed because, in the case of interaction among bins, it yielded to unnecessary communication penalties, resulting in poor code scalability.

• A much more efficient memory management exists to exploit contiguous cache memory positions along each spatial dimension. In some model configurations, this can imply a substantial gain on computing time.

• Parallel model I/O using NetCDF and parallel model pre-process. In addition, all the pre-process auxiliary programmes have been embedded within the main code (a single multipurpose executable exists in v8.0). The code can now be run to perform different tasks individually or sequentially through a single parallel workflow. As a result, all the pre-process, modelling and post-process workflows can now run as a single execution concatenating all tasks and without needing to write/read intermediate files to/from disk. In large problems, this saves substantial disk space and I/O time because the required model input data (e.g. interpolated meteorological fields) are already stored in each processor memory when running the FALL3D model task.

• A hierarchy of MPI communicators has been defined. This is actually not active yet in v8.0 but gives flexibility to extend some code functionalities in a near future with little refactoring effort. For example, plans for future versions include ensemble modelling using the Parallel Data Assimilation Framework (PDAF) or model nesting. These will require of teams of processors associated with different ensemble members or model grids, respectively.

3 Physical model

3.1 Model governing equations

In continuum mechanics, the general form describing the passive transport of a substance mixed within a fluid (air) in a domain Ω derives from mass conservation which, in conservative flux form and using the Eulerian specification, reads (see Tables 1 and 2 for list of symbols)

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial c}{\partial t}}+\mathrm{\nabla }\cdot \mathit{F}+\mathrm{\nabla }\cdot \mathit{G}+\mathrm{\nabla }\cdot \mathit{H}=S-I\text{in}\mathrm{\Omega },\end{array}$$

where F=c u is the advective flux, G=c u_s is the sedimentation flux, $\mathit{H}=-\mathbb{K}\mathrm{\nabla }c$ is the diffusive flux, and S and I are the source and sink terms, respectively. In the expressions above, t denotes time, c is concentration (in kg m⁻³), u is the fluid velocity vector (i.e. wind velocity), u_s is the terminal settling velocity of the substance, and 𝕂 is the diffusion tensor (in m² s⁻¹). Note that the definition of the diffusive flux H explicitly assumes Fick's first law of diffusion.

Boundary conditions are imposed at the Dirichlet Γ_D (inflow), Neumann Γ_N (outflow) and Robin Γ_R (ground) parts of the boundary of the computational domain Γ (with $\mathrm{\Gamma }={\mathrm{\Gamma }}_{\text{D}}\cup {\mathrm{\Gamma }}_{\text{N}}\cup {\mathrm{\Gamma }}_{\text{R}}$ and ${\mathrm{\Gamma }}_{\text{D}}\cap {\mathrm{\Gamma }}_{\text{N}}\cap {\mathrm{\Gamma }}_{\text{R}}=\mathrm{0}$) as

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{ll}c=\stackrel{\mathrm{‾}}{c}& \text{at}{\mathrm{\Gamma }}_{\text{D}};\\ \mathit{n}\cdot \mathit{H}=\mathrm{0}& \text{at}{\mathrm{\Gamma }}_{\text{N}};\\ \mathit{n}\cdot (\mathit{H}+\mathit{G})=\mathit{n}\cdot \mathit{D}& \text{at}{\mathrm{\Gamma }}_{\text{R}},\end{array}\right.\end{array}$$

where $\stackrel{\mathrm{‾}}{c}$ is the concentration prescribed at inflow (typically $\stackrel{\mathrm{‾}}{c}=\mathrm{0}$), n is the outwards unit normal vector, and D=cu_d is the ground deposition flux (u_d is the ground deposition velocity). Note that when the deposition flux D coincides with the sedimentation flux G (i.e. when u_s=u_d), the boundary condition at ground reduces to the standard free-flow condition imposed at Γ_N.

Table 1 List of Latin symbols. Bold and doubled-lined fonts indicate vector and tensor quantities, respectively.

Equation (1) is the so-called ADS equation and, in FALL3D-8.0, has been extended to handle other substances different from tephra. Note that the passive transport equation (Eq. 1) neglects inertial terms and, consequently, assumes a low (≤ 1) particle Stokes number. In a general sense, substances in FALL3D-8.0 are grouped in three broad categories: particles, aerosols and radionuclides. The category “particles” includes any inert substance characterised by a sedimentation velocity. In the context of FALL3D-8.0, the category “aerosol” refers to substances potentially non-inert (i.e. with chemical or phase changes mechanisms) and having a negligible sedimentation velocity. Finally, the category “radionuclides” refers to isotope particles subject to radioactive decay. Each of these categories admits, in turn, different subcategories or “species”, defined internally as structures of data that inherit the parent category properties (see Table 3). For example, particles can be subdivided into tephra or mineral dust; aerosol species can include H₂O, SO₂, etc.; and radionuclides can include several isotope species. Finally, each subcategory of species is tagged with an attribute name that is used for descriptive purposes only.

Depending on the species under consideration, the mass source S and sink I terms in Eq. (1) can be decomposed as

$$\begin{array}{ll}{\displaystyle }S=& {\displaystyle }{S}^{\text{e}}+S^{\text{a}}+S^r+S^{\text{c}}\\ \text{(3)}& {\displaystyle }I=& {\displaystyle }{I}^{\text{w}}+I^{\text{a}}+I^r+I^{\text{c}},\end{array}$$

where the superscripts denote emission source terms (S^e; see Sect. 3.2.3), wet deposition sinks (I^w; see Sect. 3.2.4), particle aggregation source and sinks (S^a and I^a, respectively; see Sect. 3.2.6), radioactive decay source and sinks (S^r and I^r, respectively; see Sect. 3.2.7), and chemical reactions' source and sinks (S^c and I^c). Note that FALL3D-8.0 does not account for aerosol chemistry yet. However, the code has been designed to allow incorporating this functionality in future versions in a straightforward manner.

Table 2 List of Greek symbols.

Table 3 Types of categories and related subcategories of species in FALL3D-8.0.

¹ For any species in the category particles or radionuclides, users can specify the number of effective bins from a grain size distribution. ² For tephra, the Φ number is defined as $d={\mathrm{2}}^{-\mathrm{\Phi }}$, where d is the particle diameter in millimetres. ³ SO₂ chemistry not included yet in v8.0.

Table 4 Parameterisations available in FALL3D-8.0 depending on each category and related subcategory (species). Note that, except for radioactive decay, all physical phenomena were already included in v7.x. However, several parameterisations have been updated to account for more recent developments.

¹ Only for volcanic aerosols. ² Applies only to particles/aerosols smaller than 100 µm. ³ Cut-off at 100 and 1 µm assumed for below- and in-cloud scavenging, respectively. ⁴ Not included yet in v8.0.

When the ADS equation (Eq. 1) is discretised, species in the mixture are binned in n_b discrete “classes” or bins, so that the total concentration c at any point of the domain decomposes as the sum of bin concentrations, i.e. $c={\sum }_{i=\mathrm{1}}^{n_{\text{b}}}{c}_i$. Substitution of bin discretisation in Eq. (1) yields to n_b equations (one per discrete bin), each formally identical to Eq. (1):

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial c_i}{\partial t}}+\mathrm{\nabla }\cdot {\mathit{F}}_i+\mathrm{\nabla }\cdot {\mathit{G}}_i+\mathrm{\nabla }\cdot {\mathit{H}}_i=S_i-I_{i}i=\mathrm{1}:n_{\text{b}}.\end{array}$$

Note that the n_b equations for bins can be coupled by means of the source and sink terms, which define the eventual transfer of mass among different bins, e.g. in the case of particle aggregation/disegregation, chemical reactions and formation of child radionuclides.

3.2 Model parameterisations

Parameterisations in FALL3D have been revised and updated, removing deprecated options and adding new options available from more recent studies. Table 4 shows the parameterisations implemented in FALL3D-8.0 that are described in the following subsections.

3.2.1 Diffusion tensor

The atmospheric flow is characterised by large horizontal-to-vertical aspect ratios of wind velocities and length scales, as well as by an anisotropic momentum diffusion with the horizontal diffusion coefficient being typically much larger than the vertical one (e.g. Schaefer-Rolffs and Becker, 2013). For this reason, model diffusion that accounts for subgrid-scale atmospheric eddies is typically assumed anisotropic, with two distinct eddy diffusion coefficients along the horizontal and vertical dimensions:

$$\begin{array}{}\text{(5)}& \mathbb{K}=\left(\begin{array}{ccc}{k}_{\text{h}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\text{h}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\text{v}}\end{array}\right).\end{array}$$

In the case of FALL3D-8.0, the horizontal coefficient k_h can be either assumed constant or parameterised as in Byun and Schere (2006):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{k_{\text{h}}}}={\displaystyle \frac{\mathrm{1}}{k_{\text{ht}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\text{hn}}}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{k}_{\text{ht}}& =\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}{\mathrm{\Delta }}_{\text{g}}^{\mathrm{2}}\sqrt{S_{\mathrm{\Gamma }}^{\mathrm{2}}+S_{\mathrm{\Lambda }}^{\mathrm{2}}}\\ & ={\mathit{\alpha }}^{\mathrm{2}}{\mathrm{\Delta }}_{\text{g}}^{\mathrm{2}}\sqrt{{\left({\displaystyle \frac{\partial u}{\partial x}}-{\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}}\end{array}\end{array}$$

$$\begin{array}{}\text{(8)}& k_{\text{hn}}=k_{\text{ref}}\left({\displaystyle \frac{{\mathrm{\Delta }}_{\text{ref}}}{{\mathrm{\Delta }}_{\text{g}}}}\right).\end{array}$$

In the equations above, Δ_g is a characteristic grid cell measure (e.g. the equivalent area length), α≅0.28 denotes the Smagorinsky constant, S_Γ and S_Λ are the stretching and shear strength (i.e. the two components of the bi-dimensional wind deformation), and k_ref is a reference horizontal diffusion for a reference grid cell size Δ_ref (FALL3D-8.0 considers k_ref=8000 m² s⁻¹ for Δ_ref=4 km). Equation (6) was proposed by Byun and Schere (2006) to overcome the dependency of horizontal diffusion on grid resolution and combines a Smagorinsky subgrid-scale (SGS) model giving the diffusion by transport (k_ht) with a formula that counteracts numerical over-diffusion in coarse discretisations (k_hn) so that the smaller one between k_ht and k_hn dominates. In this way, the effect of the transportive diffusion is minimised for coarse grids, whereas for fine discretisations the numerical diffusion term is reduced.

The vertical diffusion coefficient k_v can also be either assumed constant or parameterised according to the similarity theory and distinguishing among surface layer, atmospheric boundary layer (ABL) and free atmosphere (e.g. Neale et al., 2010):

$$\begin{array}{}\text{(9)}& k_v=\left\{\begin{array}{ll}\frac{\mathit{\kappa }z{u}_{*}}{{\mathit{\varphi }}_{\text{h}}}& \text{for}z\ll h_{\text{p}}\\ \frac{\mathit{\kappa }z{u}_{*}}{{\mathit{\varphi }}_{\text{h}}}{\left(\mathrm{1}-\frac{z}{h}\right)}^{\mathrm{2}}& \text{for}z<h_{\text{p}}\\ l_{\text{c}}^{\mathrm{2}}\left|\frac{\partial u_{\text{h}}}{\partial z}\right|F_{\text{c}}\left(\mathit{Ri}\right)& \text{for}z>h_{\text{p}},\end{array}\right.\end{array}$$

where κ is the von Karman constant (κ=0.4), z is the distance from the ground, u_* is wind friction velocity, ϕ_h is the atmospheric stability function for temperature, h_p is the ABL height, l_c is a characteristic length scale, u_h is the horizontal wind velocity modulus, and F_c is a stability function which depends on the Richardson number Ri. For l_c and F_c, FALL3D-8.0 adopts the relationships used by the Community Atmosphere Model (CAM) 4.0 model (Neale et al., 2010):

$$\begin{array}{}\text{(10)}& l_{\text{c}}={\left({\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }z}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\text{c}}}}\right)}^{-\mathrm{1}}\end{array}$$

$$\begin{array}{}\text{(11)}& F_{\text{c}}\left(\mathit{Ri}\right)=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{1}+\mathrm{10}\mathit{Ri}(\mathrm{1}+\mathrm{8}\mathit{Ri})}& \text{stable}(\mathit{Ri}>\mathrm{0})\\ \sqrt{\mathrm{1}-\mathrm{18}\mathit{Ri}}& \text{unstable}(\mathit{Ri}<\mathrm{0}),\end{array}\right.\end{array}$$

where λ_c is the so-called asymptotic length scale (λ_c≈30 m). The atmospheric stability function for temperature ϕ_h is calculated as

$$\begin{array}{}\text{(12)}& {\mathit{\varphi }}_{\text{h}}=\left\{\begin{array}{ll}{\mathit{\beta }}_{\text{h}}+\frac{z}{L}& z/L>\mathrm{1}\text{stable atmosphere}\\ \mathrm{1}+{\mathit{\beta }}_{\text{h}}\frac{z}{L}& \mathrm{0}\le z/L\le \mathrm{1}\text{nearly neutral}\\ {\displaystyle }{\left(\mathrm{1}-{\mathit{\gamma }}_{\text{h}}{\displaystyle \frac{z}{L}}\right)}^{-\mathrm{1}/\mathrm{2}}& z/L<\mathrm{0}\text{unstable atmosphere},\end{array}\right.\end{array}$$

where β_h=5, γ_h=15, and L is the Monin–Obukhov length, defined as

$$\begin{array}{}\text{(13)}& L={\displaystyle \frac{u_{*}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\text{v}}}{\mathit{\kappa }g{\mathit{\theta }}_{*}}},\end{array}$$

with g denoting gravity, ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\text{v}}$ the mean potential virtual temperature and θ_* the potential temperature scale. The parameters in Eq. (13) (i.e. L or $u_{*}/{\mathit{\theta }}_{*}$) are ideally furnished by the driving meteorological model. If not, and alternatively, FALL3D-8.0 estimates the friction velocity u_* and the potential temperature scale θ_* from the potential virtual temperature θ_v and the Richardson bulk number Ri_b as (Louis, 1979; Jacobson, 1999)

$$\begin{array}{}\text{(14)}& {\mathit{Ri}}_{\text{b}}={\displaystyle \frac{g\left[{\mathit{\theta }}_{\text{v}}\left(z_{\text{r}}\right)-{\mathit{\theta }}_{\text{v}}\left(z_{\text{o}}\right)\right]\left(z_{\text{r}}-z_{\text{o}}\right)}{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\text{v}}{u}_{\text{h}}(z_{\text{r}}{)}^{\mathrm{2}}}},\end{array}$$

where z_r and z_o denote the reference and the ground roughness heights, and ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\text{v}}$ denotes the average between the two vertical levels. Given Ri_b, one can estimate u_* and θ_* as

$$\begin{array}{}\text{(15)}& u_{*}\approx {\displaystyle \frac{\mathit{\kappa }{u}_{\text{h}}\left(z_{\text{r}}\right)}{\ln (z_{\text{r}}/z_{\text{o}})}}\sqrt{G_m\left({\mathit{Ri}}_{\text{b}}\right)}\end{array}$$

$$\begin{array}{}\text{(16)}& {\mathit{\theta }}_{*}\approx {\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}{u}_{\text{h}}\left[{\mathit{\theta }}_{\text{v}}\left(z_{\text{r}}\right)-{\mathit{\theta }}_{\text{v}}\left(z_{\text{o}}\right)\right]}{u_{*}{\mathit{Pr}}_{\text{t}}{\ln }^{\mathrm{2}}(z_{\text{r}}/z_{\text{o}})}}\sqrt{G_{\text{h}}\left({\mathit{Ri}}_{\text{b}}\right)},\end{array}$$

where Pr_t is the turbulent Prandtl number (Pr_t≈1) and the stability functions G_m and G_h are given by (Louis, 1979; Jacobson, 1999)

$$\begin{array}{}\text{(17)}& G_m=\left\{\begin{array}{ll}\mathrm{1}-\frac{\mathrm{9.4}{\mathit{Ri}}_{\text{b}}}{\mathrm{1}+\mathrm{70}{\mathit{\kappa }}^{\mathrm{2}}(\left|{\mathit{Ri}}_{\text{b}}\right|z_{\text{r}}/z_{\text{o}}{)}^{\mathrm{0.5}}/{\ln }^{\mathrm{2}}(z_{\text{r}}/z_{\text{o}})}& {\mathit{Ri}}_{\text{b}}\le \mathrm{0}\\ \frac{\mathrm{1}}{(\mathrm{1}+\mathrm{4.7}{\mathit{Ri}}_{\text{b}}{)}^{\mathrm{2}}}& {\mathit{Ri}}_{\text{b}}>\mathrm{0}\end{array}\right.\end{array}$$

$$\begin{array}{}\text{(18)}& G_h=\left\{\begin{array}{ll}\mathrm{1}-\frac{\mathrm{9.4}{\mathit{Ri}}_{\text{b}}}{\mathrm{1}+\mathrm{50}{\mathit{\kappa }}^{\mathrm{2}}(\left|{\mathit{Ri}}_{\text{b}}\right|z_{\text{r}}/z_{\text{o}}{)}^{\mathrm{0.5}}/{\ln }^{\mathrm{2}}(z_{\text{r}}/z_{\text{o}})}& {\mathit{Ri}}_{\text{b}}\le \mathrm{0}\\ \frac{\mathrm{1}}{(\mathrm{1}+\mathrm{4.7}{\mathit{Ri}}_{\text{b}}{)}^{\mathrm{2}}}& {\mathit{Ri}}_{\text{b}}>\mathrm{0}.\end{array}\right.\end{array}$$

3.2.2 Sedimentation velocity

Particle bins in the model are assumed to settle down with a sedimentation velocity ${\mathit{u}}_{\text{s}}=(\mathrm{0},\mathrm{0},-w_{\text{s}})$ equal to its terminal velocity:

$$\begin{array}{}\text{(19)}& w_{\text{s}}=\sqrt{{\displaystyle \frac{\mathrm{4}g({\mathit{\rho }}_{\text{p}}-{\mathit{\rho }}_{\text{a}})d}{\mathrm{3}{C}_{\text{d}}{\mathit{\rho }}_{\text{a}}}}},\end{array}$$

where ρ_a and ρ_p denote air and particle density, d is the particle equivalent diameter, and C_d is the drag coefficient that depends on the Reynolds number, $\mathit{Re}=d{u}_{\text{s}}/{\mathit{\nu }}_{\text{a}}$ (${\mathit{\nu }}_{\text{a}}={\mathit{\mu }}_{\text{a}}/{\mathit{\rho }}_{\text{a}}$ being the kinematic viscosity of air and μ_a its dynamic viscosity). For irregular particles, the drag coefficient C_d has to be obtained from experimental measurements. FALL3D-8.0 includes several parameterisations derived from laboratory results using natural and synthetic particles and that cover a wide range of particle sizes and shapes (characterised by sphericity, by circularity, or by some other model shape factor). Model options for the drag coefficient C_d include the following:

1. The GANSER model (Ganser, 1993) shows that

   $$\begin{array}{ll}\text{(20)}& {\displaystyle }{C}_{\text{d}}=& {\displaystyle \frac{\mathrm{24}}{\mathit{Re}{K}_{\mathrm{1}}}}{\displaystyle }\left\{\mathrm{1}+\mathrm{0.1118}{\left(\mathit{Re}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\right)}^{\mathrm{0.6567}}\right\}{\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{0.4305}{K}_{\mathrm{2}}}{\mathrm{1}+\frac{\mathrm{3305}}{\mathit{Re}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}}},\end{array}$$

   where $K_{\mathrm{1}}=\mathrm{3}/\left[\right(d_n/d)+\mathrm{2}{\mathrm{\Psi }}^{-\mathrm{0.5}}]$ and $K_{\mathrm{2}}={\mathrm{10}}^{\mathrm{1.8148}(-\text{Log}\mathrm{\Psi }{)}^{\mathrm{0.5743}}}$ are two shape factors, d_n is the average between the minimum and the maximum axis, and Ψ is the particle sphericity (Ψ=1 for spheres). For calculating the sphericity, it is practical to use the concepts of “operational” and “working sphericity”, Ψ_work introduced by Wadell (1933) and Aschenbrenner (1956), which are based on the determination of the volume and of the three dimensions of a particle, respectively:

   $$\begin{array}{}\text{(21)}& {\mathrm{\Psi }}_{\text{work}}=\mathrm{12.8}{\displaystyle \frac{(P^{\mathrm{2}}Q{)}^{\mathrm{1}/\mathrm{3}}}{\mathrm{1}+P(\mathrm{1}+Q)+\mathrm{6}\sqrt{\mathrm{1}+P^{\mathrm{2}}(\mathrm{1}+Q^{\mathrm{2}})}}},\end{array}$$

   with $P=d_{\mathrm{3}}/d_{\mathrm{2}}$, $Q=d_{\mathrm{2}}/d_{\mathrm{1}}$, where d₁ is the longest particle dimension, d₂ is the longest dimension perpendicular to d₁, and d₃ is the dimension perpendicular to both d₁ and d₂.

2. The PFEIFFER model (Pfeiffer et al., 2005), based on the interpolation of previous relationships by Walker et al. (1971) and Wilson and Huang (1979), shows

   $$\begin{array}{}\text{(22)}& C_{\text{d}}=\left\{\begin{array}{ll}\frac{\mathrm{24}}{\mathit{Re}}{\mathit{\phi }}^{-\mathrm{0.828}}+\mathrm{2}\sqrt{\mathrm{1}-\mathit{\phi }}& \mathit{Re}\le {\mathrm{10}}^{\mathrm{2}}\\ \mathrm{1}-\frac{\mathrm{1}-C_{\text{d}}{|}_{\mathit{Re}={\mathrm{10}}^{\mathrm{2}}}}{\mathrm{900}}({\mathrm{10}}^{\mathrm{3}}-\mathit{Re})& {\mathrm{10}}^{\mathrm{2}}\le \mathit{Re}\le {\mathrm{10}}^{\mathrm{3}}\\ \mathrm{1}& \mathit{Re}\ge {\mathrm{10}}^{\mathrm{3}},\end{array}\right.\end{array}$$

   where $\mathit{\phi }=(b+c)/\mathrm{2}a$ is the particle aspect ratio ($a\ge b\ge c$ denote the particle semi-axes).

3. The DIOGUARDI model (Dioguardi et al., 2018) demonstrates that

   $$\begin{array}{}\text{(23)}& \begin{array}{rl}{C}_{\text{d}}& ={\displaystyle \frac{\mathrm{24}}{\mathit{Re}}}{\left({\displaystyle \frac{\mathrm{1}-\mathit{\xi }}{\mathit{Re}}}+\mathrm{1}\right)}^{\mathrm{0.25}}\\ & +{\displaystyle \frac{\mathrm{24}}{\mathit{Re}}}\left(\mathrm{0.1806}{\mathit{Re}}^{\mathrm{0.6459}}\right){\mathit{\xi }}^{-{\mathit{Re}}^{\mathrm{0.08}}}{\displaystyle \frac{\mathrm{0.4251}}{\mathrm{1}+\frac{\mathrm{6880.95}}{{\mathit{Re}}^{\mathrm{2}}}{\mathit{\xi }}^{\mathrm{5.05}}}},\end{array}\end{array}$$

   where ξ is a particle shape factor (sphericity to circularity ratio), for which Dioguardi et al. (2018) suggested an empirical correlation with sphericity Ψ as ξ=0.83Ψ.

Note that, in any case, the terminal velocity w_s is defined by a triplet $(d,{\mathit{\rho }}_{\text{p}},\mathrm{\Psi })$. As a result, particles with similar values of the three parameters can be grouped within the same model bin.

3.2.3 Emissions

The emission source term for the ith bin ($S_i^{\text{e}}$ term in the bin equations; see Eq. 4) gives the mass per unit of time and volume (units of $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) released at each point (cell) of the computational domain. FALL3D-8.0 can generate and handle multiple types of emission sources, internally defined as a data structure made of n_p discrete points, each “tagged” with a time-varying position and bin emission rate (source strength) M_ip (in kg s⁻¹). As a result, $S_i^{\text{e}}$ in a model grid cell results from summing emissions from all point sources laying within the cell volume V:

$$\begin{array}{}\text{(24)}& S_i^{\text{e}}=\sum _{p=\mathrm{1}}^{n_{\text{p}}}{M}_{ip}/V.\end{array}$$

The total source strength M_o results from summing over all source points and bins, i.e.

$$\begin{array}{}\text{(25)}& M_{\text{o}}=\sum _{i=\mathrm{1}}^{n_{\text{b}}}\sum _{p=\mathrm{1}}^{n_{\text{p}}}{M}_{ip}=\sum _{i=\mathrm{1}}^{n_{\text{b}}}{M}_i.\end{array}$$

Table 5 Options available in FALL3D-8.0 for vertical distribution of mass and total source strength depending on the type of emission source S^e.

Table 5 summarises the different (exclusive) options available in FALL3D-8.0 for the emission term and related source strength.

1. The POINT option assumes that all mass is emitted from a single point (n_p=1) located at height z_t above ground level:

   $$\begin{array}{}\text{(26)}& M_{i\mathrm{1}}=\left\{\begin{array}{ll}{f}_i{M}_{\text{o}}& z=z_{\text{t}}\\ \mathrm{0}& z\ne z_{\text{t}},\end{array}\right.\end{array}$$

   where f_i is the ith bin mass fraction.

2. The HAT option defines a uniform vertical line of n_p source points spanning in height from z_b (bottom) to z_t (top) above the ground (i.e. with thickness z_t−z_b):

   $$\begin{array}{}\text{(27)}& M_{ip}=\left\{\begin{array}{ll}\frac{f_i{M}_{\text{o}}}{n_{\text{p}}}& z_{\text{b}}\le z\le z_{\text{t}}\\ \mathrm{0}& \text{otherwise}.\end{array}\right.\end{array}$$

   Note that this option includes as end-members the POINT option (if z_b=z_t) and a vertically uniform emission from ground to top (if z_b=0).

3. The SUZUKI option (Suzuki, 1983; Pfeiffer et al., 2005) assumes a mushroom-like vertical distribution of n_p emission points depending on two dimensionless parameters (A_s and λ_s):

   $$\begin{array}{}\text{(28)}& {\displaystyle }{M}_{ip}={\displaystyle \frac{f_i{M}_{\text{o}}}{n_{\text{p}}}}{\left[\left(\mathrm{1}-{\displaystyle \frac{z}{z_{\text{t}}}}\right)e^{A_{\text{s}}\left(\frac{z}{z_{\text{t}}}-\mathrm{1}\right)}\right]}^{{\mathit{\lambda }}_{\text{s}}}{\displaystyle }\mathrm{0}\le z\le z_{\text{t}}.\end{array}$$

   The Suzuki parameter A_s controls the vertical location of the maximum of the emission profile, whereas the parameter λ_s controls the distribution of the emitted mass around the maximum. When any of the previous source options are defined for volcanic plumes, it is useful to prescribe the total source strength (eruption mass flow rate) M_o in terms of the eruption column height H because this parameter is easier to be obtained from direct observations. To this purpose, FALL3D-8.0 includes two relationships that correlate M_o with H based on empirical observations and on 1-D plume model simulations, respectively. The first and simplest case considers the fit proposed by Mastin et al. (2009):

   $$\begin{array}{}\text{(29)}& M_{\text{o}}=a{H}^{\mathrm{4.15}},\end{array}$$

   where a=140.8 is a constant and H is the eruption column height expressed in kilometres above the eruptive vent. Alternatively, the 1-D model fit by Woodhouse et al. (2016) can also be used to provide M_o depending on the surrounding atmospheric conditions:

   $$\begin{array}{}\text{(30)}& M_{\text{o}}=c_{\text{o}}{N}^{\mathrm{3}}{H}^{\mathrm{4}}f\left(\stackrel{\mathrm{‾}}{W}\right),\end{array}$$

   where N is the Brunt–Väisälä frequency, c_o is a constant, and f is a function of the parameter $\stackrel{\mathrm{‾}}{W}=\mathrm{1.44}\stackrel{\mathrm{‾}}{\mathit{\gamma }}/N$ given by

   $$\begin{array}{}\text{(31)}& f\left(\stackrel{\mathrm{‾}}{W}\right)={\left({\displaystyle \frac{\mathrm{1}+b\stackrel{\mathrm{‾}}{W}+c{\stackrel{\mathrm{‾}}{W}}^{\mathrm{2}}}{\mathrm{1}+a\stackrel{\mathrm{‾}}{W}}}\right)}^{\mathrm{4}},\end{array}$$

   with the coefficients $a=\mathrm{0.87}+\mathrm{0.05}{\mathit{\beta }}_{\text{e}}/{\mathit{\alpha }}_{\text{e}}$, $b=\mathrm{1.09}+\mathrm{0.32}{\mathit{\beta }}_{\text{e}}/{\mathit{\alpha }}_{\text{e}}$, and $c=\mathrm{0.06}+\mathrm{0.03}{\mathit{\beta }}_{\text{e}}/{\mathit{\alpha }}_{\text{e}}$, with α_e and β_e the radial and cross-wind plume entrainment coefficients (Costa et al., 2016b) and $\stackrel{\mathrm{‾}}{\mathit{\gamma }}$ the mean shear rate of wind. The constant c_o in Eq. (30) depends on the conditions at the vent:

   $$\begin{array}{}\text{(32)}& c_{\text{o}}={\displaystyle \frac{\mathrm{0.35}{\mathit{\alpha }}_{\text{e}}^{\mathrm{2}}{\mathit{\rho }}_a{C}_{\text{a}}{T}_{\text{a}}}{g\left[(C_{\text{v}}{n}_{\text{o}}+C_{\text{s}}(\mathrm{1}-n_{\text{o}}\left)\right)T_{\text{o}}-C_a{T}_{\text{a}}\right]}},\end{array}$$

   with C_s, C_v and C_a the specific heat capacities at constant pressure of the solid pyroclasts, gas phase and air, respectively, T_a and T_o the air and vent magma mixture temperatures, and n₀ the gas mass fraction.

4. The PLUME option, valid only for volcanic plumes, uses the FPLUME-1.0 model (Folch et al., 2016) embedded in FALL3D-8.0. FPLUME-1.0 is a steady-state 1-D cross-section-averaged eruption column model based on the buoyant plume theory (BPT). The model accounts for plume bending by wind, entrainment of ambient moisture, effects of water phase changes, particle fallout and re-entrainment, and a model for wet aggregation of ash particles in the presence of liquid water or ice. As opposed to the previous cases, the PLUME source option automatically determines a bin-dependent vertical distribution of mass and computes height from M_o, or vice versa by solving an inverse problem.

5. The RESUSPENSION option considers the remobilisation and resuspension by wind of soil particles (e.g. mineral dust or volcanic ash previously deposited on the ground). Up to three different emission schemes are available in FALL3D-8.0 to obtain the vertical flux of suspended particles, from which M_o is obtained by multiplying by the associated surface cell area (see  Folch et al., 2014, for details). Typically, the emission schemes for mineral dust are formulated in terms of the friction velocity. For example, emission scheme 1 (Westphal et al., 1987) considers

   $$\begin{array}{}\text{(33)}& F_{\text{V}}=\left\{\begin{array}{ll}\mathrm{0}& u_{*}<u_{*t}\\ {\mathrm{10}}^{-\mathrm{5}}{u}_{*}^{\mathrm{4}}& u_{*}\ge u_{*t},\end{array}\right.\end{array}$$

   where F_V is the vertical flux (in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), occurring only above a (constant) threshold friction velocity u_*t (u_* given in m s⁻¹). An important limitation of Eq. (33) is that the vertical flux does not depend on either particle size or soil moisture. However, despite its simplicity, this parameterisation can be useful when information on soil characteristics (e.g. particle sizes and densities, moisture, roughness) is unavailable or poorly constrained. Emission scheme 2 (Marticorena and Bergametti, 1995; Marticorena et al., 1997) considers

   $$\begin{array}{}\text{(34)}& F_{\text{V}}=\left\{\begin{array}{ll}\mathrm{0}& u_{*}<u_{*t}\\ S_{\text{c}}{\displaystyle \frac{{\mathit{\rho }}_{\text{a}}}{g}}{u}_{*}^{\mathrm{3}}\left(\mathrm{1}-{\displaystyle \frac{u_{*t}^{\mathrm{2}}}{u_{*}^{\mathrm{2}}}}\right)\left(\mathrm{1}+{\displaystyle \frac{u_{*t}}{u_{*}}}\right)& u_{*}\ge u_{*t},\end{array}\right.\end{array}$$

   where the experimental coefficient S_c (in cm⁻¹) depends on the amount of available fine particles in the soil, and the threshold friction velocity is given by

   $$\begin{array}{}\text{(35)}& \begin{array}{rl}& u_{*t}=\\ & \left\{\begin{array}{ll}\frac{\mathrm{0.129}K}{(\mathrm{1.928}{\mathit{Re}}^{\mathrm{0.092}}-\mathrm{1}{)}^{\mathrm{0.5}}}& \mathrm{0.03}<\mathit{Re}\le \mathrm{10}\\ \mathrm{0.129}K(\mathrm{1}-\mathrm{0.0858}{e}^{-\mathrm{0.0617}(\mathit{Re}-\mathrm{10})})& \mathit{Re}>\mathrm{10},\end{array}\right.\end{array}\end{array}$$

   with $K=\sqrt{\frac{{\mathit{\rho }}_{\text{p}}gd}{{\mathit{\rho }}_{\text{a}}}\left(\mathrm{1}+\frac{\mathrm{0.006}}{{\mathit{\rho }}_{\text{p}}g{d}^{\mathrm{2.5}}}\right)}$ and $\mathit{Re}=\mathrm{1331}\times d^{\mathrm{1.56}}$ (the lower bound of the fit corresponds to particles of ≈ 10 µm in size). Note that in Eq. (35), ρ_p and ρ_a are particle and air densities (expressed in g cm⁻³), g is gravity (in cm s⁻²), d is the particle size (in centimetres), Re is the Reynolds number parameterised as a function of the particle size, and u_*t is given in cm s⁻¹.

   Finally, emission scheme 3 (Shao et al., 1993; Shao and Leslie, 1997; Shao and Lu, 2000) considers that the uplift from surface of the fine fraction of soil particles is controlled by the bombardment of saltating particles of larger sizes (≥ 63 µm), which breaks the cohesive forces of smaller particles.

   Based on theoretical and experimental results, Shao et al. (1993) found an expression for the vertical flux of dust particles of size d ejected by the impact of saltating particles of size d_s:

   $$\begin{array}{}\text{(36)}& F_{\text{V}}(d,d_{\text{s}})={\displaystyle \frac{{\mathit{\alpha }}_{\text{s}}(d,d_{\text{s}})}{u_{*t}^{\mathrm{2}}\left(d\right)}}{F}_{\text{H}}\left(d_{\text{s}}\right),\end{array}$$

   where α_s (in m s⁻²) is the coefficient of sandblasting efficiency determined experimentally (Shao and Leslie, 1997) and F_H is the horizontal flux (in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) of saltating particles:

   $$\begin{array}{}\text{(37)}& F_{\text{H}}\left(d_{\text{s}}\right)=\left\{\begin{array}{ll}\mathrm{0}& u_{*}<u_{*t}\left(d_{\text{s}}\right)\\ c_{\text{o}}\frac{{\mathit{\rho }}_{\text{a}}{u}_{*}^{\mathrm{3}}}{g}\left(\mathrm{1}-\frac{u_{*t}^{\mathrm{2}}\left(d_{\text{s}}\right)}{u_{*}^{\mathrm{2}}}\right)& u_{*}\ge u_{*t}\left(d_{\text{s}}\right),\end{array}\right.\end{array}$$

   where c_o is an empirical dimensionless constant close to 1. In this scheme, the threshold friction velocity u_*t(d) is given by

   $$\begin{array}{}\text{(38)}& u_{*ts}=\sqrt{\mathrm{0.0123}\left({\displaystyle \frac{{\mathit{\rho }}_{\text{p}}gd}{{\mathit{\rho }}_{\text{a}}}}+{\displaystyle \frac{\mathit{\gamma }}{{\mathit{\rho }}_{\text{a}}d}}\right)},\end{array}$$

   where γ is an experimental parameter ranging between $\mathrm{1.65}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ kg s² (a value of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ kg s² is assumed in FALL3D-8.0).

3.2.4 Deposition mechanisms

In FALL3D-8.0, dry and wet deposition mechanisms can be activated for any type of bin below a certain particle/aerosol size.

Dry deposition on the ground is imposed prescribing the deposition velocity through a Robin boundary condition in Eq. (2). FALL3D-8.0 admits two dry deposition parameterisations which describe the vertical depositional fluxes by Brownian diffusion and inertial impaction, parameterised through the Schmidt and the Stokes numbers, respectively. The first option considers the mass-consistent formulation proposed by Venkatram and Pleim (1999):

$$\begin{array}{}\text{(39)}& \left|{\mathit{u}}_{\text{d}}\right|=w_{\text{s}}+{\displaystyle \frac{w_{\text{s}}}{\mathrm{1}-e^{-(r_{\text{a}}+r_{\text{b}})w_{\text{s}}}}}\approx w_{\text{s}}+{\displaystyle \frac{\mathrm{1}}{r_{\text{a}}+r_{\text{b}}}},\end{array}$$

where r_a describes the effects of aerodynamic resistance and r_b the quasi-laminar resistance (e.g. Brandt et al., 2002, and references therein). The aerodynamic resistance r_a can be calculated as

$$\begin{array}{}\text{(40)}& r_{\text{a}}={\displaystyle \frac{\mathrm{1}}{k{u}_{*}}}\left[\ln \left({\displaystyle \frac{z}{z_{\text{o}}}}\right)-{\mathit{\varphi }}_{\text{h}}\left({\displaystyle \frac{z}{L}}\right)\right],\end{array}$$

with z_o denoting the ground roughness height and ϕ_h the atmospheric stability function for temperature given by Eq. (12). The quasi-laminar resistance r_b can be expressed in terms of the Schmidt number $Sc=\mathit{\nu }/D$ and Stokes number $\mathit{St}=w_s{u}_{*}^{\mathrm{2}}/\left(g\mathit{\nu }\right)$ (with ν kinematic viscosity of air and D molecular diffusivity of particles) (e.g. Brandt et al., 2002):

$$\begin{array}{}\text{(41)}& r_{\text{b}}={\displaystyle \frac{\mathrm{1}}{u_{*}\left(S{c}^{-\mathrm{2}/\mathrm{3}}+{\mathrm{10}}^{-\mathrm{3}/\mathit{St}}\right)}}.\end{array}$$

The second option is that proposed by Feng (2008), which essentially differs from Eq. (39) in the estimation of r_b:

$$\begin{array}{}\text{(42)}& \left|{\mathit{u}}_{\text{d}}\right|=w_{\text{s}}+{\displaystyle \frac{\mathrm{1}}{r_{\text{a}}+\mathrm{1}/(u_{*}{c}_{\mathrm{1}}{e}^{-\mathrm{0.5}{\left[({\mathit{Re}}^{*}-c_{\mathrm{2}})/c_{\mathrm{3}}\right]}^{\mathrm{2}}}+a{u}_{*}^b)}},\end{array}$$

where c₁=0.0226, c₂=40300 and c₃=15330 are dimensionless parameterisation constants, Re^* is the Reynolds number (computed with the friction velocity u_*), and a and b are coefficients that depend on the particle size and surface characteristics. Note that Feng (2008) gives a and b best-fit values for seven land use categories and four aerosol size modes: nuclei (up to 0.1 µm), accumulation (up to 2.5 µm), coarse (up to 10 µm) and giant (up to 100 µm). A cut-off is assumed above this size because the sedimentation velocity term w_s dominates, and therefore the dry deposition contribution can be neglected.

Wet deposition mechanisms in FALL3D-8.0 are assumed to occur only within the ABL, and the corresponding sink term in Eq. (3) is parameterised as

$$\begin{array}{}\text{(43)}& I^{\text{w}}=\mathrm{\Lambda }c,\end{array}$$

where Λ differs for in-cloud (ic) and below-cloud (bc) sinks. For below-cloud scavenging (precipitation), Λ_bc is estimated from the total precipitation rate as (e.g. Brandt et al., 2002; Jung and Shao, 2006)

$$\begin{array}{}\text{(44)}& {\mathrm{\Lambda }}_{\text{bc}}=a{P}^b,\end{array}$$

where P is the precipitation rate (in mm h⁻¹), and $a=\mathrm{8.4}\times {\mathrm{10}}^{-\mathrm{5}}$ and b=0.79 are two empirical constants. For in-cloud scavenging (rainout), the model considers a parameterisation based on the atmospheric relative humidity (RH; in %) as in Brandt et al. (2002):

$$\begin{array}{}\text{(45)}& {\mathrm{\Lambda }}_{ic}=\left\{\begin{array}{ll}\mathrm{0}& \text{for}\mathrm{RH}<{\mathrm{RH}}_{\text{t}}\\ A_{\text{RH}}\frac{\mathrm{RH}-{\mathrm{RH}}_{\text{t}}}{{\mathrm{RH}}_{\text{s}}-{\mathrm{RH}}_{\text{t}}}& \text{for}\mathrm{RH}\ge {\mathrm{RH}}_{\text{t}},\end{array}\right.\end{array}$$

with $A_{\text{RH}}=\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{5}}$, RH_t=80 % (threshold value) and RH_s=100 % (saturation value). Two critical particle cut-off sizes of 100 and 1 µm are assumed for below- and in-cloud scavenging, respectively.

3.2.5 Gravity spreading of the umbrella region

Large explosive volcanic eruptions can generate gravity-driven transport mechanisms that dominate over passive transport close to the vent and cause a radial spreading of the cloud (e.g. Woods and Kienle, 1994; Sparks et al., 1997). In order to simulate this mechanism, FALL3D-8.0 includes a gravity current model (see Costa et al., 2013, and the erratum published in June 2019). This option consists on adding a radial velocity field to the background wind, so that contributions from both passive and density-driven mechanisms are accounted for. The added radial wind is centred above the eruptive vent in the umbrella region and extended up to a radius R given by

$$\begin{array}{}\text{(46)}& R={\left({\displaystyle \frac{\mathrm{3}{\mathit{\lambda }}_{\text{g}}Nq}{\mathrm{2}\mathit{\pi }}}\right)}^{\mathrm{1}/\mathrm{3}}{t}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where t is time since eruption onset, λ_g is an empirical constant constrained to ≈0.2 from direct numerical simulations (Suzuki and Koyaguchi, 2009), N is the Brunt–Väisälä frequency, and q is the volumetric flow rate into the umbrella region, estimated as (Morton et al., 1956; Suzuki and Koyaguchi, 2009; Costa et al., 2013, as correct in the erratum published in 2019)

$$\begin{array}{}\text{(47)}& q={\displaystyle \frac{c{k}_{\text{e}}^{\mathrm{1}/\mathrm{2}}{M}_{\text{o}}^{\mathrm{3}/\mathrm{4}}}{N^{\mathrm{5}/\mathrm{4}}}},\end{array}$$

where M_o is the total source strength (i.e. mass eruption rate), k_e is the air entrainment coefficient, and c is a constant that from varies from tropical to midlatitude/polar locations:

$$\begin{array}{}\text{(48)}& c=\left\{\begin{array}{ll}\mathrm{0.43}& {\mathrm{m}}^{\mathrm{3}}{\mathrm{kg}}^{-\mathrm{3}/\mathrm{4}}{\mathrm{s}}^{-\mathrm{3}/\mathrm{2}}\text{for tropical}\\ \mathrm{0.87}& {\mathrm{m}}^{\mathrm{3}}{\mathrm{kg}}^{-\mathrm{3}/\mathrm{4}}{\mathrm{s}}^{-\mathrm{3}/\mathrm{2}}\text{midlatitude/polar}.\end{array}\right.\end{array}$$

Given the radius R, the radial velocity field as a function of distance r is calculated as (Costa et al., 2013)

$$\begin{array}{}\text{(49)}& u_{\text{r}}\left(r\right)={\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}{u}_{\text{r}}\left(R\right){\displaystyle \frac{R}{r}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\displaystyle \frac{r^{\mathrm{2}}}{R^{\mathrm{2}}}}\right)(\mathrm{0}\le r\le R),\end{array}$$

where u_r(R) is the front velocity:

$$\begin{array}{}\text{(50)}& u_{\text{r}}\left(R\right)={\left({\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\text{g}}Nq}{\mathrm{3}\mathit{\pi }}}\right)}^{\mathrm{1}/\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\sqrt{R}}}.\end{array}$$

In order to avoid sudden jumps at the gravity current front, FALL3D-8.0 interpolates the front velocity u_r(R) with far field wind velocity using an exponential decay function of the cloud thickness h_c as

$$\begin{array}{}\text{(51)}& \exp [-r/(\mathrm{4}{h}_{\text{c}}\left)\right],\end{array}$$

where r is the distance from the current front.

Table 6 List of radionuclides implemented in the model. Table shows half life, decay rate (in s⁻¹) and resulting child product.

3.2.6 Aggregation

Aggregation of tephra particles can occur inside the eruptive columns or even downwind in ash clouds during atmospheric dispersion, thereby affecting the sedimentation dynamics and deposition of volcanic ash. FALL3D-8.0 includes some simple a priori aggregation options and a wet aggregation model (Costa et al., 2010; Folch et al., 2010) that can be activated for tephra bins. The a priori options consist of user-defined or empirically based predefined fractions of aggregating classes being transferred to one or more class of aggregates at the source points (i.e. aggregation is performed before transport). In contrast, in the wet aggregation model, ash particles aggregate on a single effective class of diameter d_A; i.e. aggregation only affects tephra bins with diameter smaller than d_A, typically in the range 100–300 µm. This option can run embedded in FPLUME-1.0 or as a stand-alone module. Consider a tephra grain size distribution in which k particle bins can aggregate. Then, the aggregation model defines the source (S^a) and sink (I^a) bin terms for the corresponding k+1 bins as

$$\begin{array}{}\text{(52)}& \left\{\begin{array}{ll}{S}_{k+\mathrm{1}}^{\text{a}}={\sum }_{j=\mathrm{1}}^k{I}_j^{\text{a}}& \\ I_j^{\text{a}}=\frac{\mathit{\pi }{d}_j{\mathit{\rho }}_j}{\mathrm{6}}{\dot{n}}_j& j=\mathrm{1}:k,\end{array}\right.\end{array}$$

where d_j(<d_A) and ρ_j are, respectively, the diameter and density of particles in bin j, and ${\dot{n}}_j$ is the number of particles per unit volume and time that aggregate. The model assumes that this is proportional to the total particle decay per unit volume and time ${\dot{n}}_{\text{tot}}$, i.e.

$$\begin{array}{}\text{(53)}& {\dot{n}}_j\approx {\displaystyle \frac{N_j}{{\sum }_{i=\mathrm{1}}^k{N}_i}}{\dot{n}}_{\text{tot}},\end{array}$$

where

$$\begin{array}{}\text{(54)}& N_j=k_{\text{f}}{\left({\displaystyle \frac{d_A}{d_j}}\right)}^{d_f}\end{array}$$

is the number of primary particles of diameter d_j in an aggregate of diameter d_A, k_f is a fractal pre-factor (k_f≈1), and D_f is the fractal exponent (D_f≤3). The model estimates the total particle decay per unit time ${\dot{n}}_{\text{tot}}$ integrating the coagulation kernel over all particle sizes, depending on the sticking efficiency times a collision frequency function which accounts for Brownian motion, collision due to turbulence as a result of inertial effects, laminar and turbulent fluid shear, and differential sedimentation (see Costa et al., 2010; Folch et al., 2016, for details).

3.2.7 Radioactive decay

FALL3D-8.0 can handle the fate of radioactive material dispersed from accidental releases (e.g. Brandt et al., 2002; Leelössya et al., 2018). Five common species of radionuclides have been implemented: cesium (¹³⁴Cs, ¹³⁷Cs) and iodine (¹³¹I), which decay to stable isotopes, and ⁹⁰Sr, which decays to the unstable isotope ⁹⁰Y (see Table 6).

Radionuclide species need to specify the source ($S_n^{\text{r}}$) and the sink ($I_n^{\text{r}}$) terms in Eq. (3) associated with the radioactive production or decay of the isotope n, respectively. The radioactive decay term indicates the mass per unit volume of the isotopes of type n that decays per unit time:

$$\begin{array}{}\text{(55)}& I_n^{\text{r}}=k_{\text{r}}{c}_n,\end{array}$$

where c_n is concentration (expressed in kg m⁻³) and k_r a constant specific of each isotope (decay rate) that can be calculated from the radioactive element half life t_1∕2 as $k_{\text{r}}=\ln \left(\mathrm{2}\right)/t_{\mathrm{1}/\mathrm{2}}$. Values of t_1∕2 and k_r for common radionuclides are reported in Table 6. Note that the decay term is more relevant for isotopes with short half lives, e.g. for ¹³¹I, which has $t_{\mathrm{1}/\mathrm{2}}\simeq \mathrm{8}$ d (Brandt et al., 2002; Leelössya et al., 2018).

The radioactive decay term of the isotope n constitutes a sink $I_n^{\text{r}}$ for the isotope itself. However, the decay of the isotope m, father of isotope n, constitutes a source $S_n^{\text{r}}$ for the isotope n:

$$\begin{array}{}\text{(56)}& S_n^{\text{r}}=p_{mn}{I}_m^{\text{r}},\end{array}$$

where p_mn is the relative probability of decay of the isotope m to the isotope n. If the isotope m has only one child n, the relative probability of the branch m↦n is p_mn=1. Note from Table 6 that ¹³⁴Cs, ¹³⁷Cs and ¹³¹I decay to stable isotopes, whereas ⁹⁰Sr decays to ⁹⁰Y which, in turn, is unstable and decays to the stable isotope ⁹⁰Zr. The production rate $S_n^{\text{r}}$ of ⁹⁰Y is therefore equivalent to the decay rate of ⁹⁰Sr.

In FALL3D-8.0, the radioactive decay is implemented by firstly transporting the radionuclides for a time step Δt and then by evaluating the decay during the same time step. For radionuclides that decay to a stable isotope (¹³⁴Cs, ¹³⁷Cs and ¹³¹I), it is considered that after a time step Δt the concentration decreases as

$$\begin{array}{}\text{(57)}& c(t+\mathrm{\Delta }t)=c\left(t\right)e^{-k_{\text{r}}\mathrm{\Delta }t},\end{array}$$

whereas for decay to an unstable isotope (yttrium), the concentration varies as

$$\begin{array}{}\text{(58)}& c_{\text{Y}}(t+\mathrm{\Delta }t)=\left[c_{\text{Sr}}\left(t\right)(\mathrm{1}-e^{-k_{\text{Sr}}\mathrm{\Delta }t})+c_{\text{Y}}\left(t\right)\right]e^{-k_{\text{Y}}\mathrm{\Delta }t},\end{array}$$

where c_Y and c_Sr are, respectively, the concentrations of ⁹⁰Y and ⁹⁰Sr, and k_Y and k_Sr are the corresponding decay rates.

Table 7 Characteristics of sensors for ash and SO₂ detection aboard the new generation of geostationary satellites. Table courtesy of Andrew Prata.

3.3 Data insertion

Instrumentation aboard the new generation of geostationary satellites provides an unprecedented level of spatial resolution and temporal frequency (2 to 4 km pixel size and 10 to 15 min observation period; see Table 7), yielding a quasi-global coverage considering the overlap of different existing platforms. This is very suitable for high-resolution model data assimilation and related uncertainty quantification, as well as to implement ensemble-based dispersal forecast systems. These aspects are still under development and hopefully will be part of next FALL3D model distributions. However, FALL3D-8.0 already includes the possibility of initialising a model run from satellite retrievals. This option, known as data insertion, is typically used in dispersal of volcanic ash and aerosols (mainly SO₂) in order to reduce model uncertainties coming from the eruption source term.

Satellite retrievals giving cloud column mass of fine volcanic ash and aerosols can be furnished to the model together with values of cloud thickness, with the latter needed in order to compute initial concentration (in kg m⁻³) from column mass (in kg m⁻²). In the model initialisation step, gridded satellite data are interpolated into the model grid imposing conservation of mass when concentration values are computed for each model grid cell, i.e. ensuring that the resulting column mass in the model (computed concentration multiplied by cloud thickness) equals that of satellite data over the same cell area. Examples showing how data insertion improves model accuracy are given in the companion paper (Prata et al., 2020).

4 Numerical implementation

4.1 Coordinate mappings and scaling

Consider the ADS equation (Eq. 1) written in a Cartesian system of coordinates $(x,y,z)$ assuming a diffusion tensor as in Eq. (5) and a sedimentation velocity ${\mathit{u}}_{\text{s}}=(\mathrm{0},\mathrm{0},-w_{\text{s}})$ aligned with the vertical coordinate z:

$$\begin{array}{}\text{(59)}& \begin{array}{rl}& {\displaystyle \frac{\partial c}{\partial t}}+{\displaystyle \frac{\partial \left(cu\right)}{\partial x}}+{\displaystyle \frac{\partial \left(cv\right)}{\partial y}}+{\displaystyle \frac{\partial \left(cw\right)}{\partial z}}-{\displaystyle \frac{\partial \left(c{w}_{\text{s}}\right)}{\partial z}}\\ & -{\displaystyle \frac{\partial }{\partial x}}\left(k_{\text{h}}{\displaystyle \frac{\partial c}{\partial x}}\right)-{\displaystyle \frac{\partial }{\partial y}}\left(k_{\text{h}}{\displaystyle \frac{\partial c}{\partial y}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left(k_{\text{v}}{\displaystyle \frac{\partial c}{\partial z}}\right)=S-I.\end{array}\end{array}$$

It is straightforward to discretise the above equation in a “brick-like” computational domain Ω_c using a structured regular (i.e. equally spaced) mesh, although the regularity condition is typically relaxed across the vertical direction so that the vertical grid resolution increases close to ground, where higher gradients are expected. In order to use other coordinate systems, Eq. (59) can be written on a generalised orthogonal system of coordinates $(X_{\mathrm{1}},X_{\mathrm{2}},X_{\mathrm{3}})$ (e.g. Toon et al., 1988; Byun and Schere, 2005):

$$\begin{array}{}\text{(60)}& \begin{array}{rl}& {\displaystyle \frac{\partial C}{\partial t}}+{\displaystyle \frac{\partial \left(CU\right)}{\partial X_{\mathrm{1}}}}+{\displaystyle \frac{\partial \left(CV\right)}{\partial X_{\mathrm{2}}}}+{\displaystyle \frac{\partial \left(CW\right)}{\partial X_{\mathrm{3}}}}-{\displaystyle \frac{\partial \left(C{W}_{\text{s}}\right)}{\partial X_{\mathrm{3}}}}\\ & -{\displaystyle \frac{\partial }{\partial X_{\mathrm{1}}}}\left(K_{\mathrm{1}}{\displaystyle \frac{\partial C}{\partial X_{\mathrm{1}}}}\right)-{\displaystyle \frac{\partial }{\partial X_{\mathrm{2}}}}\left(K_{\mathrm{2}}{\displaystyle \frac{\partial C}{\partial X_{\mathrm{2}}}}\right)-{\displaystyle \frac{\partial }{\partial X_{\mathrm{3}}}}\left(K_{\mathrm{3}}{\displaystyle \frac{\partial C}{\partial X_{\mathrm{3}}}}\right)\\ & =S^{*}-I^{*},\end{array}\end{array}$$

where C is the scaled concentration, $(U,V,W)$ are the scaled wind components, $(K_{\mathrm{1}},K_{\mathrm{2}},K_{\mathrm{3}})$ are the scaled diffusion coefficients, and S^* and I^* are the scaled source and sink terms. The implementation of a generalised equation like Eq. (60) presents two major advantages. On one hand, the generalised equation reads formally equal to that in Cartesian coordinates, so that little computational penalty exists to map physical domains (e.g. accounting for Earth's curvature and topography) to a “brick-like” computational domain (see Fig. 1) by using coordinate-dependent horizontal and vertical mappings. On the other hand, a generalised form simplifies the structure and implementation of the code because the model can be solved on various horizontal (Cartesian, spherical, Mercator, polar stereographic, etc.) and vertical (terrain-following, σ-coordinate, etc.) coordinate systems using only one solving routine. To this purpose, one needs first to scale the model coordinates and some terms in the equation using adequate mapping and scaling factors, then solve for the scaled concentration C in the regular computational domain Ω_c (as in Cartesian coordinates) and, finally, transform the scaled concentration back to the original one.

Figure 1 Mapping of coordinates from the physical domain (a) to a brick-like computational domain (b) along the horizontal and vertical dimensions. Horizontal mapping corrects for Earth's curvature; vertical mapping accounts for topography.

In general and given two orthogonal coordinate systems, $(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})$ and $(X_{\mathrm{1}},X_{\mathrm{2}},X_{\mathrm{3}})$, coordinate mapping factors are given by the terms m_ij of the Jacobian transformation matrix 𝕄:

$$\begin{array}{}\text{(61)}& m_{ij}={\displaystyle \frac{\partial x_i}{\partial X_j}},\end{array}$$

but, for the transformations considered here, 𝕄 will always be diagonal with three non-zero components (m₁, m₂ and m₃). For example, in the horizontal transformation to spherical Earth surface coordinates (λ,ϕ), one has

$$\begin{array}{}\text{(62)}& \begin{array}{rlrl}& x=R_{\text{e}}\sin \mathit{\gamma }\mathit{\lambda }& \equiv & \sin \mathit{\gamma }{X}_{\mathrm{1}}\\ & y=R_{\text{e}}\mathit{\varphi }& \equiv & X_{\mathrm{2}},\end{array}\end{array}$$

where R_e is the radius of the Earth, and λ, ϕ and γ are the longitude, latitude and colatitude, respectively (in Rad). Trivially, this transformation yields to m₁=sin γ and m₂=1. Table 8 gives horizontal mapping factors for different coordinate systems (Toon et al., 1988). In most practical cases, FALL3D-8.0 simulations only consider spherical coordinates, but other options are in principle possible. For vertical transformations, FALL3D-8.0 incorporates a new σ-coordinate system with linear decay (Gal-Chen and Somerville, 1975) in which

$$\begin{array}{}\text{(63)}& z={\displaystyle \frac{H_{\text{T}}-h}{H_{\text{T}}}}{X}_{\mathrm{3}}+h{x}_{\mathrm{3}}\in [\mathrm{0},H_{\text{T}}],\end{array}$$

where h(x,y) is the terrain height and H_T the height of the top of the computational domain. In the σ-coordinate system, the influence of the terrain decreases linearly with height, from terrain-following at the surface (X₃=0) to a rigid lid at the top of the computational domain (X₃=H_T). This option has been added to partially correct numerical oscillations in the previous terrain-following model mapping ($z=X_{\mathrm{3}}+h$), which can appear near the surface in the case of flows over mountain ranges and propagate upwards (Schar et al., 2002).

Table 8 Horizontal mapping factors $(x,y)\leftarrow (X_{\mathrm{1}},X_{\mathrm{2}})$ for different coordinate systems where λ is longitude, ϕ latitude, R the radius of the Earth, γ the colatitude, and ϕ_o the latitude at which the projection is true. Mercator and polar stereographic projections use Cartesian coordinates projected into the Earth's surface, i.e. account for curvature. Regular Cartesian coordinates should be used only for local domains, where the Earth's curvature can be neglected.

Table 9 Vertical mapping factors z←X₃ for different coordinate systems, where h(x,y) is topography and H is the top of the computational domain.

Table 10 Scaling factors for the different terms in the generalised coordinate ADS equation (Eq. 60).

Table 9 gives the vertical mapping factors for the different coordinate systems available in FALL3D-8.0. Once defined, these coordinate mapping factors are used to scale the variables and parameters that appear in the generalised equation (Eq. 60). The scaling of scalar quantities is straightforward since it only involves the Jacobian determinant of the transformation, e.g. $C=\mid \mathbb{M}\mid c=m_{\mathrm{1}}{m}_{\mathrm{2}}{m}_{\mathrm{3}}c$ for concentration and so on. Note that for a volume one has dV=m₁m₂m₃dv, so that the mass comprised in a cell dX₁dX₂dX₃ of the computational domain is equal to that in the transformed cell of the physical domain. The horizontal velocity components are trivially scaled as

$$\begin{array}{}\text{(64)}& \left\{\begin{array}{l}U=\frac{\text{d}{X}_{\mathrm{1}}}{\text{d}t}=u\frac{\partial X_{\mathrm{1}}}{\partial x}+v\frac{\partial X_{\mathrm{1}}}{\partial y}+w\frac{\partial X_{\mathrm{1}}}{\partial z}=u/m_{\mathrm{1}}\\ V=\frac{\text{d}{X}_{\mathrm{2}}}{\text{d}t}=u\frac{\partial X_{\mathrm{2}}}{\partial x}+v\frac{\partial X_{\mathrm{2}}}{\partial y}+w\frac{\partial X_{\mathrm{2}}}{\partial z}=v/m_{\mathrm{2}},\end{array}\right.\end{array}$$

whereas for the vertical component one has to consider that, in terrain-following coordinate systems, the coordinate X₃ depends also on (x,y) through the terrain elevation h(x,y):

$$\begin{array}{}\text{(65)}& {\displaystyle }W& {\displaystyle }={\displaystyle \frac{\text{d}{X}_{\mathrm{3}}}{\text{d}t}}=u{\displaystyle \frac{\partial X_{\mathrm{3}}}{\partial x}}+v{\displaystyle \frac{\partial X_{\mathrm{3}}}{\partial y}}+w{\displaystyle \frac{\partial X_{\mathrm{3}}}{\partial z}}\text{(66)}& {\displaystyle }& {\displaystyle }=u{\displaystyle \frac{\partial h}{\partial x}}{\displaystyle \frac{\partial X_{\mathrm{3}}}{\partial h}}+v{\displaystyle \frac{\partial h}{\partial y}}{\displaystyle \frac{\partial X_{\mathrm{3}}}{\partial h}}+w/m_{\mathrm{3}}.\end{array}$$

Basic algebra manipulation yields to the scaling factors shown in Table 10 for the different vertical coordinate systems. Note that the expression above implicitly contains the correction for topography in the vertical velocity. Finally, scaling factors for diffusion coefficients can also be obtained after some manipulation (Toon et al., 1988; Byun and Schere, 2006).

4.2 Discretisation and solving algorithm

FALL3D solves for each model bin the 3-D generalised equation (Eq. 60) using a fractional step method, which splits the equation along each spatial direction (Folch et al., 2009):

$$\begin{array}{}\text{(67)}& \begin{array}{rl}& {\displaystyle \frac{\partial C}{\partial t}}=S^{*}-I^{*}\\ & {\displaystyle \frac{\partial C}{\partial t}}+{\displaystyle \frac{\partial \left(CU\right)}{\partial X_{\mathrm{1}}}}-{\displaystyle \frac{\partial }{\partial X_{\mathrm{1}}}}\left(K_{\mathrm{1}}{\displaystyle \frac{\partial C}{\partial X_{\mathrm{1}}}}\right)=\mathrm{0}\\ & {\displaystyle \frac{\partial C}{\partial t}}+{\displaystyle \frac{\partial \left(CV\right)}{\partial X_{\mathrm{2}}}}-{\displaystyle \frac{\partial }{\partial X_{\mathrm{2}}}}\left(K_{\mathrm{2}}{\displaystyle \frac{\partial C}{\partial X_{\mathrm{2}}}}\right)=\mathrm{0}\\ & {\displaystyle \frac{\partial C}{\partial t}}+{\displaystyle \frac{\partial \left[C(W-W_{\text{s}})\right]}{\partial X_{\mathrm{3}}}}-{\displaystyle \frac{\partial }{\partial X_{\mathrm{3}}}}\left(K_{\mathrm{3}}{\displaystyle \frac{\partial C}{\partial X_{\mathrm{3}}}}\right)=\mathrm{0},\end{array}\end{array}$$

with the solving order of each resulting one-dimensional ADS equation being permuted in each successive time step in order to avoid any privileged direction. One advantage of this splitting strategy is that, even if each one-dimensional equation is solved explicitly in time, the result is a semi-implicit scheme that adds stability.

Figure 2 Arakawa D grid on a 2-D computational domain limited by the bold line. Scalar quantities are stored/computed at the centre of the cells (empty circles), whereas the u and v velocity components are staggered at its respective cell faces. One row of ghost cells is shown in grey for reference; the actual number of ghost cells needed depends on the numerical stencil (order of the solving algorithm). The 3-D grid is formed as a succession of 2-D layers, with the w velocity components at the bottom/top faces of the cell.

In FALL3D-8.0, the “brick-like” computational domain is discretised using a variation of the staggered Arakawa D grid, in which the wind velocity components are evaluated at the respective cell faces and the rest of scalar quantities at the cell centres (Fig. 2). Note that this configuration is very convenient for solving the 3-D equation in a fractional manner because, when solving each one-dimensional case, wind velocities are already aligned with the boundaries of the corresponding one-dimensional cells. The previous versions of FALL3D used the classical LW central-difference scheme for solving the resulting one-dimensional ADS equations in Eq. (67), combined with a slope limiter to reduce numerical over-/undershooting near discontinuities. This resulted in a second-order accuracy except near sharp concentration gradients where, nonetheless, accuracy remained higher than in a first-order upwind method. The main advantage of using the LW scheme was its simplicity but, in contrast, it is well known that it introduces numerical dissipation that prevents accurate resolution of discontinuities, leading to over-diffusive results. In order to circumvent this drawback, FALL3D-8.0 uses instead a high-resolution Kurganov–Tadmor (KT) scheme that can be combined either with a fourth-order explicit Runge–Kutta or with a first-order Euler time-marching method. The latter option, even if less accurate, is still supported for computational efficiency reasons because it implies 4 times less solver calculations.

Consider the general one-dimensional advection-diffusion equation for a scalar variable c(x,t) on its conservative form:

$$\begin{array}{}\text{(68)}& {\displaystyle \frac{\partial c}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}F\left(c\right)={\displaystyle \frac{\partial }{\partial x}}G\left(c,{\displaystyle \frac{\partial c}{\partial x}}\right),\end{array}$$

where, in our particular case, F=c u is the advective flux and $G=k\partial c/\partial x$ is the diffusive flux (as already introduced, u(x,t) and k(x,t) are the velocity and diffusivity, respectively). Consider also a 1-D computational domain discretised as in Fig. 2, where c is computed at cell centres (“mass” points) and u is stored at staggered cell boundaries, which do not need to be equally spaced. The semi-discrete form of the KT scheme can be written at centre of each cell i in terms of fluxes at boundaries $i\pm \mathrm{1}/\mathrm{2}$ as (Kurganov and Tadmor, 2000):

$$\begin{array}{ll}{\displaystyle \frac{\partial c_i}{\partial t}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{\stackrel{\mathrm{‾}}{x}}_i}}\left[F_{i+\mathrm{1}/\mathrm{2}}^{*}-F_{i-\mathrm{1}/\mathrm{2}}^{*}\right]+{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{\stackrel{\mathrm{‾}}{x}}_i}}\left[G_{i+\mathrm{1}/\mathrm{2}}^{*}-G_{i-\mathrm{1}/\mathrm{2}}^{*}\right]\\ \text{(69)}& {\displaystyle }& {\displaystyle }=f(t,c),\end{array}$$

where $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{x}}_i$ is the ith cell width and

$$\begin{array}{ll}{\displaystyle }{F}_{i+\mathrm{1}/\mathrm{2}}^{*}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[F\left(c_{i+\mathrm{1}/\mathrm{2}}^{\text{r}}\right)+F\left(c_{i+\mathrm{1}/\mathrm{2}}^{\text{l}}\right)\right]\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_{i+\mathrm{1}/\mathrm{2}}\left(c_{i+\mathrm{1}/\mathrm{2}}^{\text{r}}-c_{i+\mathrm{1}/\mathrm{2}}^{\text{l}}\right)\\ {\displaystyle }{F}_{i-\mathrm{1}/\mathrm{2}}^{*}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[F\left(c_{i-\mathrm{1}/\mathrm{2}}^{\text{r}}\right)+F\left(c_{i-\mathrm{1}/\mathrm{2}}^{\text{l}}\right)\right]\\ \text{(70)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_{i-\mathrm{1}/\mathrm{2}}\left(c_{i-\mathrm{1}/\mathrm{2}}^{\text{r}}-c_{i-\mathrm{1}/\mathrm{2}}^{\text{l}}\right)\end{array}$$

$$\begin{array}{ll}{\displaystyle }{G}_{i+\mathrm{1}/\mathrm{2}}^{*}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[G\left(c_i,{\displaystyle \frac{c_{i+\mathrm{1}}-c_i}{\mathrm{\Delta }{x}_i}}\right)+G\left(c_{i+\mathrm{1}},{\displaystyle \frac{c_{i+\mathrm{1}}-c_i}{\mathrm{\Delta }{x}_i}}\right)\right]\\ {\displaystyle }=& {\displaystyle }G\left({\displaystyle \frac{c_{i+\mathrm{1}}-c_i}{\mathrm{\Delta }{x}_i}}\right)\\ {\displaystyle }{G}_{i-\mathrm{1}/\mathrm{2}}^{*}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[G\left(c_{i-\mathrm{1}},{\displaystyle \frac{c_i-c_{i-\mathrm{1}}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}}}\right)+G\left(c_i,{\displaystyle \frac{c_i-c_{i-\mathrm{1}}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}}}\right)\right]\\ \text{(71)}& {\displaystyle }=& {\displaystyle }G\left({\displaystyle \frac{c_i-c_{i-\mathrm{1}}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}}}\right).\end{array}$$

Note that, in the expression above, the last equality holds because our flux G depends only on the gradient of c, i.e. $G=G(\partial c/\partial x)$. In Eq. (70), $a_{i\pm \mathrm{1}/\mathrm{2}}$ is the maximum absolute value of the eigenvalue of the Jacobian matrix of F (in our particular case, it reduces to $a_{i\pm \mathrm{1}/\mathrm{2}}=\left|u_{i\pm \mathrm{1}/\mathrm{2}}\right|$), and c^r and c^l are, respectively, right and left values at a cell boundary, computed as

$$\begin{array}{ll}\text{(72)}& {\displaystyle }{c}_{i+\mathrm{1}/\mathrm{2}}^{\text{r}}& {\displaystyle }=c_{i+\mathrm{1}}-\mathrm{0.5}\mathit{\varphi }\left(r_{i+\mathrm{1}}\right)\left(c_{i+\mathrm{1}}-c_i\right){\displaystyle }{c}_{i-\mathrm{1}/\mathrm{2}}^{\text{r}}& {\displaystyle }=c_i-\mathrm{0.5}\mathit{\varphi }\left(r_i\right)\left(c_{i+\mathrm{1}}-c_i\right)\\ {\displaystyle }{c}_{i+\mathrm{1}/\mathrm{2}}^{\text{l}}& {\displaystyle }=c_i+\mathrm{0.5}\mathit{\varphi }\left(r_i\right)\left(c_i-c_{i-\mathrm{1}}\right)\\ {\displaystyle }{c}_{i-\mathrm{1}/\mathrm{2}}^{\text{l}}& {\displaystyle }=c_{i-\mathrm{1}}+\mathrm{0.5}\mathit{\varphi }\left(r_{i-\mathrm{1}}\right)\left(c_i-c_{i-\mathrm{1}}\right),\end{array}$$

where

$$\begin{array}{}\text{(73)}& r_i={\displaystyle \frac{c_i-c_{i-\mathrm{1}}}{c_{i+\mathrm{1}}-c_i}},\end{array}$$

and ϕ(r) is a flux limiter function (e.g. Sweby, 1984). Options available in FALL3D-8.0 are the well-known superbee ϕ_s(r) and minmod ϕ_m(r) (Roe, 1986):

$$\begin{array}{ll}{\displaystyle }{\mathit{\varphi }}_{\text{s}}\left(r\right)=& {\displaystyle }max(\mathrm{0},min(\mathrm{1},\mathrm{2}r),min(\mathrm{2},r\left)\right)\\ \text{(74)}& {\displaystyle }{\mathit{\varphi }}_{\text{m}}\left(r\right)=& {\displaystyle }max(\mathrm{0},min(\mathrm{1},r\left)\right).\end{array}$$

Note that c is needed at two extra mass points in order to evaluate $c_{i-\mathrm{1}/\mathrm{2}}^{\text{l}}$ and $c_{i+\mathrm{1}/\mathrm{2}}^{\text{r}}$ at the left/right cell boundaries, respectively. In other words, the “stencil” of the KT scheme needs two ghost nodes at the boundaries of the computational domain or a two-point halo for internal domains in the case of parallel domain decomposition. Note also that the solving strategy derives from a one-dimensional finite-volume formulation that, in our one-dimensional case, is also in practice equivalent to using linear finite elements.

Figure 3 Benchmark 1 results. Step-like discontinuity given by Eq. (80) moving right at constant velocity u=1 and without diffusion k=0 on a domain $x\in [-\mathrm{1},\mathrm{1}]$ with periodic boundary conditions. The solid black line shows the analytic solution after each periodic cycle, which coincides with the initial condition (t=0). Red dots show the KT + RK4 numerical solution after simulating 10 cycles (at t=20). Blue dots and blue circles show the LW + EU1 results after 1 cycle (t=2) and 10 cycles (t=20), respectively. Note that mass (area below curves) is conserved in all the cases, but the KT + RK4 scheme adds almost no numerical diffusion and preserves discontinuities. Results use 200 equally spaced grid cells.

Time marching from tⁿ to $t^{n+\mathrm{1}}=t^n+\mathrm{\Delta }t$ in Eq. (69) can be performed with the explicit first-order-in-time Euler method (EU1):

$$\begin{array}{}\text{(75)}& c^{n+\mathrm{1}}=c^n+\mathrm{\Delta }tf\left(t^n,c^n\right),\end{array}$$

or, alternatively, using the classical fourth-order Runge–Kutta method (RK4), in which

$$\begin{array}{}\text{(76)}& c^{n+\mathrm{1}}=c^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{6}}}(k_{\mathrm{1}}+\mathrm{2}{k}_{\mathrm{2}}+\mathrm{2}{k}_{\mathrm{3}}+k_{\mathrm{4}}),\end{array}$$

with

$$\begin{array}{ll}\text{(77)}& {\displaystyle }{k}_{\mathrm{1}}& {\displaystyle }=f\left(t^n,c^n\right){\displaystyle }{k}_{\mathrm{2}}& {\displaystyle }=f\left(t^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}},c^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{k}_{\mathrm{1}}\right)\\ {\displaystyle }{k}_{\mathrm{3}}& {\displaystyle }=f\left(t^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}},c^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{k}_{\mathrm{2}}\right)\\ \text{(78)}& {\displaystyle }{k}_{\mathrm{4}}& {\displaystyle }=f\left(t^n+\mathrm{\Delta }t,c^n+\mathrm{\Delta }t{k}_{\mathrm{3}}\right),\end{array}$$

where the function f(t,c) is given on the right-hand side of Eq. (69). In any case, the Courant–Friedrichs–Lewy (CFL) condition can be imposed to guarantee convergence in time integration (e.g. Hindmarsh et al., 1984) along each one-dimensional problem:

$$\begin{array}{}\text{(79)}& \mathrm{\Delta }t\le min\left({\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{2}{k}_x}{\mathrm{\Delta }{x}^{\mathrm{2}}}+\frac{u}{\mathrm{\Delta }y}}},{\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{2}{k}_y}{\mathrm{\Delta }{y}^{\mathrm{2}}}+\frac{v}{\mathrm{\Delta }y}}},{\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{2}{k}_z}{\mathrm{\Delta }{z}^{\mathrm{2}}}+\frac{w}{\mathrm{\Delta }z}}}\right),\end{array}$$

multiplied by a user-defined safety factor. This factor should theoretically be lower than 1 in fully explicit cases but, given the semi-implicit nature of the splitting algorithm, slightly larger values can also yield to stability.

4.3 Algorithm benchmarks

Three benchmark cases serve us to illustrate the gains of the KT + RK4 numerical scheme with respect the former LW + EU1 implemented in the previous versions of the code.

Example 1 considers the pure advection (k=0) of a step-like discontinuity. Consider a 1-D domain $x\in [-\mathrm{1},\mathrm{1}]$ with an initial concentration of

$$\begin{array}{}\text{(80)}& {\displaystyle }\left\{\begin{array}{l}c(t=\mathrm{0})=\mathrm{1}\left|x\right|\le \mathrm{0.5}\\ c(t=\mathrm{0})=\mathrm{0}\left|x\right|>\mathrm{0.5}\end{array}\right.\end{array}$$

that is advected by a uniform velocity field u=1. Periodic conditions are imposed at the boundaries, so that “mass” leaving the computational domain at x=1 is re-injected at $x=-\mathrm{1}$. As a result, the initial condition is periodically recovered after each cycle with a period t=2. Results are shown in Fig. 3. Note how, as opposed to LW, the KT scheme adds almost no numerical diffusion and preserves discontinuities. In addition, because of its higher order in time, the KT + RK4 scheme preserves the solution, whereas for LW + EU1 accuracy deteriorates with time (compare the accuracy of the LW + EU1 solutions after 1 and 10 cycles in Fig. 3).

Figure 4 Benchmark 2 results. Steady-state solutions for different Péclet numbers (Pe) on a domain $x\in [-\mathrm{1},\mathrm{1}]$ with boundary conditions c=0 at $x=-\mathrm{1}$ and c=1 at x=1. Black lines show the respective analytic solutions given by Eq. (81). Red and blue dots show, respectively, the KT + RK4 and LW + EU1 numerical solutions. The case of Pe=0 (pure diffusion problem) has a linear solution and is exactly matched by both schemes. For the rest of the cases, KT + RK4 outperforms LW + EU1. Results use 200 equally spaced grid cells.

Example 2 considers the classical 1-D advection–diffusion problem for the onset of numerical instability in a domain $x\in [-\mathrm{1},\mathrm{1}]$ subject to the boundary conditions c=0 at $x=-\mathrm{1}$ and c=1 at x=1. The problem has a steady-state analytic solution given by

$$\begin{array}{}\text{(81)}& c(t\to \mathrm{\infty })=\left(e^{\mathit{Pe}(x+\mathrm{1})}-\mathrm{1}\right)/\left(e^{\mathrm{2}\mathit{Pe}}-\mathrm{1}\right),\end{array}$$

where $\mathit{Pe}=ul/\mathrm{2}k=u/k$ is the Péclet number. Figure 4 shows the steady-state solutions for different Péclet numbers, illustrating also how KT + RK4 outperforms LW + EU1 as the Péclet number increases.

Example 3 considers a case with pure advection (k=0) on a 2-D domain $(x,y)\in [-\mathrm{1},\mathrm{1}]\times [-\mathrm{1},\mathrm{1}]$ with an initial condition at t=0 given by a conic concentration distribution with a unit (c=1) peak concentration centred at $(x_{\text{c}},y_{\text{c}})=(\mathrm{0},\mathrm{0.695})$ and having a radius r=0.1 (Fig. 5a). The cone is advected by a rotating clock-wise velocity field centred at (0,0):

$$\begin{array}{ll}{\displaystyle }u& {\displaystyle }=\mathrm{\Omega }y\\ \text{(82)}& {\displaystyle }v& {\displaystyle }=-\mathrm{\Omega }x,\end{array}$$

with an angular velocity Ω=π (in Rad s⁻¹), so that each cycle is repeated with a period of t=2. Concentration profiles along two transects A (y=y_c) and B (x=0) after two cycles (i.e. at t=4) are shown in Fig. 5b and c, respectively. Note again the substantial improvement in the KT + RK4 scheme, which performs well even in this numerically challenging test.

Figure 5 Benchmark 3 results. (a) A cone centred at $(x_{\text{c}},y_{\text{c}})=(\mathrm{0},\mathrm{0.695})$ with radius r=0.1 (shaded circle) is advected by a rotating velocity field with angular velocity Ω=π in a computational domain $(x,y)\in [-\mathrm{1},\mathrm{1}]\times [-\mathrm{1},\mathrm{1}]$. Plots in panels (b) and (c) show concentration profiles after two cycles along lines A (y=y_c) and B (x=0). The solid black line shows the analytic solution after each cycle. Red and blue dots show the KT + RK4 and LW + EU1 numerical solutions, respectively. Results use 200 equally spaced cells along each direction.

Table 11 Summary of the FALL3D-8.0 model tasks.

¹ Model input file (same for all tasks). ² Number of processors along each spatial dimension in domain decomposition. If not given, the execution is serial.

Figure 6 Pre-process and execution workflow and associated model tasks. Task All runs all tasks in a single (parallel) execution. In this case, the write/read of the intermediate files can be omitted. All tasks share a unique model input file (name.inp).

5 Model execution workflow

In FALL3D-8.0, the pre-process auxiliary programmes have been parallelised and embedded in the code, so that a single executable file exists for all the pre-process steps and execution workflow (see Fig. 6). These formerly independent programmes can still be run individually as model tasks (specified by a programme call argument) or, alternatively, concatenated with the model in a single execution. In the first case, pre-process tasks generate output files that are later given as inputs to the model task. This is similar to what occurred in the previous version (v7.x) but with the difference of a parallel pre-process. In contrast, the second option does not require intermediate file writing/reading and therefore saves disk space and overall computing time. In any case, all tasks share a unique model input file and generate their own log file to track execution and report eventual warnings and errors. Possible task options are summarised in Table 11 and include the following:

1. Task SetTgsd can generate Gaussian and bi-Gaussian particle grain size distributions in Φ (log-normal in diameter) or, alternatively, Weibull and bi-Weibull distributions (Costa et al., 2016a), and assumes a linear variation of density and particle shape factor between two specified cut-offs. For other kinds of grain size distributions, the user must provide a total grain size distribution file (name.tgsd).

2. Task SetDbs interpolates all the required meteorological data from the original grid of the driving meteorological models to the computational domain. Table 12 summarises the different meteorological drivers available in FALL3D-8.0. Global datasets include ERA5, the new reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hersbach and Dee, 2016), the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) and Global Data Assimilation System (GDAS) final analysis. Regional models supported include the Advanced Research Weather Research and Forecasting (WRF) (ARW) core of the WRF model (Skamarock et al., 2008), the mesoscale models HARMONIE-AROME (Bengtsson et al., 2017) and the COSMO-LAMI, run by the Italian Regional Environmental Protection Agency (ARPA). Note that some former dataset options have been deprecated in v8.x. The FALL3D-8.0 distribution package provides a set of utilities to download and pre-process meteorological data for the SetDbs model task. Python scripts are provided to download and crop the variables required by the model. ERA5 can be obtained on either model levels (137 vertical levels) or pressure levels (37 vertical levels) via the Climate Data Store (CDS) infrastructure. GFS datasets can be accessed using the online archive of real-time weather model output from the National Operational Model Archive and Distribution System (NOMADS) (Rutledge et al., 2006). The NCEP FNL (final) operational global analysis and forecast data from the Global Data Assimilation System (GDAS) can be accessed using the Open-source Project for a Network Data Access Protocol (OPeNDAP) through the Thematic Real-time Environmental Distributed Data Services (THREDDS) Data Server (TDS) offered by the National Center for Atmospheric Research (NCAR) Research Data Archive (RDA).

3. Task SetSrc generates different emission source terms (see Sect. 3.2.3 and Table 4), including the PLUME option based on the FPLUME-1.0 model (Folch et al., 2016). If necessary, it also performs a priori tephra particle aggregation (Sect. 3.2.6) and a TGSD cut-off in order to select the effective bins considered in the atmospheric transport.

4. Task FALL3D runs the FALL3D-8.0 model itself.

5. Task All runs all previous tasks consecutively as a single parallel execution.

Table 12 Summary of the meteorological drivers available in FALL3D-8.0.

6 Parallelisation and performance

Parallelisation in FALL3D-8.0 considers a 3-D domain decomposition, with freedom for user to choose the number of processors along each direction. In contrast, previous versions considered two different levels of parallelisation: one on particle bins and another on domain but only along the vertical dimension. Parallelisation on bins was convenient in v7.x because no interaction among bins existed, but such a form of trivial parallelism has been deprecated given that it would now yield to unnecessary communication penalties. Note also that the full 3-D domain decomposition allows solving on much larger grid sizes before reaching hardware memory limits.

In terms of runtime performance, it is important to recall that the splitting algorithm combined with the RK4 time marching implies solving a series of 1-D equations four times along each spatial dimension, contrasting with a single solution for the EU1 case. In other words, each time integration step of Eq. (67) using the RK4 scheme along X₁ implies solving n_y×n_z one-dimensional problems four times, the solution along X₂ implies solving n_x×n_z problems four times and so on. In order to minimise the number of communications between neighbouring processors, it is more convenient to compute first each of the four partial increments in Eq. (78) for all the 1-D problems on a given dimension. This optimises swapping communications among domain partitions because a single swapping request (i.e. one MPI send/receive call) communicates all the data necessary to compute the next RK4 increment for the whole mesh. In contrast, if all partial increments in Eq. (78) were computed for each 1-D problem, it would require as many swapping requests as 1-D problems involved. Clearly, this solving approach is convenient but it presents two drawbacks. Firstly, it obviously increases the amount of memory required and, secondly, it poses an issue on memory data access, with potential increase of cache memory misses given the more frequent access to data. Moreover, when one solves for the dimensions stored on array lowest axes, memory data access is scattered and therefore less efficient. This issue has been addressed in FALL3D-8.0 by rearranging the components of the velocity and diffusion arrays, with the dimension to solve being stored always in the fastest axis, and also by transposing the concentration array that has to be updated at the first partial increment. Such a memory management strategy allows exploiting always contiguous cache memory positions on any spatial dimension.

Figure 7 Simulation results for the 2011 Cordón Caulle eruption. The model was configured as shown in Table 13 and solved on the $\mathrm{500}\times \mathrm{500}\times \mathrm{100}$ grid cells used in Sect. 6 for the scalability analysis. (a) PM₁₀ ash cloud column mass contours (in kg m⁻²) on 5 June 2011 at 16:00 UTC, just after satellite data insertion. (b) Same on 7 June 2011 at 01:00 UTC after model evolution assuming a constant column height of 10 km a.v.l.

Figure 8 Strong scaling results on MareNostrum 4 supercomputer for the Cordón Caulle simulation using a $\mathrm{500}\times \mathrm{500}\times \mathrm{100}$ cell grid. (a) Scaling curves up to 2048 processors for the RK4 (red) and EU1 (blue) time integration schemes. Results using code v7.3.4 are also shown with up to 64 processors only (pink). Ideal strong scaling behaviour is indicated by the solid black line. Note that the scaling curves refer to the total computing time and therefore include I/O operations, not just computing time. (b) Parallel efficiency (Eq. 83) depending on the number of computation units.

Figure 9 Computing time ratio between code v7.3.4 and v8 depending on the number of processors. Results for EU1 (blue bars) and RK4 (red bars) including I/O operations.

Costa et al. (2016a)Byun and Schere (2006)

Table 13 FALL3D-8.0 model configuration for the Cordón Caulle simulation example shown in Sect. 6.

All the aspects discussed above have improved substantially the performance and scalability of the code. As an example of strong scaling (i.e. time to solution for a fixed problem size), let us consider a real-case ash dispersal simulation from the 2011 Cordón Caulle eruption (e.g. Collini et al., 2013). The model was configured as shown in Table 13 and solved on a typical grid size of $\mathrm{500}\times \mathrm{500}\times \mathrm{100}$ computational cells (horizontal model resolution of $\approx \mathrm{0.03}{}^{\circ }$). For illustrative purposes, Fig. 7 shows model snapshots of ash cloud column mass at two different times. Figure 8a shows strong scaling results (speedup) up to 2048 processors obtained on the MareNostrum 4 supercomputer, composed by general-purpose nodes with 48 Intel Xeon Platinum processors interconnected by a 100 GB Intel Omni-Path full-fat tree. For comparison, this figure shows also the strong scaling curve obtained with latest code (v7.3.4), in this case limited to 64 processors (larger values were not possible on this grid using v7.3.4 given the former one-dimensional domain decomposition along z). Figure 8b plots the corresponding parallel efficiencies, defined as

$$\begin{array}{}\text{(83)}& PE=\mathrm{100}\times {\displaystyle \frac{t_{\mathrm{1}}}{N{t}_{\text{N}}}},\end{array}$$

where t₁ is the total computing time with one processing unit, and t_N is the time with N processing units. Note that ideal strong scaling implies a parallel efficiency of 100 %. Clearly, v8.x improves notably the scalability of the code, with a trend close to that of perfect scaling up to ≈ 100 processors (parallel efficiency ≥ 90 %). Depending on the time integration scheme, values of parallel efficiency above 50 % are obtained with up to 512 and 1024 processors for EU1 and RK4, respectively. This is in striking contrast with v7.3.4, for which parallel efficiency already drops below 50 % with only 16 processors (Fig. 8b).

Much better v8.x code performance is also observed not only in terms of code scalability but also in terms of total computing time. Figure 9 shows the computing time ratio (total elapsed time) between v7.3.4 and v8 depending on the number of processors and the time integration scheme. The EU1 case (i.e. equal order of accuracy in time to that in v7.3.4), shows a similar serial performance, only with an insignificant computation overhead probably due to algorithmic differences. However, this overhead is rapidly balanced with only four processors, and v8.x with EU1 is already up to 4 times faster than v7.3.4 with 64 processors. The RK4 case (i.e. fourth order of accuracy in time) is around 4 times slower than v7.3.4 for the serial execution. This is a logical result given that the RK4 scheme needs to iterate 4 times more to solve in time. However, this penalty is also rapidly balanced as the number of processors increases, and RK4 already outperforms v7.3.4 in absolute terms with only a few tens of processors. In summary, v8.x clearly shows a much better scalability and performance (absolute computing time to solution for a given problem and computational resources) than v7.3.4. The companion paper (Prata et al., 2020) contains a detailed model validation and shows that this is also true in terms of model accuracy.

7 Conclusions

After 15+ years, the atmospheric transport model FALL3D has been completely rewritten and modernised to overcome legacy constraints in the former release (v7.x) that precluded the introduction of new functionalities and seriously limited the scalability and performance of the code on hundreds or thousands of processors. With this, FALL3D-8.0 can be considered as a baseline for the successive optimisations and preparation for exascale computing. However, as detailed in the paper, v8.x already contains remarkable improvements and updates on model physics, numerics and performance. In particular, the code has been prepared to deal also with particles different from tephra, aerosols and radionuclides, includes new coordinate mapping options, a more efficient and less diffusive solving algorithm (KT) that can be combined with a high-order-in-time solver (RK4) and a better memory management and parallelisation strategy based on a full 3-D domain decomposition. Strong scaling results have shown perfect scaling with a few hundred processors and a parallel efficiency above 50 % with 1024 processors. This remarkable improvement is also true in terms of performance (total computing time), with v8.x outperforming the previous release by a factor of 4 with only 64 processors. Further expected improvements in the preparation towards exascale include memory optimisation, introduction of thread parallelism (OpenMP), code vectorisation, porting to accelerators (GPUs), performance portability, load balance, asynchronous I/O and preparation for emerging heterogeneous architectures (exascale hardware prototypes).