New model for dispersion of volcanic ash and dust in the troposphere

Jonas Eliasson^1,2

¹Earthquake Research Institute, University of Iceland, Selfoss, Iceland.
²Disaster Prevention Research Institute, Kyoto University, Kyoto, Japan.

DOI: 10.4236/ijg.2020.118029                

Abstract

Dispersion of volcanic ash and dust is traditionally modeled as advection and Gaussian diffusion. This is the tradition in treating smoke stack plumes. About 100 meters above earth the velocity profile may disintegrate, diffusion coefficients become rather unpredictable and stratified flow occur. It is suggested that gravitational flattening may be the main cause of dispersion in dust plumes above the turbulent boundary layer. A dust plume in between two layers of small temperature difference has a certain carrying capacity of dust. The corresponding mass loading can be estimated from the temperature difference between the layers above and beneath the plume. Such dust plumes will be forced to jettison a load they may have in excess of this carrying capacity; this may be seen as streak fallout from the plume. In the same time, the plume will be subjected to gravitational flattening to the sides, in addition to any diffusion if there is any. The plume width resulting from the flattening may be estimated from the temperature difference. This can explain the behavior of plumes like the plume from the Eyjafjallaj?kull 2010 in absence of diffusion. In the long run diffusion and gravitational flattening will cause different developments of the plume width. Gravitational flattening and streak fallouts are important elements from plume physics not included in most plume models. It is concluded that modelling dust plumes with diffusion and ordinary fallout only; can cause serious errors in the model, the simulated plumes will become too big. To avoid them, the new model should be included in dust models in the same manner as the turbulent diffusion, i.e. as a sub grid model. Then, the plume model only needs to include horizontal turbulent diffusion of the same order of magnitude as the vertical one.

Keywords

Dust Plume, Gravitational Dispersion, Load Capacity, Streak Fallout

1. Introduction

The scope of this article is the science of the advection and dispersion of buoyant plumes and the physics governing their migration and the modeling of dust plumes. The traditional technique of mathematical modeling of volcanic plumes uses the theory of atmospheric dispersion utilized in ambient air quality research of air pollutants from smoke stacks. Smoke plumes are almost all in the PBL (Planetary Boundary Layer) where the diffusion is governed by the Gaussian air pollutant dispersion [1] equation and the mixing force is the shear turbulence created by the vertical gradient of the horizontal wind velocity profile close to the surface as originally proposed by the early theories on the structure of turbulence, now accessible in textbooks on Fluid Dynamics.

Today, valuable information on dust plumes can be acquired from satellite images. In [2], Table 2 shows that horizontal dispersion coefficient, or K value, of 1300 - 1400 m²/s is needed to explain the sideways expansion of the plume from Mt. Etna. So high K values are doubtful here on Earth, but something similar has been suggested from observations of dust storms on the planet Mars [3].

Diffusion in dust plumes can be studied using traditional aerosol measurement techniques [4]. This technology can be useful if installed in airplanes. The articles [5] [6] [7] demonstrates how airborne measurement techniques were used on the famous Eyjafjallajökull plume in 2010. Various other observations from Europe and Japan point clearly to, that dispersion of dust plumes cannot be explained by turbulent diffusion only.

In this paper, it is suggested that gravitational flattening of dust clouds is responsible for the dispersion rather than diffusion, and the gravitational flow of non-buoyant parts out of the cloud—streak fallout—is what controls the mass flux in dust plumes in the long run, rather than ordinary fallout or what is blown up in the air at the source. This is further explained in Chapter 3, and a new formula developed for the maximum carrying capacity of a neutrally buoyant plume.

Research has established the existence of streak fallouts. They are vertical density currents carrying large quantities of dust to the ground. They are clearly visible in most pictures taken of weak volcanic plumes close to the source (wind-bent plumes). The reason for this is the large mass carried up by the heavily buoyant vertical hot plume at the source. It makes the density so large that it exceeds the carrying capacity of a neutrally buoyant plume in horizontal flight. The streak fallouts will quickly jettison the load exceeding the carrying capacity. Thus, it may be considered unimportant for big neutrally buoyant plumes, how large the mass flux at the source really is.

The major disruption of air transport in Europe during the Eyjafjallajökull volcanic eruption in 2010, caused a huge economic loss to the world, US$ 4 - 6 billion according to various estimates, mostly incurred by private firms. If this figure is taken at face value, this eruption is the costliest in history. However, it was a small eruption and local damage was almost nil compared to the total economic loss of the aviation industry. They suffered the greatest loss (about US$ 1.7 billion according to various estimates) as a consequence of a much-criticized flight ban, that caused thousands of cancellations.

The information used by the ANSP (Aeronautical Service Provider) to direct the air traffic during the Eyjafjallajökull volcanic eruption in 2010 were simulations that provided maps of the ash cloud. They showed too big clouds, [8] Figure 1, and the question is, do we have improved models today, 10 years on. An inspection of the most recent publications using the same model (NAME model), [9], show that neither gravitational dispersion, maximum load or streak fallouts are included in the new versions of the model. Nor are they included in [10], another new publication where an extended version of the WRF model is used to simulate the Holuhraun eruption in Iceland 2014.

The distance from Iceland to Europe is 1000 - 2000 km. Ash plumes cannot travel this distance without being neutrally buoyant. Volcanic plumes have a very small content of ash measured in terms of volume ratio; hence, it is not easy to see how the plumes can get that far and still be visible from the ground. But then it must be considered that the aviation authorities used maps of simulated ash clouds, not direct observations.

