The Development of Volcanic Ash Cloud Layers over Hours to Days Due to Atmospheric Turbulence Layering

: Volcanic ash clouds often become multilayered and thin with distance from the vent. We explore one mechanism for the development of this layered structure. We review data on the characteristics of turbulence layering in the free atmosphere, as well as examples of observations of layered clouds both near-vent and distally. We then explore dispersion models that explicitly use the observed layered structure of atmospheric turbulence. The results suggest that the alternation of turbulent and quiescent atmospheric layers provides one mechanism for the development of multilayered ash clouds by modulating vertical particle motion. The largest particles, generally > 100 µ m, are little affected by turbulence. For particles in which both settling and turbulent diffusion are important to vertical motion, mostly in the range of 10–100 µ m, the greater turbulence intensity and more rapid turbulent diffusion in some layers causes these particles to spend greater time in the more turbulent layers, leading to a layering of concentration. The results may have important implications for ash cloud forecasting and aviation safety.


Introduction
Volcanic ash is a multi-billion dollar economic hazard to aviation, as shown during the 2010 eruptions of Eyjafjallajökull, Iceland [1,2]. It is also a risk to flight safety, with hundreds of encounters of varying severity recorded, and several instances of multiple engine flame-out in flight. The International Airways Volcano Watch (IAVW), which seeks to safely separate aircraft from volcanic ash in flight, relies on detecting areas of ash and forecasting its future movement [3]. The IAVW acts through nine Volcanic Ash Advisory Centers (VAACs). In case of eruption that potentially affects flight safety, the VAACs issue both a written Volcanic Ash Advisory (VAA) and a Volcanic Ash Graphic (VAG), which state and show the analysis of the current position of an ash cloud and the forecast position to 18 h from present. Ash cloud positions are rendered as polygons for separate aircraft flight levels or altitude ranges, and flight paths are adjusted accordingly. However, the model forecasting of ash presence and concentration is generally poorly resolved vertically, although there is some progress in this direction, for example, Ref. [4,5]. Aircraft flying in a supposedly ash-contaminated region at a particular altitude may encounter no ash or significant and potentially damaging amounts, due to the high degree of ash stratification with altitude. The improved understanding and forecasting of stratification would assist enormously in managing the hazard and support the continuing development of the IAVW.
Photography and satellite imagery of numerous volcanic eruptions show that a single or multiple layered structure is a fundamental aspect of near-source volcanic cloud development [6,7], with a single, volcanic umbrella cloud layer being dominant in most of Figure 1. Aerial as well as ground-based photography have been particularly useful near the volcano to elucidate this layering. Lidar backscatter data have been key in defining the more complicated layered structure in distal regions, particularly those from the CALIOP (Cloud-Aerosol LIdar with Orthogonal Polarization) instrument ( Figure 2). Volcanic layers can be stratospheric as well as tropospheric [8]. Aqua-MODIS RGB false color image [9] of North Atlantic capturing this ash cloud (pink hues). In distal regions, volcanic clouds tend to be thin, less concentrated, and often broken into multiple bodies.
Multiple layers are thought to form and separate by numerous processes, some of these unique to volcanic clouds. The primary volcanic cloud layer near the vent, such as the volcanic umbrella cloud or anvil cloud, arises from the driving of hot eruptive gas and ash parcels outward around their equilibrium level, or neutral buoyancy height [10]. As suggested by brightness temperatures over the surface of near-vent clouds and ground based photography, umbrella clouds can be solitary, accompanied by a single lower intrusion resulting from re-entrainment and column-edge downflow [11,12] or accompanied by lower level skirt clouds [13], which may or may not contain ash. Ash accretion, ash re-entrainment and source variability-injection of ash at different altitudes with changing eruption rate and wind fields, cause the development of multiple layers [7,14,15]. Gas-ash separation of particulate-laden and gaseous volcanic cloud components  [9] of North Atlantic capturing this ash cloud (pink hues). In distal regions, volcanic clouds tend to be thin, less concentrated, and often broken into multiple bodies.
Multiple layers are thought to form and separate by numerous processes, some of these unique to volcanic clouds. The primary volcanic cloud layer near the vent, such as the volcanic umbrella cloud or anvil cloud, arises from the driving of hot eruptive gas and ash parcels outward around their equilibrium level, or neutral buoyancy height [10]. As suggested by brightness temperatures over the surface of near-vent clouds and ground based photography, umbrella clouds can be solitary, accompanied by a single lower intrusion resulting from re-entrainment and column-edge downflow [11,12] or accompanied by lower level skirt clouds [13], which may or may not contain ash. Ash accretion, ash re-entrainment and source variability-injection of ash at different altitudes with changing eruption rate and wind fields, cause the development of multiple layers [7,14,15]. Gas-ash separation of particulate-laden and gaseous volcanic cloud components has been noted as a further cause of volcanic layer formation [16,17], perhaps enhanced by gravity current slumping of the particle laden component [18]. Double diffusion and convective sediment flux to a single [11,12,[19][20][21] and multiple [6] levels by descending fingers that intrude below the level of a major volcanic cloud layer have been observed, and recreated in the laboratory.
Further from the vent, ∼>500 km and hours in transport time, downward looking satellite sensors are used for cloud tracking, but these penetrate only partially into the highest layer of an optically thick cloud, and satellite remote sensing algorithms are most sensitive to the column integrated ash properties, even when multiple optically thin layers are present. Although horizontal, planview resolution can thus be good, satellite data and often volcanic ash transport and dispersal models (VATDs), used to forecast future position, have difficulty in reproducing the vertical structure of distal volcanic clouds [4,22,23]. Numerical inversion of VATDs has been used to map distal cloud positions seen in CALIOP and other data to time-varying release height from the source vent [5,[24][25][26]. Dacre et al. [27] used a new scheme to calculate clear-air turbulence (CAT) in the NAME VATD model. Numerous distal clouds from the Eyjafjallajökull eruption were analyzed and thought to be controlled in thickness by an interaction between wind shear, which acted as thin layers, and turbulent diffusion, acting to thicken them [27].

Problem Statement
A correct understanding of layer formation and morphology is critical for any attempt to construct VATD models capable of producing output consistent with observations of vertical ash distribution, yet the formation mechanisms of atmospheric ash cloud layers are not fully understood. One potential mechanism for their formation and characteristics is the subject of the present contribution. Our working hypothesis is that, because particle settling and atmospheric turbulence act to varying degrees in a layered atmosphere, the turbulence structure of the atmosphere can cause ash layer formation, through the process of enhanced suspension in vertically restricted regions of high turbulence.

Atmosphere
Ash transport occurs in the large-scale structures of the wind field, while ash dispersion or spread is caused by the small-scale turbulence structures and eddies [12,28]. The free atmosphere itself, with or without volcanic clouds, generates small-scale turbulent eddies. Eddies are generally created in the troposphere by Rayleigh-Taylor (convective) instabilities associated with clouds, and can be found 33% of the time [29]. In the stratosphere, clear-air turbulence (CAT) is primarily thought to be generated by breaking of upward-propagated gravity waves from the troposphere or tropopause, i.e., the Kelvin-Helmholtz mechanism, acting during episodes of vertical wind shear in the horizontal wind components [30][31][32]. In volcanic clouds near the vent, turbulence is created by both the Rayleigh-Taylor and Kelvin-Helmholtz mechanisms, as ash clouds intrude into the atmosphere as gravity currents. Kelvin-Helmholtz instability is driven by the shear between the intruding cloud and the atmosphere [33]. Rayleigh-Taylor instability is driven by cloud top convective instabilities [29], convective sedimentation, fingering and local, eddy scale density reversal [11,16,34].
In the free atmosphere, parameters such as moisture content and temperature do not change monotonically with height; there are layers of relatively homogeneous air, separated by regions in which parameters vary rapidly [29,35]. Both the stratosphere and the troposphere are layered in this way on scales of O[0.1-1 km] [36,37], although the layers in the troposphere tend to be more transient and discontinuous [38]. Information of sufficient resolution in the vertical direction to discern and characterize the turbulence layered structure is obtained from airborne measurement campaigns or rawinsonde balloon releases [39][40][41]. Several methods have been developed to derive turbulence from rawinsonde and other high-resolution data. Vasseur and Vanhoenacker [29] measured changes in the refractive index structure parameter for radio waves, as turbulence causes changes in the refractive index, based on rawinsonde pressure, temperature, humidity, wind speed and wind direction data. Clayson and Kantha [42] used variations in the potential temperature profile from an idealized profile to calculate the Thorpe scale, and derive turbulent dissipation and diffusivity. These estimates can be made for single rawinsonde profiles with simple calculations. Marlton et al. [43] have recently developed an accelerometer that can be added to the sonde instrument package and provide measures of instrument motion transformable into turbulence measures. In the troposphere, layers of constant relative humidity (RH) or mixing ratio, q, can act as tracers for high turbulence [40]. This is because as turbulence intensifies, mixing intensifies, and tracers become more evenly distributed throughout a layer. The macroscopic mixing caused by turbulence acts as a diffusive process, which is characterized by a (turbulent) eddy diffusivity, κ. (Atmospheric moisture is also important in aiding plume lift, especially in plumes from weak sources, or at low latitude, where the moisture content is high [10,44]. High RH layers might therefore indicate not only turbulence, but also particularly intense reservoirs of energy available for plume lift.) As a result of the layering, large-volume or bulk turbulence, and the eddy diffusivity, κ, that expresses its intensity, is highly anisotropic (κ h κ z ), and only within thin, well-defined layers is it approximately isotropic (κ h,local ≈ κ z,local ) [38].
There are drawbacks to extrapolating high-resolution or point data to a regional scale because of spatial and temporal inhomogeneity; the troposphere is highly transient and spatially variable [42,45]; tropospheric isobaric surfaces are not necessarily parallel to the earth's surface, especially at fronts and in mountain waves [35]. Fronts are associated with tropopause folds, non-horizontal isobaric surfaces separating cold from warm air, and turbulence in folds is generated by local dynamic and convective instabilities. Mountain waves form as the density stratified atmosphere flows past the lee side of a mountain or mountain range. These waves can break, resulting in local turbulence concentrated in non-horizontal layers.
We will use the high-resolution atmospheric data of Cho et al. [40] as the basis for testing a model of ash cloud layer development from the turbulence structure of the atmosphere. To explore a case study of the Pinatubo 1991 volcanic clouds, we will speculate on potential atmospheric layers based on the considerations outlined in this section.

Ash Clouds
Data are available from a number of sources on the shape and structure of both nearvent and distal ash clouds. In the near-vent region, data tend to be more limited, due to the lower optical depth. Nevertheless, cloud brightness temperature (BT) as measured from nadir-looking geostationary and low-earth orbiting satellites provides useful information, as the topography of the top of the near-vent clouds can be quite variable, with a distinct high point or swell above the central vent that might be many kilometers above the top of the main umbrella cloud [17]. In addition, airborne and ground-based photography and videography have provided extensive data on the features at the base of the main umbrella or anvil, and within the underlying cloud layers. Visible satellite imaging of the near vent cloud top consistently reveals strong, well-defined three-dimensional vortex structures above the vent, which evolve to smooth, somewhat more diffuse structures in the umbrella cloud [46].
Although the air can be choked with opaque, diffuse ash bodies that extend to ground level near vent, and although gravitational intrusions, such as umbrella clouds, are wedgeshaped by nature and therefore of variable depth, measurements have been made of the ∆BT between the highest and lowest points on the cloud top, as indicated by the lowest and highest BT, respectively. The results suggest that umbrella clouds at the vent typically encompass depths of O [5] km (Table 1), which makes a large mass of ash available for transport at multiple levels. Some of the lower near-vent clouds also become distal clouds. In some cases, it might be possible to find the entire troposphere and even lower stratosphere charged with ash, or only a distinct layer or two of O[1-10] km depth. In the distal region, airborne lidar, EARLINET-AERONET (European Aerosol Research LIdar NETtwork-AErosol RObotics NETwork) and CALIOP data typically show thinner, more discontinuous cloud structures ( Figure 2) than is commonly seen near-source ( Figure 1). Three separate tabulations of distal ash cloud layer data for the Eyjafjallajökull plume have been published [49][50][51], and the EARLINET data have been summarized and modeled [27]. The data suggest that distal clouds from this tropospheric eruption were typically 0.3-3 km thick, made up of 2-3 layers, with individual layers of 0.3-1.4 km depth, and maximum age of 129 h (<1 week) ( Table 2). Dacre et al. [27] estimated a mean and standard deviation of the measured depths of 1.2 ± 0.9 km. The number density of particles in the clouds as they propagated over Europe generally peaked in the micron to submicron range, while the mass density peaked near 10 µm [51], although the larger size modes may be underrepresented [52].
Vernier et al. [53] using CALIOP data, discerned two or more well-defined layers in the cloud from Puyehue-Cordon Caulle three weeks after the eruption. Some of the layers showed fold or wrap-around structures ( Figure 1a in Vernier et al. [53]), perhaps related to vertical-plane chaotic mixing [54]. Clouds were up to 3 km thick, with individual layers of 0.1-2 km depth, in the upper troposphere-lower stratosphere (UTLS), centered on the tropopause at 8-14 km altitude.
On 12-13 July 1991, 26 days after the last major eruption of Pinatubo, a lidar flight noted numerous stratospheric layers [8]. The data showed a number of well-defined layers of 0.5-1 km depth between about 14 and 25 km altitude (Figure 1 in Winker and Osborn [8]). Along much of the line of flight there were two layers (22 and 25 km), but in places there were up to five. Within the first day, modeling by Fero et al. [17] is consistent with a mean grain size of 90 µm in the 22-km layer. Even after nearly a month, high backscatter lidar depolarization ratios at the base of the 22-km layer were still inexplicably consistent with particles in the 10-100 µm range [8]. Although depolarization ratios in the 25-km layer were consistent with mostly sulfur aerosols, bothWinker and Osborn [8] and Fero et al. [17] mention that fine silicate particles occurred with the aerosols.
In all studies cited above, distal ash layers were horizontal or tilted relative to the horizon, and had extents in the cross-transport direction of hundreds of km in the troposphere (Eyjafjallajökull), to thousands of km in the stratosphere (Puyehue, Pinatubo). Other consistent features are thinning, more sharply defining, and breaking into multiple, more-distal clouds of near-vent volcanic clouds.
Back trajectories of distal ash clouds for Eyjafjallajökull and Puyehue-Cordon Caullé are generally consistent with theoretically possible cloud heights at the source [49,53]. However, it is not clear that ash was injected at these altitudes by the eruptions, given uncertainties in vertical parcel motion and settling speed, or lack of incorporation thereof in the models [53,55]. This observation suggests that settling might not be the sole, or even most important, factor in accounting for vertical positions of ash clouds in the distal region. Dacre et al. [27] suggested that an increase and decrease in ash cloud thickness are likely due to a balance between vertical wind shear and turbulent diffusion for distal clouds comprised of mostly micron and submicron particles.

Methods
The atmospheric observations summarized in Section 2.1.1 show that both the troposphere and stratosphere are layered in turbulence intensity on a scale of fractions of a kilometer to several kilometers, often due to both cloud top turbulence and vertical wind shear generation. The ash cloud observations summarized in Section 2.1.2 show that more distal, distinct ash clouds evolve from more proximal, diffuse and thicker clouds.
In many VATDs, sub-grid dispersal in the vertical direction is usually described by a single vertical diffusivity, κ z , which may take on the same value as the horizontal diffusivity, κ h ; the result is then uniform or isotropic ash dispersion. Our goal here is to contrast the physical behavior of an ash cloud under conditions of isotropic (κ h = κ z = const.) or constant vertical turbulence (κ z = κ z (z)), and layered turbulence. We base our methodology on analytical and numerical models using synthetic atmospheres with and without multiple turbulent layers separated by relatively quiescent air.

Eulerian Analytical Formulation
In the case of isotropic turbulence, we begin by assuming Cartesian coordinates, (x, y, z), with velocity components, (u, v, w). The three components of the turbulent diffusivity, (κ x , κ y , κ z ) are the same, κ. The concentration of particles in the i-size fraction, C i , can be written for time, t and space as: where Φ represents the source/sink function, which in the case of ash clouds is mostly represented by aggregation and disaggregation of small particles. In the present case, such processes are set to zero. To derive an instructive solution, we assume a two-dimensional system in x and z (z upward, and motion in the y direction can be ignored) with a pointsource in time and space, w = w s , the settling speed is negative, and that, following a streamtube, the motion of the volcanic cloud can be characterized by a single downwind coordinate direction s-for which the axis is everywhere tangent to the plume centerline, for example, Ref. [56,57]-and speed U in that direction. Under these assumptions, the advection-diffusion equation becomes: with the well-known solution for the impulse initial condition [58,59]: It is reasonably clear that the solution is a Gaussian in (s, z), in which ash in the mean is blown downwind at speed U and falls at speed w s , while spreading from the mean at rate 4κt. The second case is for layered turbulence, in which κ z = κ z (z), and generally κ z = κ h . Since a 1D model is instructive in this system, we investigate vertical motions alone. In this system, there is a discontinuity in the diffusive flux at the lower boundary of a more turbulent layer because κ z decreases suddenly, and in this case, the one-dimensional advection-diffusion equation: becomes, in the region of the lower boundary of the upper layer: for κ z,upper κ z,lower . A step-like concentration gradient develops at the base of the upper layer, as some particles settle from the layer, while others are swept back into it by turbulence rather than settling through the boundary. The solution is well-known in volcanology and sedimentation research generally [21,60,61]: where t 0 is some initial time at which sedimentation starts in a given control volume, and w s is negative or downward.
In a more quiescent layer below the boundary, particles only settle and are advected downwind, there is little turbulence to enhance persistence within the layer. Thus, for systems in which particle motion is controlled by both turbulent diffusion and settling, turbulent layers can retain particles longer than do quiescent layers because of continuing re-entrainment in eddies (Equation (6)). Particles fall relatively rapidly through the quiescent layers because of unhindered settling (Equation (3)), sometimes even enhanced by the effects of convective sedimentation [19], which is not included in the present model.

Lagrangian Formulation
We can gain additional insight by looking at the Lagrangian formulation of the problem. In a Lagrangian framework, common to many ash dispersion models, such as HYSPLIT, movement is given by: where r is the position vector, v is the velocity vector and j is a time index, for which ∆t = t j+1 − t j . For the vertical, z component of the position: If turbulence is added, and assuming no mean vertical flow, this becomes: where w is the random vertical component of an instantaneous turbulence velocity. Note that the direction (up or down) and magnitude of w change with time, and therefore the effective settling speed, w s + w , can be positive or negative (rising parcel) and of widely varying magnitude, depending on the ratio of w s to w . Furthermore, w can be related to the turbulent diffusion by: where L o is the length scale of the largest eddies, the Ozmidov scale, which must be less than the layer thickness, i.e., L o < h, otherwise layers would be eroded by turbulence.
One of the important features in a layered Lagrangian model is the treatment of the boundary between regions of different turbulence intensity [62] (cf. Equation (6)). For settling particles, Webster and Thomson [63] showed that, when particles are incident upon a thin region around the boundary, probabilities of transmission and reflection must be calculated. We have implemented the conceptualization of Thomson et al. [62] and Webster and Thomson [63] by numerical experimentation with heuristic, stochastic variation in thickness of a boundary region, wherein particles entering from the more turbulent side have a greater probability of being reflected back into the high-turbulence region. We chose outputs that conformed to the ideas and results of these previous workers, but detailed discussion of the different schemes is beyond the scope of the present contribution.

Similarity Theory
The velocity scales for turbulence in the z-direction and settling are, respectively, |w |, the absolute value of the time-mean vertical turbulence velocity, and w s . The ratio of these forms the simple dimensionless group, Π 1 : Likewise, the timescales for the processes can be used to examine the conditions under which diffusion or settling dominates. From Equation (3), the timescale of vertical diffusion, τ 1 , through a layer of depth, h, is given as τ 1 = h 2 κ . From Equation (6), the timescale of settling through the same layer, τ 2 , is τ 2 = h w s . The ratio of the two timescales indicates domination of particle transport through the layer by settling or turbulent diffusion in the vertical direction. The ratio is given by the dimensionless group, Π 2 : Potentially important or critical values of particle size and settling speed, w s , turbulence intensity, |w |, and eddy diffusivity, κ z are given in Table 3. In previous work describing the generation of ash cloud layering, values of κ have ranged from κ h = 10,000 m 2 /s and κ z = 10 m 2 /s [17] using the Puff VATD model, to 10 −5 ≤ κ z ≤ 9, with a mean of 1 m 2 /s, using the NAME model [27]. Note that in these models, larger-scale eddies are assumed to be characterized by the wind field as described in the NWP model, so that κ represents only the subgrid eddy diffusion. Thus, the effects of diffusion on eddies of the scale of the grid might be poorly resolved. Fero et al. [17] used NWP grids with a resolution of c. 200 km in the horizontal, and 1 km in the vertical, while Dacre et al. [27] used a grid with a resolution of c. 50 km in the horizontal, and 100 m in the vertical. In the Lagrangian model used in the present contribution, we have explored results for values of effective κ z ranging from ∼3 ×10 −6 to ∼6 ×10 3 m 2 /s, in constructing a three-dimensional, unsteady velocity and diffusivity field following the spectral technique discussed in detail in Jenkins [64] and Bursik [12]. To obtain turbulence layering in this spectral model, the Fourier series describing the vertical dispersion is truncated at different wavenumbers [64] in different layers, resulting in an alternation of high-and low-eddy diffusivity or turbulence intensity (Table 3). Table 3. Values of dimensionless groups for different particle sizes and layer thicknesses. Note that the boundary between Π 1 or Π 2 and unity tends to fall in 10-250 µm particle size range.

Numerical Analysis
To illustrate effects of turbulence layering on particle settling and motion, numerical experiments were performed using both the VATD model Ash3D [67] and a Lagrangian model (Section 2.2.2). For a thorough description of the Lagrangian model and the synthetic atmosphere in which it is run, see Jenkins [64] and Bursik [12]. Experimental parameter values of importance used in the models are shown in Table 4. In Section 3, we explore results from similarity analysis and numerical solutions based on the theoretical development. Numerical solutions are provided for atmospheres both layered and non-layered with respect to turbulence (Table 4).

Results
In this section, we first investigate expected behavior of ash particles given a layering in atmospheric turbulence structure using similarity theory. We then look at numerical results showing what the typical particle vertical speed should be, given a synthetic atmospheric turbulence layering (Table 4). Output from several simulations for a synthetic, turbulencelayered atmosphere using different particle size and settling speed are then shown to investigate the potential effects of turbulence layering on particles of different sizes. We then discuss possible turbulence layering in the atmosphere associated with the climactic 1991 Pinatubo eruption, and run simulations in the resulting synthetic atmosphere.
We gain insight into expected behavior using results of the similarity analysis, and refer to Equations (11) and (12) to explore asymptotic behavior. In layers for which Π 1 1 or Π 2 1, the mean turbulence speed is low relative to the settling speed, the timescale of diffusion is therefore long, hence motion is controlled by settling. In such a layer, mean concentration decreases linearly with time, as all particles exit the layer at approximately their settling speed. In layers for which Π 1 1 or Π 2 1, diffusion is rapid relative to settling, i.e., the timescale of settling is long, hence motion is dominated by diffusion and continuous suspension. For layers in which Π 1 ∼ 1 or Π 2 ∼ 1, Equation (6) holds, settling and turbulent transport are both important, and although particle concentration decreases exponentially, some particles persist in the layer a relatively long time due to turbulent suspension. Note that as layer thickness, h, increases, the timescale of diffusion increases faster than does that for settling, meaning it becomes more likely the particles will exit a layer by settling than by diffusion. For a typical diffusivity of κ = 1 m 2 /s in the UTLS [68] at 10 km altitude, and layer of depth h = 1 km, the critical settling speed, w s,crit , dividing settling from diffusion dominated motion is c. 1 m/s, which would correspond to a pumice particle of diameter c. 100 µm at about 2400 kg/cu m (e.g., [69]).
Numerical results for a layered system (Figure 3) are shown in Figure 4. Figure 3 is based on the observations of Cho et al. [40], who point out two layers of especially striking turbulence from 2-2.7 km and 3.8-4.2 km (Figure 3), but do not give specific values for turbulence intensity. We therefore apply high turbulence in these two layers through a simple Lagrangian random walk, and in other layers, low turbulence. The largest particles, in this case >1 mm, are dominated by settling and show little sensitivity to turbulence (Figure 4a). Intermediate sized particles, ∼10-250 µm (Figure 4b,c), in the more turbulent layers initiate a random walk, being "stuck" within the eddies and subjected to turbulencehindered settling. In these cases, consider the lower boundary of a layer with strong turbulence. All the particles there are subject to a random walk. They have a 50% percent probability of going up, and a 50% probability of going down. Those particles sent above the lower boundary due to turbulence are sent to a position higher above the boundary than their original position at the boundary. This will give them a greater chance to spend a longer time in the turbulent layer, whether or not one considers settling. Thus, for those layers dominated by the random walk and dispersion, Π 1 ≤ 1, Equation (6) holds, and the behavior seen in Figure 4b,c occurs. The motion of the smallest particles, here c. 1 µm, is dominated by turbulent diffusion. The particles spread slowly through the less turbulent layers, and rapidly in the more turbulent layers.
A of specific humidity squared q (z)] 2 /(Áz) 2 ; the BruntÀVäis (g/ " q)[q(z + Áz) À q(z)]/Áz, acceleration, q is the potenti mean potential temperature; number, Ri = N 2 /(dU/dz) 2 . thermodynamic effects of hu puted N v 2 and Ri v where the replaced by the virtual poten urated conditions and by the ature q e for saturated conditio [20] For statistics compari quantities, we interpolated a grid used above in order to h

Single-Parameter Prob Functions (PDFs)
[21] Figure 1 indicates the files used in this paper. One se were taken over water. [22] Figure 2d gives the tion (PDF) of the number of to latitude. The PDFs in Figu the free troposphere (solid li (dashed lines). One can see were taken in the lower to humidity PDFs show a ve troposphere (but with a nonand a very wet mode fo surprising since almost all th The broad tail in the free tr PDF should give us a sign which to examine the potenti generation. [23] Before we go on to look at an example profile. profiles of temperature (T), For vertical gradient quantit For the calculation Át = 0 lines at log = À3 indicate potentially unstable layers.   , high turbulent energy dissipation rate (log ε), and bounded by high total shear ((dv/dz) 2 ). Modified from Cho et al. [40]. Shaded layers are high-diffusivity, unshaded are low-diffusivity, used in simplified layer model (Table 4). (b) Rolling mean of absolute value of instantaneous vertical speed for an example model 30 µm particle (w S = 0.01m/s). Note that speed tends to be higher in more turbulent layers (shaded in (a)).
In contrast, and following from Equation (3), spread from a point source in a VATD model, with isotropic turbulence and a wind of constant speed with height, is shown in Figure 5, together with spread in a layered atmosphere from the Lagrangian model. Ash diffuses and progressively spreads from the source as the center of mass descends at the settling speed. Using larger particles and a higher settling speed, the rate at which the center of mass descends would increase, but the rate at which the particles disperse from the center of mass would not be affected by particle size. Thus, at any one height below the source, particles with a higher settling speed would be spread less distance from the source.   Figure 3. See Table 4. Most particles colored semitransparent blue, so deeper hue means more particles. A few separate particles are colored differently to show individual paths. In (b) and (c), the black particle path is for particle falling at effective settling speed. (a) 4 mm particles not affected by turbulence (Table 3) (Table 4). Instantaneous source and no wind shear, showing advection, settling and dispersal of 30µm ash after 30 h. Note in the Lagrangian layered model, particles are concentrated in lower, high-turbulence layer.

Application to 1991 Pinatubo Eruption
There are no direct observations of turbulence values associated with the climactic 1991 eruption of Pinatubo. As noted in Section 2.2.3, Fero et al. [17] explored values of κ z = 10 and κ h = 10 4 m 2 /s for the eruption. Given the understanding of turbulence generation, we can nevertheless speculate with some assurance that Rayleigh-Taylor cloud top turbulence generation was probably strong above the Pinatubo umbrella cloud [29], and that turbulence generation by both wind and umbrella flow near the tropopause at 16-17 km was probably strong due to the intrusion of the umbrella cloud and shift from southwesterly to northeasterly winds at that altitude [17]. This assessment is consistent with the general observation that the most intense level of clear-air turbulence development in the lower atmosphere is the tropopause [42]. We therefore construct a speculative turbulence layering for the eruption based on these considerations.
The main umbrella cloud of the climactic 1991 Pinatubo eruption was centered between 24 to 26 km height and ∼3-6 km thick, although ash was injected as low as 17 km and as high as 40 km [7,17]. The umbrella cloud acted as an intrusive gravity current for the first 4-5 h of eruption, and continued to be partly driven by buoyancy for an additional 5-6 h [7,17,46]. From the initial umbrella cloud, the development and separation of layers rapidly resulted in two main clouds, the first rich in larger particles (∼90 µm; [17]) and centered at the tropopause (∼16-17 km); the second, a higher cloud rich in SO 2 and fine particles that remained centered near the neutral buoyancy height of 25 km [17]. These two main ash transport regions persisted, centered around 14-16 and 22-25 km [8] until at least about a month after the eruption. Based on possible/likely generation of stronger turbulence by cloud-top convective and intrusion, and wind shear mechanisms, respectively, two regions centered around 17 and 25 km are given higher values of turbulence intensity in layered, Lagrangian model runs (Table 4).
With the atmosphere having these speculative intensely turbulent layers, we modeled the settling of 1 to 4000 µm particles from the Pinatubo cloud for periods up to a week. Although particles as large as 1000 µm were affected by the turbulence layering [70], we focus on the results for representative 10 and 30 µm particles at 18 h, as these sizes were dominant in parts of the cloud for periods up to nearly a month (Section 2.1.2) [8], as their settling and dispersal behavior were similar to that for particles <10 and >30 µm, respectively, and as the 18-h time window is particularly important to VAAC ash cloud forecasting. Eighteen hours is the length of the forecast window typically issued in volcanic ash advisories (VAA) and shown in volcanic ash graphics (VAG) issued by the VAACs. Also, by 18 h, the motion was dominated by the atmosphere [7,17]. The results for model 10 and 30 µm particle sizes are consistent with observations of cloud heights, particle sizes (Section 2.1.2) and speculative layering within them ( Figure 6). Given particles released from heights ranging from 23 to 27 km, after 18 h, 10 µm particles are concentrated in and around the 24-25 km high-turbulence layers, whereas 30 µm particles are mostly spread throughout the 14.5 to 19 km high-turbulence layer. This layer would not have persisted through time except for the continued suspension by the higher level of turbulence ( Figure 7). In other numerical experiments, particles <10 µm are distributed similarly to the 10 µm particles, while particles up to ∼250 µm are distributed similarly to the 30 µm particles, so these two sizes give a generally good picture of fine and coarse particle distributions at this time.   Table 4 and Figure 6 for Pinatubo. Note 18 h = 64,800 s. Existence of more turbulent layers enhances persistence of particles.

Discussion and Conclusions
In the present work, we have presented data and models on near-vent and distal volcanic cloud morphology and ash loading. We have performed numerical experiments comparing dispersal in an atmosphere with constant κ and variable κ z with height. The observational data suggest that distal clouds of depth O[0.1 to ∼5] km develop from nearvent clouds of O[1-10] km depth. The downwind clouds occur at heights consistent with the original eruption column heights for both tropospheric and stratospheric eruptions. The depth range of the distal layers, being generally less than the near-vent depth range, and the common stacking of distal layers, suggest that their development is controlled by atmospheric processes. The observations from the numerical test case are consistent with the working hypothesis that the layering of the atmosphere in turbulence intensity, generally causing alternating suspension and settling dominated behavior of appropriately sized particles (in the present experiments and as observed, those of O[10-100] µm) is a cause of distal layer morphology. Particles larger than those in this size range are dominated by settling and show little sensitivity to variations in turbulence intensity. The results for the Pinatubo test case show that the turbulence structure of the atmosphere can cause ash layer formation or increase ash layer persistence, through the process of enhanced suspension in vertically restricted regions of high turbulence. The numerical results for both the test and Pinatubo cases are consistent with and explain in more detail those of Dacre et al. [27] that the motion and layering of O[1-30] µm particle clouds is driven by a balance between wind shear and turbulent diffusion.
Because of the complex interaction of the values of turbulence velocity, layer thickness, Ozmidov scale, treatment of the boundary between regions of different turbulence intensity, and settling speed or particle size, the modeling needs to be carefully constructed. Errors in the present work could have arisen especially due to poor characterization of the turbulence field and boundary effects between the different layers. The results suggest that the particle size affected by the processes explored here is not only sensitive to model conditions but also to the true physical state of the atmosphere, thus were one to measure turbulence intensity in detail, it might be possible to invert that information to obtain particle sizes persisting in more turbulent layers.
Thus, in addition to the near-vent volcanic processes modifying ash cloud layering, atmospheric processes can produce distal multilayered ash clouds, due to multiple, alternating layers of turbulent and quiescent air. The distal ash layers scale to the depth of the alternating turbulent and quiescent layers in the atmosphere, which are O[0.1-1] km deep. The relevant model outputs (Figures 6 and 7) are not inconsistent with those of VATD and backtrajectory models, in which homogeneous turbulence can be associated with layers, given pre-existing ash cloud layering at the source vent, or wind shear [4,22,23,49,53]. It is likely that shear together with layering in turbulent diffusivity controls the thickness of distal layers for periods of days to perhaps a month by hindering settling of particles in the 1-1000 µm range [27]; note that the wind shear and turbulence layering processes are not mutually exclusive [40].
The results suggest that to better forecast the position and morphology of ash clouds for aviation safety and other purposes in VATDs, the vertical characteristics of the atmosphere need to be better resolved and characterized, particularly in the 18-h time window critical to forecasting in VAAs and VAGs. Because of the importance of turbulence and moisture to layer formation and detection, it is critical that these two parameters especially be estimated well, and at as high a vertical resolution as possible. Specifically, if appropriate sonde or other high-vertical resolution data are available, they should be used to measure heights of potentially elevated κ, to be used as input into high-resolution NWP models to provide better vertical turbulence parametrization for VATDs, and a better understanding of heights at which ash should occur in higher concentration. One perhaps important ramification of the present work is that it might be possible to determine a priori an ash-height injection time history, which should be critical to accurate forecasting of the height of ash, a priori, rather than find it a posteriori as currently done through inversion [5,24]. Being able to provide such information a priori could have a strong positive impact on forecasting the future heights and thicknesses of ash clouds for aviation using VATDs. Acknowledgments: The helpful comments of two anonymous reviewers on an earlier ms, and four anonymous reviewers on the current ms are gratefully acknowledged. Simon Carn is thanked for initial interest and encouragement in this work. Bruce Pitman is thanked for useful conversation on this topic. Larry Mastin gave valuable feedback on an earlier version of the ms regarding related work. Q. Yang was funded by the Earth Observatory of Singapore. We thank the creators of Ash3D, Pandas, Matplotlib and Seaborn for excellent software tools enabling this research.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: