New Model for Dispersion of Volcanic Ash and Dust in the Troposphere

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 occurs. The mo-tivation of this paper is to improve the existing models for dispersion of dust plumes. The traditional models using Gaussian diffusion do consider vertical diffusion negligible, but include coefficients for horizontal diffusion large enough to explain the sideways dispersion, so clearly seen on satellite images. It is found that gravitational flattening caused by atmospheric density stratification may be the main cause of dispersion in dust plumes above the turbulent boundary layer. Additionally, this causes a dust plume in between two layers of small temperature difference to have a certain carrying capacity for dust load. The corresponding mass loading can be estimated from the temperature difference between the layers above and beneath the plume. Dust plumes having a mass load in excess of this carrying capacity will be forced to jettison the extra load. This may be seen as streak fallout from the plume. In the same time, the plume will be subjected to flattening to the sides, caused by density-controlled pressure gradients, in addition to any diffusion if there is any. The length scale of the plume width resulting from this 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 de-velopments 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.


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

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 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 10 /U * = 15 is a good example [11] and we have ( ) ( ) 10 1 0.167 ln 10 U U y = + .
In the elevation 80 -100 meters this means about 35 % higher wind velocity than U 10 . 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 2 /s in the elevation range y = 10 -100 m for U 10 = 10 m/s. If this is extended up to 10 kilometers the velocity will be 115% higher and ε will become 2500 m 2 /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: where y and z are the horizontal and vertical coordinates perpendicular to the direction of the advection (x-axis ). C 0 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.
( ) It shows that 2 σ π denotes the effective width and height of the plume with respect to C 0 . Associated with σ is the coefficient of dispersion K. It is really a vector, but for demonstration purposes it is assumed a constant For times longer than T T , the integral time scale of the turbulence, we have And similar for the other directions of the flow. The  (4) and (5). When the mean velocities in y and z are V = W = 0, the non-zero velocity fluctuations , u v ′ ′ and w′ are all related through the continuity equation, is a little unrealistic. High up one would expect isotropic turbulence with 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′ 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.

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 ( ) 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  [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 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.

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. 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 2 /s up to 20 km, and only occasionally above 10. In [18]  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.

Settling Velocity
Most dust particles consist of silica-rich minerals with density of around 2500 kg/m 3 . 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 2 S V . 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 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 ( ) 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.

Fallout
But whatever the foregoing result is, we still have a V S for each grain size fraction International Journal of Geosciences and the concentration in each vertical is determined by the differential equation for the particle flow q(z) inside the plume, 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)  If K is very big dC/dz will be very small. If we have as initial condition an average concentration C 0 at t = 0 in the vertical, we will have a mass loading of dust M C = C 0 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 B V 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 Equation (11) supposes a concentration profile ( ) ( ) 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 0 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)  Due to the large difference of the densities of the dust and the clean air, the plume's density is ( ) in any point, ρ 0 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 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 ρ 0 < ρ 2 . 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 ρ 0 = ρ 2 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

The Unstable Plume: Streak Fallouts
When T p is considerably higher than the temperature of the overlying air (ρ 0 < ρ 2 < ρ 1 ), 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 mp C C ≈ , 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 (ρ 0 = ρ 2 ) 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 mp C C ≤ .

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 (ρ = ρ 2 ) 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  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)

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 1 . On top International Journal of Geosciences 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.  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 2 σ π . 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. ∆ (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.

Comparison of Models
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   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 2 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 2 /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 2 /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.

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 International Journal of Geosciences 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.

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.