Figure 1. Cross-section in a volcanic plume drifting in a density stratification in the atmosphere. The y and z are the symmetry axis of the plume. (a) is the pressure distribution in the center (y = 0) of the plume, p₀. (b) is p₀ overlaid on the pressure at the ends p_L, in x = L and x = −L. (c) The horizontal arrows denote the direction of the pressure gradient and the velocity of the gravitational spreading V. In x = L and x = −L, V = V_L.

2. Conditions in Different Altitudes

2.1. An Example of Temperature and Wind Close to the Earth’s Surface

The most effective dispersion is advection by wind. Close to the surface the logarithmic velocity profile is very popular

$\left(U-U_0\right)/U^{*}=2.5\ln \left(y/y_0\right)$ (1)

U is the wind speed, suffix 0 refers to a known elevation and U^* is the shear velocity defining the surface friction. To know it we must have another measurement at another elevation, let’s say higher up than elevation 0, which is commonly 2 or 10 meters above ground level.

Estimates of U^* are difficult. In rugged terrain, landscape generated turbulence overtakes the skin friction effect and measurements have to be made higher up to find a stable U^* value. Over the oceans there is no landscape except the waves, in this conditions the relation U₁₀/U^* = 15 is a good example [11] and we have $U=U_{10}\left(1+0.167\ln \left(y/10\right)\right)$ .

In the elevation 80 - 100 meters this means about 35 % higher wind velocity than U₁₀. This is the original plausible assumption of turbulence theory, a constant shear velocity U^* from 10 to 100 meters. This leads to an estimate of an eddy viscosity, ε, of 2.5 - 25 m²/s in the elevation range y = 10 - 100 m for U₁₀ = 10 m/s. If this is extended up to 10 kilometers the velocity will be 115% higher and ε will become 2500 m²/s. ε is comparable with K, but there is no evidence that supports the extension of the velocity profile from 100 m up to 10000 m.

The logarithmic wind velocity profile can last about 100 m up, but it becomes somewhat unpredictable above that. Wind turbines reach to this height, and estimating the wind profile for designing them may become difficult [12].

Smoke stacks are mostly below 100 meters, to estimate the diffusion of pollution from them is quite a complicated process based on Fickian diffusion theory mathematics [13]. In a stationary wind U in the x direction, the concentration of a polluting agent C, from a smoke stack H meters high, is usually estimated by:

$\begin{array}{c}C\left(x,y,z\right)=C_0\left(x\right)\exp \left(-y^2/2{\sigma }_y^2\right)(\exp \left(-\left(z^2-H\right)/2{\sigma }_z^2\right)\\ +\exp \left(-\left(z^2+H\right)/2{\sigma }_z^2\right))\end{array}$ (2)

where y and z are the horizontal and vertical coordinates perpendicular to the direction of the advection (x-axis ). C₀ is the centerline concentration. The σ_yσ_z denote the standard deviations of the normal distribution; they vary with x but not with time. If Q is the steady pollution flux in the plume we have by integrating (2) in y and z.

$Q/U=C_0\left(x\right)2\pi {\sigma }_y{\sigma }_z$ (3)

It shows that $\sqrt{2\pi }\sigma$ denotes the effective width and height of the plume with respect to C₀. Associated with σ is the coefficient of dispersion K. It is really a vector, but for demonstration purposes it is assumed a constant $K=K_x=K_y=K_z$ ,

${\sigma }^2=2Kt$ (4)

For times longer than T_T, the integral time scale of the turbulence, we have

${\sigma }_x^2=2E\left[{u^{\prime }}^2\right]T_{T}t\Rightarrow K_x=E\left[{u^{\prime }}^2\right]T_T$ (5)

And similar for the other directions of the flow. The $E\left[{u^{\prime }}^2\right]$ stands for the expected value (long time average) or variance, of the $u^{\prime }$ (x-axis) component of the turbulent velocity fluctuation $\left(u^{\prime },v^{\prime },w^{\prime }\right)$ . In stationary flow with $U\gg 0$ this component of dispersion is often disregarded because of the much stronger advection; this is done in Equation (2). In stationary models where turbulence does not exist except in a subgrid model this means that t = x/U in Equations (4) and (5). When the mean velocities in y and z are V = W = 0, the non-zero velocity fluctuations $u^{\prime },v^{\prime }$ and $w^{\prime }$ are all related through the continuity equation, so a $v^{\prime }\gg w^{\prime }$ is a little unrealistic. High up one would expect isotropic turbulence with $E\left[{u^{\prime }}^2\right]=E\left[{v^{\prime }}^2\right]=E\left[{w^{\prime }}^2\right]$ but $u^{\prime },v^{\prime }$ and $w^{\prime }$ uncorrelated in time. In shear turbulence with velocity profile like Equation (1) we have

$U^{\ast 2}=E\left[u^{\prime }{w}^{\prime }\right]=E\left[{u^{\prime }}^2\right]c_{uw2}\approx 0.01U_{10}^2\Rightarrow K_x=0.0044U_{10}^2{T}_T/c_{uw}$ (6)

For the correlation coefficient we have c_uv < 1. Equation (6) may also be related to the eddy viscosity that can be deducted from Equation (1), but only in the range where the shear velocity is reasonably constant. Equation (6) shows that K may be estimated through measurements of the velocity profile, but the more correct way is to record the time history of $u^{\prime }$ and use Equation (5).

The logarithmic velocity profile can break down for several reasons, in the temperate zone its thickness may be around 100 meters or less. The reason may be non-isenthalphic temperature stratification or another type of density stratification.

2.2. Stratified Flow

In [14], the surface layer closest to the earth is named the inertial layer and above that is the Ekman layer. Here internal waves in stratified flow may occur with interfacial friction and mixing and occasionally entrainment. In stratified flows the temperature gradient fluctuates. In tables for standard atmosphere, a temperature lapse rate of 6.5˚C/km is assumed. This is about 30 % lower than the adiabatic lapse rate, giving rise to convection and the buildup of layers of slightly different temperature that may be temporarily stable.

Internal waves in a stable stratification cause internal friction and creation of turbulent energy, that again causes mixing of the two layers. The mixing can be very slow however. If a hotter layer is flowing on a colder with a higher velocity U, the stability of the interface is controlled by the densimetric Froude number defined as $F_{R\Delta }={\left(\Delta g\rho H\right)}^{1/2}$ , where

$\Delta =2\left({\rho }_1-{\rho }_2\right)/\left({\rho }_1+{\rho }_2\right)$ (7)

is the relative density ratio of the layers, g the acceleration of gravity and H the depth of the flow. For high F_RΔ the stability is low and large internal waves can develop with eventual breaking and intensive mixing. For lower F_RΔ values entrainment of the slower fluid into the faster fluid will occur, but this can be seriously hampered if the Richardson number, especially the bulk flux Richardsons number [15], is high. The reason is that mixing of denser fluid into a less dense fluid creates potential energy that diminishes the production of turbulent energy. Temperature differences between layers thus means lower interfacial shear stress and slower mixing of the two fluids and dispersion of contaminants.

The mathematics behind both the flow and the mixing processes are very difficult. They are known, [16] but the data necessary to use them, very accurate temperature and velocity profiles, are in general unavailable.

We are thus left with qualitative analysis, lacking data for quantitative analysis. Some conclusions can be drawn however, e.g. if we have three layers in the above example a high layer above, a deep layer beneath and a layer of thickness H in the middle, H = 100 m say. If the difference in average velocity is U between the lowest and the highest layers, and the difference in potential temperature is a few degrees ΔT, the densimetric Froude number will, by virtue of the equation of state for clean air, be $F_{R\Delta }={\left(U/\left(Hg\Delta T/T\right)\right)}^{1/2}\approx 0.5{\left(U/\Delta T\right)}^{1/2}$ . When ΔT is only a few degrees, U must be small to keep F_RΔ under 1 and the flow stable, otherwise the 100 m thick middle layer starts to undulate and mix into the other two layers. In addition, the average velocity in the middle layer, relative to the others, will be only U/2 so production of turbulent energy is low, mixing and dispersion slow, though it may be significant around the level z = H/2. But in the same time, we can have K ≈ 0 in the mixing levels z = 0 and z = H, if the flow is stable and the density gradients are stable, then they slow down the mixing, this is an effect of the bulk flux Richardson number.

Where such stratified flow exists, it creates conditions that can carry contamination in the atmosphere a long way inside the middle layer without it being mixed up or down as it would be in the inertial layer below. One can assume, that it is something like this, that is the idea behind the assumption K_z = 0 in the existing dust models. But then gravitational flattening will be active as long as there exist any horizontal pressure gradients, and that will be the case, if H is not constant but varies from place to place which it does if a more or less round plume invades the middle layer. Also, when the density inside the middle layer is constant. To take a simple example, an oil drop on a water surface flattens out to the sides and ends up as a flat oil slick, without any turbulent mixing involved.

2.3. Mixing Conditions in the Troposphere above the Logarithmic Vellocity Profile

In the standard atmosphere the average lapse rate of 6.5˚C/km is normally assumed to reach to 12 km. Taylors hypothesis, the frozen turbulence hypothesis, is widely accepted. It means that the spatial characteristics of the turbulence are unchanged in the direction of the wind as long as the wind does not change. This gives local measurements regional validity, i.e. we can assume K to be a constant in large regions.

In stably stratified flow Taylors hypothesis can fail. After a thorough analysis it is concluded [16] that turbulence in stable conditions is not in equilibrium with the nonturbulent (averaged) flow. The turbulence has a non-stationary structure that sometimes is wavy with more complex signatures. In other words, we are left with the rather pessimistic view, that if K is not zero, the possibilities to find it using relations like Equation (6) are practically nil.

There are certain possibilities to obtain information on K using measurements to estimate it from temperature and velocity data via the Richardsons number [17]. Their results show K that is always below 100 m²/s up to 20 km, and only occasionally above 10. In [18] radar measurements show similar results, average K's in the lowest 20 km are even lower, but the elevations 84 - 86 km show K-values around 200 m²/s. But this altitude is too high up for having practical interest for the dispersion of dust.

The dispersion characteristics will thus depend heavily on the stability. Nothing new in that of course, but stability depends heavily upon the density stratification, and density stratification can change other components of the flow than just the turbulence.

3. Dust Particles in Stably Stratified Flow

3.1. Settling Velocity

Most dust particles consist of silica-rich minerals with density of around 2500 kg/m³. They fall through still air with a terminal settling velocity, V_S, that can be related to the effective diameter of the grains. That can be the diameter of a sphere with a mass equal to the grain, or an equal V_S. The V_S depends on the flow resistance, which is laminar for small grains, or linearly dependent on V_S in still air. For heavier grains, or higher velocities, flow resistance depends on $V_S^2$ . In the laminar range the problem is simple, due to mass flow continuity we will have a flow of grains, or a fallout from the dust cloud, equal to

$\text{d}M/\text{d}t=C\left(z\right)V_S$ (8)

where M is the mass loading above the latitude level z. But turbulence changes all that. If there is a quadradic term in the resistance and we have a velocity fluctuation w' in a vertical direction, the averaged flow resistance of a particle in a steady fall will depend on

$E\left[{\left(V_S+w^{\prime }\right)}^2\right]=V_S^2+{w^{\prime }}^2$ (9)

And we can easily have w' > V_S. This can slow the fallout tremendously, especially the fallout of big grains. Accounting for velocity variations of the grain and the other fluctuation components u' and v', makes Equation (9) so complex we cannot find the V_S without knowing the full spectrum of the turbulence in details.

3.2. Fallout

But whatever the foregoing result is, we still have a V_S for each grain size fraction and the concentration in each vertical is determined by the differential equation for the particle flow q(z) inside the plume,

$q\left(z\right)=-K\text{d}C/\text{d}z-C{V}_S;\text{d}C/\text{d}t=\text{d}q/\text{d}z$ (10)

If we can assume that the density differences are large enough to create stable conditions in the interfaces of the plume and the clean air above and below the interesting possibility arises that we can have K ≈ 0 in the top and the bottom of the plume (z = 0 or H). Several turbulence models are available for the vertical distribution of K in other respect, they have the largest K in the middle, and in between we can have K constant equal to the largest K (isotropic case), changing linearly (constant shear turbulence) or parabolically (linear shear stress variation) down to the zero at the boundary. Then we can solve (10) with a constant C as an initial condition, C = 0 for z = H and q(0) = −C∙V_S for z = 0 as boundary conditions. But we need assumptions like that to reach a solution of the fallout problem. Without them the boundary conditions, especially the bottom condition q(0) = −C∙V_S, one of the main principles in all dust cloud modelling, are not valid.

If K is very big dC/dz will be very small. If we have as initial condition an average concentration C₀ at t = 0 in the vertical, we will have a mass loading of dust M_C = C₀H. The fallout according to Equation (10) will quickly create a concentration profile, the highest value C_B the bottom z = 0, monotonously decreasing down to C = 0 at the top. The fallout q(0) = −C_BV_S will deplete the mass loading according to the value of C_B at any time. From the time C_B starts depleting we may have for a particle fraction with the terminal fall velocity V_S

$C_B=C_0{\text{e}}^{-V_{S}t/H}$ (11)

Equation (11) supposes a concentration profile $C\left(z\right)/C_B=f\left(z/H\right)$ where the f is a near-parabolic curve. This is a very slow process for long times t. In contrast to this is the situation K = 0, or zero turbulence everywhere. Then the whole mass loading of the particular particle fraction will fall out of the cloud like a stone in the time t = H/V_S with C_B = C₀ until the entire particle fraction is cleared out of the plume. As an example, the 20 μ fraction will be cleaned out of a 1 km thick cloud in 5 days, but if there is turbulence there will still be about 35% of the original mass loading left.

Equation (11) explains how big grains survive the transport over so long distances that they should have disappeared from the cloud. Examples of this are coarse grains from Iceland found in deposits in Europe [19].

It can be inferred from Equation (10) and Equation (11) that the lower half of the plume will always be heavier than the upper half. But quantitative estimations of this effect may lead to misleading results. The real V_S values are unpredictable in turbulent air inside the cloud.

It can also be inferred, that the ordinary fallout rule used in plume models, q(0) = −C₀V_S, assumes K = 0 in the bottom of the plume, but in the same time a K value up in the plume that is high enough to facilitate sideways dispersion due to turbulent mixing. Such K values are of the order of magnitude several hundreds or thousands of m²/s. This is rather inconsistent, as already pointed out in section 2.1. K = 0 in the bottom should mean K values inside the plume of the order shown in [17] and similar papers. These values are commonly in the order of magnitude K = 10 - 20 m²/s.

4. Dust Plume in a Stable Inversion

4.1. The Stable Plume: Maximum Dust Loading

Consider a plume within a density inversion, it will have the upper part of the plume in the upper layer (ρ = ρ₂) and a lower part immersed in the lower layer below (ρ = ρ₁). If neutrally stable the plume will sail on without floating up or sinking down. This means that the mass loading in any vertical is in equilibrium with the surrounding clean air. If it is not, the plume is not neutrally buoyant and will either float up or sink down.

Due to the large difference of the densities of the dust and the clean air, the plume’s density is

$\rho ={\rho }_0\left(1+C\right)$ (12)

in any point, ρ₀ being the clean air density for (C = 0). We assume for simplicity that the K is high enough to keep C fairly constant within the plumes vertical so we can disregard the small density differences within the plume. The air above the inversion must be ΔT ˚K hotter than the air underneath to keep the inversion stable. A mass loading equal to the total buoyancy gives

${\rho }_0\left(1+C\right)=\left(1-\alpha \right){\rho }_1+\alpha {\rho }_2$ (13)

Here, we have 0 < α < 1 as a measure of the vertical position of the plume relative to the inversion plane. In the case C = 0 and ${\rho }_0={\rho }_2$ or ${\rho }_0={\rho }_0\left(1+C\right)$ we have the stable inversion and no plume. With ${\rho }_2<{\rho }_0<{\rho }_1$ the plume can carry some dust C > 0, up to the density limit ${\rho }_0\left(1+C\right)={\rho }_1$ . This gives the maximum concentration the inversion can hold with ${\rho }_0={\rho }_2$ ,

$C_{mp}=\Delta T/T_p$ (14)

T_p being the temperature of the plume in ˚Kelvin. At larger concentrations the plume is too heavy and sinks down in the lower layer. Laboratory experiments show this [20]. The underside of the plume will contain the heaviest dust load, the upper layers will have C < C_mp. This complicates matter considerably, see next Section 4.2.

An exception occurs if T_p is higher than the temperature of the overlying clean air or ρ₀ < ρ₂. Volcanic plumes may start out this way, but they will radiate heat to both sides and cool down quickly, and sink down to the inversion level. With ρ₀ = ρ₂ they have the same possibility to survive as the stable inversion without the dust. However, if the concentration in the plume is everywhere like Equation (14), the plume will be totally immersed in the lower layer which is unrealistic because the lower half will be the denser. An average concentration closer to C_mp/2 is more realistic for a neutrally buoyant plume, then it will float half immersed in the lower layer and have the same degree of stability as the inversion itself.

4.2. The Unstable Plume: Streak Fallouts

When T_p is considerably higher than the temperature of the overlying air (ρ₀ < ρ₂ < ρ₁), Equation (14) is still valid, but the plume hotter than all the surrounding air. This is common for volcanic plumes that can start out much hotter than the environment. In the case $C\underset{\_}{\approx }{C}_{mp}$ , the cooling of the plume will deplete the stability steadily and rapidly. The cooling will result in streak fallouts from the heavier lower half of the plume. [21] Figure 9, shows a measurement of a streak fallout. The streak fallouts are dust laden density currents of vertical air flow, subjected to entrainment of cold air (ρ₀ = ρ₂) that steadily increases the air density in them. They can run in a steady flow with a fixed densimetric Froude number (F_RΔ value) [15].

Streak fallouts are very common in the start of the neutrally buoyant flight of volcanic plumes, they can be seen in almost any picture of such a plume, see e.g. Figure 3. They are much more effective in removing mass load from the plume than ordinary fallout and they contain all grain sizes, not only the coarse grains.

The onset of dense gravity currents like the streaks, is an instability process that cannot be predicted. As the streak fallouts stop, the plume has become stable and Equation (14) is valid. Plumes will therefore start their neutrally buoyant flight with $C\le C_{mp}$ .

4.3. Gravitational Flattening

We can assume α = 1/2 in Equation (13) without any loss of generality. Then the neutrally buoyant plume is flowing with the upper half above the inversion (ρ = ρ₂) and the lower half below as before. The diffusivity in the plume may be K = 0 but the plume will spread to the sides never the less as the following analysis shows.

Figure 1 shows a density stratification in a stable plume drifting along with the wind velocity U coming out towards the reader of the paper. There is an overpressure in the center of the plume with respect to the ends as illustrated by the triangle in Figure 1(b). Using Navier-Stokes equation with Boussinesque’s approximation, (leave the density difference out everywhere except in the gravity term) the sideways flow, widening the plume, can be estimated. The half width, L in Figure 1, will increase with time from an initial value R₀ when downwind distance increases as x = Ut. Here, and in forthcoming text, x means the downwind distance, a meaning different from Figure 2. We have from Eliasson et al. 2014,

$L/R_0={\left[1.5\frac{s/U{T}_{pf}+0.733}{1+0.3{\left(L/R_0\right)}^{-4}}\right]}^{2/3}$ (15)

Figure 2. Two Lagrangian particle tracking dispersion models. Upper panel: Traditional diffusion model, no pancaking and high K. Lower panel: New model with pancaking and low K. Black lines: Downwind coordinate. Colored lines: The track of each particle.

$T_{pf}={\left(R_0/\Delta g\right)}^{1/2}$ (16)

T_pf is a time constant of the spreading, sometimes called pancaking, of the plume, but s is the downstream coordinate from the point where L = R₀, not to be mixed up with the location of the source at x = 0. In that situation the velocities in the deformation of the plume will be $\text{d}L/\text{d}t=V_L=V_H=R_0/T_{pf}=\sqrt{\Delta g{R}_0}$ , both sideways and up in the bottom and down in the top. This is the deformation velocity in the initial condition L = R₀, it can be compared to $V_L=\sqrt{\Delta gh}$ for gravitational flow in [22], and it shows that the two expressions compare.

Further explanations of gravitational deformation are in [21] and the assumptions made to derive Equation (15). The plume will spread sideways and be compressed in the vertical direction. The compression will not change the dust concentration. The spreading is the work of the pressure gradient in Figure 1 alone, as mixing and fallout is disregarded. This assumption is justified if streamwise derivatives (direction of U) within the plume are small compared to derivatives in the other directions. It will hold to a sufficient degree of approximation for fine dust plumes with low K and ordinary fallout, but not for coarse particle plumes and heavy streak fallout. In using Equations (15) and (16) as a subgrid model, T_pf and L/R₀ are recalculated in each time step.

Finally, the use of Bernoulli’s equation assumes no flow resistance. It is therefore suggested to introduce a correction factor, 0 < B < 1, into the equation for T_p and use BΔ instead of Δ in Equation (10).

5. Dispersion Models

5.1. The Models and Their Plume Physics, Mathematics and Computational Techniques

There are numerous computer models designed for simulating the dispersion of dust in the atmosphere, usually called atmospheric dispersion models¹. On top of that NWP (Numerical Weather Prediction) models are beginning to include diffusion. The plume physics these models are based on, are almost without exception Gaussian diffusion resulting in Equations (2) - (5). In steady flow Equation (10) is the governing differential equation of the dispersion, it is sometimes extended to all three spatial dimensions. This kind of modelling can involve very difficult mathematics and solution procedures.

There are two main systems for model simulating dispersion, Eulerian models and Lagrangian models. When the full flow equation (Navier-Stokes) is used, the dust concentration may be included in the gravitational term in the flow equations as is the case in Equation (12). Eulerian systems thus can include the gravitational term that varies with dust concentration. Then the buoyancy Equation (15) will in theory be correctly adjusted in the computational routine solving for the velocity of the spreading.

However, this kind of modelling involves great difficulties. Very small time steps are necessary to obtain reasonable accuracy. In [10] a time step of 0.25 seconds is used and it is found necessary to use an associated time step of 0.025 seconds. These time steps and associated sizes of the computational grid, are too small for a model used to simulate clouds maps for the North Atlantic region, but also too large to simulate streak fallouts. Measurements show that dust clouds are very patchy and dust concentrations fluctuate heavily, [6] [8]. It is for all practical reasons not possible to make the Eulerian computational nets fine enough to simulate plumes that are kept afloat for several days in a run of several hundreds of kilometers and can produce the same turbulent fluctuations of C as the measurements show. Nor is it possible to compute the horizontal pressure gradients from the center plume and outwards, responsible for the gravitational flattening, with acceptable accuracy due to the limitations of the spatial grid. The common approach to solve these problems is to keep all fluctuations in a subgrid model utilizing suitable turbulence closure rules. For the same reason the gravitational flattening has to be in a subgrid model as well, utilizing Equations (15) and (16).

Another difficulty is that NWP models are actually simple CFD (Computational Fluid Dynamics) models. It is a well-known fact that they are all unable to compute turbulent fluctuations correctly. This has to affect the estimation of the fallout velocities seriously, and thereby the fallout estimations as Equation (9) shows. In the long run, the computational errors in the concentrations of the model plume will increase in each time step and effect the subsequent changes in the mass flux of the plume.

Many models use the wind vectors available in the international meteorological databases, then there is no transverse velocity from gravitational flow. This should be included; it may be seen from Equation (12). In steady flow, the differences in C in horizontal direction will cause density gradients that induce a flow the models cannot predict. Fluid layers of different density are in equilibrium in horizontal layers only with the less dense fluid on top of the denser. This principle is common knowledge of course, but an investigation is needed to clarify what errors it causes, not to include gravitational flow caused by differences in dust concentrations in the modelling.

Lagrangian models, or puff models, use the NWP wind vectors and particle tracking mathematics. The mixing effect of the turbulence can be modeled with a displacement vector alluding random walk around the NWP streamline. It is added to the displacement in each time step. After a number of time steps the resulting summed-up displacement vector of each particle, will take on the Gaussian distribution. This is a consequence of the central limit theorem, no matter what is the original statistical distribution behind the individual displacements in each time step.

Some pollution dispersion models make use of this property and simply predict the pollution to be Gaussian distributed around the NWP streamline through the source. Then Equation (2) and Equation (4) are used to map the resulting pollution concentrations.

5.2. Comparison of Models

Figure 2 explains the difference between the Gaussian process and the gravitational dispersion (pancaking). In this figure, a Lagrangian particle tracking model with a large dispersion coefficient and no gravitational flattening is used above the centerline, but a small dispersion with normal gravitational flattening is assumed in the figure below. The average rate of the widening of the dust cloud is almost the same in both pictures.

In Figure 2 the effective width, derived from Equation (3) is made to match the L of Equation (15) as well as possible. Figure 2 demonstrates nicely that this can be done for a while, how long it works depends on the length scales UT_pf and $\sqrt{\text{2}\pi }\sigma$ . The match is difficult to see in Figure 2 as K is about 20 times bigger in the upper panel flow. In the lower panel the weak undulations of the boundaries indicate a mixing layer between the dust plume and the clean air. In the long run the matching will break. Taking the limit for large downstream distances s, we will find this relation.

$\sigma /L~K^{1/2}{U}^{1/6}{R}_0^{1/6}{\Delta }^{-1/3}{g}^{-1/3}{s}^{-1/6}$ (17)

This shows that the σ/L ratio has no fixed limit, even in steady wind (U constant) it goes down with Δ and s. Gravitational plumes will therefore spread more quickly than diffusion plumes, and in the same time get thin. As fallout will be largest from the middle part, this can result in a rather strange hole in the middle of the plume. Further pancaking can therefore lead to formation of patches. Undulations in the wind direction and wavy flow will have similar effect.

Another effect that may be derived from Figure 2 is the particle tracks that will be created if there is no turbulence. In that case the undulations of the lines will disappear and the tracks will become straight. In the upper panel the particle tracks will become straight lines parallel to the wind direction, there is no widening of the plume. In the lower panel the widening will be mostly unaffected when the turbulence disappears and the particle tracks become almost straight lines, the only curvature being the result of the application of Equations (15) and (16).

Another possibility is to make the upper panel undulation the same order of magnitude as the lower panel. This corresponds to making the horizontal diffusion coefficient the same in both panels of the picture. Then the dispersion of the plume is clearly underestimated in the upper panel, as the result of not using the gravitational dispersion Equations (15) and (16) represent.

6. An Example: Eyjafjallajokull 2010

6.1. The Development of the Plume

The eruption in Eyjafjalljokull 2010 was not very big, but it had great consequences and an uncountable number of papers have been published about it, many satellite photos of it exist² and many simulations of the dust cloud have been made. It is therefore a good example to take. Figure 3 shows a picture of the plume 12.5. 2010. The streak fallouts are clearly visible. The wind conditions in Figure 3 lasted for some time. Figure 4 is a MODIS picture of the eruption May 10^th 2010. The rather sharp boundary between the plume and the clean air wavy, indicating penetrative convection [15] of the dusty air into the clean air as the plume’s boundary is being pushed outwards from the middle.

Figure 3. Streak fallouts.

Figure 4. Eyjafjallajokull May 10^th 2010 (MODIS).

There are also waves on the top of the plume. This could very well be the Kelvin-Helmholzt vortices one can see on close up pictures [7] on top of the plume. Later in the flight, where the plume begins to sail under the white rain clouds, the tendency to develop a thin section or a hole can be seen, and in the same time a thinning out of the boundary layer at the edge of the plume as in Figure 2. All together the plume can be traced about 200 - 300 km out in the sea, but then it disappears. The longest run of the Eyjafjallajökull plume was about 1000 km, but then the plume was very thin and patchy. Satellite photos from April 15 on display on² all show this. The models used by the VAAC system to simulate volcanic plumes have improved a little but the results still show too big clouds. This may be seen in [9]. But their result is probably the best that can be done without using gravitational flattening, maximum carrying capacity and streak fallouts

This example indicates gravitational dispersion rather than diffusion until the boundary layer around the plume has grown enough to reach vertically through the plume. As this takes several hundred kilometers, the diffusion coefficient K is very low, order of magnitude 10 - 20 m²/s. The gravitational deformation can explain the sideways spreading in Figure 4. This is shown in [21] citing [23] a paper in Icelandic. That investigation reveals that a ΔT of few degrees is sufficient to explain the gravitational deformation. A diffusion with the same effect would need a K value of around 500 m²/s.

A closer investigation of Figure 4 shows the plume to be remarkably stable, the white cover of normal clouds does not indicate similar stability, probably because they are not trapped in an inversion. Never the less, it is quite clear that diffusion will overtake the dispersion process in the long run. But little may be left of the plume when that happens.

6.2. Smoke Stack Plumes

Models of smoke stack plume models utilize the system Equations (2) - (4) with a large cross-stream K in the horizontal plane, resulting in a larger σ_y than and σ_z pancaking the plume. They also use parameters to define the stability of the air, it is sufficient to refer to the Briggs plume model as an example of how this kind of modelling works [24] [25].

Stack plumes will be in the shear layer where the logarithmic velocity profile dominates, but that does not explain a large cross-stream K. But if it is there, by some reason or another, e.g. the wake of the stack can cause that effect, the diffusion will result in a cross-stream gradient of C, dC/dy from the center of the plume outwards. If the pressure in the vertical is hydrostatic, as all models do assume for a neutrally bouyant plume, then Equation (12) gives rise to an outwards pressure gradient if there is dust in the plume. If there are just gasses in the plume we will get the same result for a plume that starts out hot and thins down in the downstream direction.

Smoke stack plumes will thus be subject to gravitational flattening, acting on top of the diffusion. Tropospheric plumes are larger, so gravitational forces are playing a bigger part because of the size, so diffusion can be quite absent and the plume will spread to the sides anyway.

7. Conclusions

The foregoing physical evidence and argumentation leads to the following conclusions.

Gravitational dispersion, sometimes called pancaking, is a sideways dispersion of dust plumes, controlled by gravitational forces acting as a pressure gradient in the cross-stream direction in the vertical plane away from the center line of the plume. Gravitational dispersion can flatten plumes without mixing. In the case of no diffusion, such a plume that is initially circular, and subjected to advection by a steady wind, will become deformed to an elliptic form with the long axis in the horizontal plane, perpendicular to the wind vector. In the case of diffusion, gravitational dispersion is also active.

There will be a mixing layer between the clean air and the dusty air. It will spread to the main plume eventually, but that can take a long time while the plume migrates hundreds of kilometers.

Diffusion only is possible, but not in dust plumes, only if the pollution is a passive substance with the same density as the ambient air. When different mixing ratios of the pollution cause different densities, the density gradients will cause a flow in the horizontal plane.

Streak fallouts are vertical currents of dust laden air, visible in many pictures of migrating plumes. In young dust plumes, carried upwards in buoyant flow of warm air, streak fallouts will be active from the start. When neutral buoyancy is reached, it is possible to estimate the dust carrying capacity of the plume, provided the temperatures of the plume and the ambient atmosphere are known from measurements. The carrying capacity is the maximum concentration of dust in the plume that can be sustained by the temperature differences of the plume and the clean air. If dust concentrations higher than the capacity are carried up to the neutral buoyancy level, streak fallouts from the underside of the plume will scavenge them out, taking a lot of the fine grains with them and reducing the mass flux of the plume.

Ordinary fallout scavenges coarse grains out of the plume in a steady stream controlled by the terminal velocity of the grain fractions. The terminal velocity of the coarse grains is influenced by the turbulence, so estimation of ordinary fallout, based on terminal velocity in still air, may be very inaccurate.

Models that account only for the mixing process of the turbulent diffusion will inevitably contain too much dust, become too big and cover too large an area. Inclusion of the carrying capacity, streak fallouts and gravitational deformation is a much-needed improvement for such models. Without it, models of tropospheric plumes use the same plume physics as smoke stack models and the plumes become too big.

Acknowledgements

Thanks are due to my colleagues in the Disaster Prevention Institute, Kyoto University, Japan; University of Applied Sciences in Düsseldorf, Germany; Department of Engineering, Reykjavik University, Iceland. Furthermore, thanks are due to Dr. Ólafur Rögnvaldsson in Belgingur Ltd., Dr. Gylfi Árnason and others for their invaluable help in measuring and modeling plumes with me.

NOTES

¹https://en.wikipedia.org/wiki/Atmospheric_dispersion_modeling.

²https://www.flickr.com/photos/gsfc/albums/72157623862023918/.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Taylor, G.I. (1921) Diffusion by Continuous Movements. Proceedings of the London Mathematical Society, Series A, 20, 196-211.
[2]	Tiesi, A., Villani, M.G., D’Isidoro, M., Prata, A.J., Maurizi, A. and Tampieri, F. (2006) Estimation of Dispersion Coefficient in the Troposphere from Satellite Images of Volcanic Plumes: Application to Mt. Etna, Italy. Atmospheric Environment, 40, 628-638. https://doi.org/10.1016/j.atmosenv.2005.09.079
[3]	Toon, O.B., Pollack, J.B. and Sagan, C. (1977) Physical Properties of the Particles Composing the Martian Dust Storm of 1971-1972. Icarus, 30, 663-696. https://doi.org/10.1016/0019-1035(77)90088-4
[4]	Kulkarni, P., Baron, P.A. and Willeke, K. (2011) Aerosol Measurement: Principles, Techniques, and Applications. John Wiley & Sons, Hoboken. https://doi.org/10.1002/9781118001684
[5]	Mackie, S., Cashman, K., Ricketts, H., Rust, A. and Watson, M. (2016) Volcanic Ash: Hazard Observation. Elsevier, Amsterdam.
[6]	Weber, K., Eliasson, J., Vogel, A., Fischer, C., Pohl, T., Van Haren, G., Meier, M., Grobéty, B. and Dahmann, D. (2012) Airborne In-Situ Investigations of the EyjafjallajÖkull Volcanic Ash Plume on Iceland and over North-Western Germany with Light Aircrafts and Optical Particle Counters. Atmospheric Environment, 48, 9-21. https://doi.org/10.1016/j.atmosenv.2011.10.030
[7]	Schumann, U., Weinzierl, B., Reitebuch, O., Schlager, H., Minikin, A., Forster, C., Baumann, R., Sailer, T., Graf, K., Mannstein, H. and Voigt, C. (2011) Airborne Observations of the Eyjafjalla Volcano Ash Cloud over Europe during Air Space Closure in April and May 2010. Atmospheric Chemistry and Physics, 11, 2245-2279.
[8]	Eliasson, J., Yoshitani, J., Miki, D., Weber, K., Boelke, C. and Scharifi, E. (2016) Measurements of Particle Distribution and Ash Fluxes in the Plume of Sakurajima Volcano with Optical Particle Counter. Journal of Disaster Research, 11, 85-95. https://doi.org/10.20965/jdr.2016.p0085
[9]	Osman, S., Beckett, F., Rust, A. and Snee, E. (2020) Sensitivity of Volcanic Ash Dispersion Modelling to Input Grain Size Distribution Based on Hydromagmatic and Magmatic Deposits. Atmosphere, 11, 567. https://doi.org/10.3390/atmos11060567
[10]	Burton, R.R., Woodhouse, M.J., Gadian, A.M. and Mobbs, S.D. (2020) The Use of a Numerical Weather Prediction Model to Simulate Near-Field Volcanic Plumes. Atmosphere, 11, 594. https://doi.org/10.3390/atmos11060594
[11]	Zeng, Z., Wang, Y., Duan, Y., Chen, L. and Gao, Z. (2010) On Sea Surface Roughness Parameterization and Its Effect on Tropical Cyclone Structure and Intensity. Advances in Atmospheric Sciences, 27, 337-355. https://doi.org/10.1007/s00376-009-8209-1
[12]	Emeis, S. (2014) Current Issues in Wind Energy Meteorology. Meteorological Applications, 21, 803-819. https://doi.org/10.1002/met.1472
[13]	Crank, J. (1979) The Mathematics of Diffusion. Oxford University Press, Oxford, p. 572.
[14]	Garratt, J.R. (1994) Review: The Atmospheric Boundary Layer. Earth-Science Reviews, 37, 89-134. https://doi.org/10.1016/0012-8252(94)90026-4
[15]	Pedersen, F.B. (1986) Environmental Hydraulics: Stratified Flows. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-86600-5
[16]	Mahrt, L. (2014) Stably Stratified Atmospheric Boundary Layers. Annual Review of Fluid Mechanics, 46, 23-45. https://doi.org/10.1146/annurev-fluid-010313-141354
[17]	Clayson, C.A. and Kantha, L. (2008) On Turbulence and Mixing in the Free Atmosphere Inferred from High-Resolution Soundings. Journal of Atmospheric and Oceanic Technology, 25, 833-852. https://doi.org/10.1175/2007JTECHA992.1
[18]	Wilson, R. (2004) Turbulent Diffusivity in the Free Atmosphere Inferred from MST Radar Measurements: A Review. Annals of Geophysics, 22, 3869-3887. https://doi.org/10.5194/angeo-22-3869-2004
[19]	Stevenson, J.A., Millington, S.C., Beckett, F.M., Swindles, G.T. and Thordarson, T. (2015) Big Grains Go Far: Understanding the Discrepancy between Tephrochronology and Satellite Infrared Measurements of Volcanic Ash. Atmospheric Measurement Techniques, 8, 2069-2091. https://doi.org/10.5194/amt-8-2069-2015
[20]	Manzella, I., Bonadonna, C., Phillips, J.C. and Monnard, H. (2015) The Role of Gravitational Instabilities in Deposition of Volcanic Ash. Geology, 43, 211-214. https://doi.org/10.1130/G36252.1
[21]	Eliasson, J., Yoshitani, J., Weber, K., Yasuda, N., Iguchi, M. and Vogel, A. (2014) Airborne Measurement in the Ash Plume from Mount Sakurajima: Analysis of Gravitational Effects on Dispersion and Fallout. International Journal of Atmospheric Sciences, 2014, Article ID: 372135. https://doi.org/10.1155/2014/372135
[22]	Bonadonna, C. and Phillips, J.C. (2003) Sedimentation from Strong Volcanic Plumes. Journal of Geophysical Research, 108, 2340. https://doi.org/10.1029/2002JB002034
[23]	Andradóttir, H.G. (2010) Lárétt útbreiesla gosstróka EyjafjallajÖkuls reiknue frá gervihnattamyndum (Horizontal Extent of Volcanic Plumes from EyjafjallajÖkull Inferred from Satellite Photographs), 22 árbók VFí/TFí 2010 (in Icelandic with English Abstract).
[24]	Briggs, G.A. (1972) Discussion: Chimney Plumes in Neutral and Stable Surroundings. Atmospheric Environment, 6, 507-510. https://doi.org/10.1016/0004-6981(72)90120-5
[25]	Beychok, M.R. (2005) Fundamentals of Stack Gas Dispersion. 4th Edition.