Natural Hazards and Earth System Sciences

Abstract. The impact of ash-fall on people, buildings, crops, water resources, and infrastructure depends on several factors such as the thickness of the deposits, grain size distribution and others. Preparedness against tephra falls over large regions around an active volcano requires an understanding of all processes controlling those factors, and a working model capable of predicting at least some of them. However, the complexity of tephra dispersion and sedimentation makes the search of an integral solution an almost unapproachable problem in the absence of highly efficient computing facilities due to the large number of equations and unknown parameters that control the process. An alternative attempt is made here to address the problem of modeling the thickness of ash deposits as a primary impact factor that can be easily communicated to the public and decision-makers. We develop a semi-empirical inversion model to estimate the thickness of non-compacted deposits produced by an explosive eruption around a volcano in the distance range 4–150 km from the eruptive source. The model was elaborated from the analysis of the geometric distribution of deposit thickness of 14 world-wide well-documented eruptions. The model was initially developed to depict deposits of potential eruptions of Popocatepetl and Colima volcanoes in Mexico, but it can be applied to any volcano. It has been designed to provide planners and Civil Protection authorities of an accurate perception of the ash-fall deposit thickness that may be expected for different eruption scenarios. The model needs to be fed with a few easy-to-obtain parameters, namely, height of the eruptive column, duration of the explosive phase, and wind speed and direction, and its simplicity allows it to run in any platform, including a personal computers and even a notebook. The results may be represented as tables, two dimensional thickness-distance plots, or isopach maps using any available graphic interface. The model has been tested, with available data from some recent eruptions in Mexico, and permits to generate ash-fall deposit scenarios from new situations, or to recreate past situations, or to superimpose scenarios from eruptions of other volcanoes. The results may be displayed as thickness vs. distance plots, or as deposit-thickness scenarios superimposed on a regional map by means of a visual computer simulator based on a user-friendly built-in computer graphic interface.


Introduction
Dispersion and sedimentation of airborne tephra are the volcanic phenomena that affect the largest areas around erupting volcanoes.The impact of ash-fall on people, buildings, crops, water resources, and infrastructure depend on many factors, such as the duration of the fallout, thickness of the deposits, granulometric distribution and composition of the falling tephra and its leachates, among others.
When a volcano begins an episode of activity after a long period of quietness, authorities, media and the public usually do not have a clear perception of such an impact, and find it difficult to plan and prepare against ash falls.The complexity of tephra dispersion and sedimentation makes it difficult any attempt to look for a consistent perception of this phenomenon among different sectors of a population.On the other hand, preparedness against tephra falls over large regions around an active volcano requires an understanding of the eruptive process, and a working model capable of predicting most of its effects.
An attempt is made here to address one of those factors, namely the thickness of the ash deposits.A semi-empirical inversion model is used to estimate the thickness of noncompact deposits produced by an explosive eruption around a volcano in the distance range 4-150 km from the eruptive source.A graphic interface can be used to picture the expected deposits of a potential explosive eruption.The model was initially developed to depict deposits of a potential eruption of Popocatépetl volcano in central Mexico, and it was later applied to Colima volcano.This model seeks to assist planners and Civil Protection authorities, as well as other public sectors to obtain a more precise perception of what sort of ash-fall deposit thickness may be expected for different eruption scenarios.The model requires to be fed with a few, relatively easy to obtain parameters, namely the height of the eruptive column, the duration of the explosive phase, and the wind speed and direction.The computer graphic interface is user-friendly and has minimum CPU requirements.
The model was elaborated from the analysis of the ashfall deposits of 14 well-documented eruptions.In all of them, a basic set of data was available: H , height of the eruptive column; τ , duration of the explosive phase, and Ū , the wind vector (speed and direction).The model has been tested with available data from some of the recent eruptions Popocatépetl and Colima volcanoes.The model thus permits to generate ash-fall deposit scenarios from new situations, or to recreate past situations observed at any given volcano, or to superimpose scenarios from eruptions of other volcanoes.

Popocatépetl and Colima volcanoes
Popocatépetl, a 5452 m high volcano, is located in one of the most densely populated areas in the world.For over two millennia, its low eruption rate, the fertile volcanic soil, and the mild climate of the high plateau have favored population growth on hazardous areas around the volcano.Currently, over one hundred thousand people live in areas that may be vulnerable to the direct effects of a major eruption, and nearly 20 million live within a radius of 100 km, where they may be threatened by ash falls.Popocatépetl reawakened in 1993 with increased fumarolic and seismic activity, after nearly 70 years of quiescence.By October 1994, internal seismic activity increased significantly.This unrest culminated in the early hours of 21 December 1994, with a series of small explosions in the crater.Those explosions were followed by emissions of ash, which fell on several towns to the East and Northeast of the volcano and on the city of Puebla (Fig. 1).That night, the Mexican government evacuated about 75 000 people living in the most vulnerable towns.Afterwards, the activity decreased, and took the form of frequent short-duration ash emissions having a characteristic emerging seismic signature.These "exhalations" persisted with a decreasing rate through 1995.
In March 1996, a lava dome was observed growing on the floor of the crater.Since then, the activity turned more explosive.Several dome construction and destruction episodes have produced explosions of different intensities.On 30 June 1997, a relatively large eruption and a northwesterly wind produced a light ash-fall in Mexico City (with a metropolitan area population near 19 million) that shocked the public opinion.Some of the main urban regions that have been affected by ashfall from Popocatépetl are shown in Fig. 1.Although the nature and intensity of the current activity is similar to the previous episode (1919)(1920)(1921)(1922)(1923)(1924)(1925)(1926)(1927), this new activity and the eruptive history of Popocatépetl raised a generalized concern.
Geologic evidence sets the origin of the modern Popocatépetl after a cataclysmic eruption ca.23 000 years ago, which destroyed a previous volcano.Since then, many small eruptions and at least seven major eruptions have produced several cubic kilometers of ash and pumice.The most recent of the explosive eruptions occurred ca.3000 BC, 200 BC and 800 AD, as evidenced by archaeological remains buried by ash fall beds and pottery shards incorporated by mudflows (Siebe et al., 1996;De la Cruz-Reyna and Siebe, 1997;Siebe and Macías, 2004;De la Cruz-Reyna and Tilling, 2008).
Currently, the activity of the volcano seems to be declining, however, whether the current evolution will end without exceeding the levels of explosivity of the last few hundred years, or whether it may conduct to an eruption as large as the 1200 y.b.p. event, or even larger, is a question that cannot be answered in advance.
The last large eruption (VEI 4 ∼ 5) in 1913 left a deep 450 m diameter crater.Lava started to slowly fill that crater, reaching its rim around 1960.A succession of dome growth and destruction events in the summit crater followed during the last three decades, causing volcanic earthquake swarms, explosions and block-and-ash flows (Rodríguez-Elizarrarás et al., 1991;Núñez-Cornú et al., 1994;Saucedo et al., 2002;Bretón et al., 2002;Luhr et al., 2010;Saucedo et al., 2010).Unlike Popocatépetl, its high eruption rate has prevented large settlements around the volcano.However, about 300 000 live close enough to be seriously affected by ashfall.Such was the case of Ciudad Guzmán (85 000), 26 km NE of the volcano, covered by about 15 cm of ash during the 1913 eruption.Colima city (130 000), 31 km away, may be similarly affected in the event of southwards winds during a large eruption (Fig. 2).
Since large eruption cannot be ruled out in either volcano, the need of scenarios that may help planners and decisionmaking authorities to grasp the problem of ash-fall deposits and to assess the risks it poses is one of the main motivations of this study.Attempting to provide at least a partial answer to this request, a semi-empirical model for ash deposits is proposed here.
Particle-tracking are Eulerian (observer in a fixed location) or Lagrangian (observer moving along with the fluid parcels) models that can forecast volcanic-cloud positions at specific times.Lately, they have been used largely for aviation-safety purposes.Among them we can find ATHAM (Active Tracer High-resolution Atmospheric Model;Oberhuber et al., 1998); CANERM (CANadian Emergency Response Model used by the Canadian Meteorological Centre; D'Amours, 1998); MEDIA (Eulerian Model for DIspersion in the Atmosphere, used by the Toulouse VAAC); PUFF (particle tracking model used by the US National Weather Service, Anchorage, Alaska; Searcy, 1998); VAFTAD (Volcanic Ash Forecast Transport And Dispersion model developed by the NOAA Air Resources Laboratory; Heffter and Stunder, 1993); VOL-CALPUFF (a 3-D time-dependent model that couples an Eulerian approach for the plume rise with a Lagrangian representation of particle dispersal, Barsotti et al., 2008).Recently, a Lagrangian particle tracking module to model near-vent depositions, the Large Pyroclast Module or LPM has been incorporated to the ATHAM models (Kobs, 2009).
Advection-diffusion models are Eulerian representations of the solutions of particle diffusion, transport and sedimentation, and can forecast material accumulation on ground relative to a particle-release source.These models are mostly used for civil protection purposes, such as public warnings and the planning of mitigation measures.Among the latest are: ASHFALL (Hurst and Turner, 1999); FALL3D (Costa et al., 2006;Folch et al., 2009); HAZMAP (Armienti et al., 1988;Macedonio et al., 1988;Barberi et al., 1990;Macedonio et al., 2005); TEPHRA and TEPHRA2 (Connor et al., 2001;Bonadonna et al., 2005).These advectiondiffusion models are based on the analytical solutions of Suzuki et al., 1983.Applications to hazard zonation of different advection-diffusion models may be found for example in some maps computed by FALL3D simulating the 2001 Etna eruption, Italy (Costa et al., 2006) and the 2008 Chaitén eruption, Chile (Folch et al., 2008); in probability maps computed using HAZMAP for the tephra hazard assessment of Montserrat, West Indies; and Tarawera, New Zealand (Bonadonna et al., 2002(Bonadonna et al., , 2005)).Most of these models are physically complex and require sophisticated numerical methods derived from the complexity of the plume dynamics, the particle sedimentation, and the granulometric distribution of the deposits.
Not attempting to replace the physical models, the semiempirical inversion approach developed here searches for a simple-to use tool that may be used for training purposes, generating realistic scenarios for diverse eruption styles and intensities that can run in any platform, including small PC's and notebooks.Furthermore, this tool should be capable of rendering quick results, for specialists and non-specialists required to have an understanding of the ash-fall and its effects on the ground during an ongoing eruption.Thus, the model should be of simple use, yet have a degree of reliability be-yond its possible use as a training tool, assuring that no incorrect results may lead to wrong decisions of Civil Protection or Civil Defence officials.
Three possible lines were recognized in the process of developing ash-deposit thickness models: Numerical modeling, when fluid-dynamics equations and other phenomenological equations describing diffusion, sedimentation and meteorological aspects are solved with the appropriate boundary conditions in order to model the complex process of ash fall and deposit on physical grounds.These models usually require a significant computational power, and deposit thickness and granulometric distributions cannot be treated separately.See the website of the IAVCEI Tephra Group (http://www.ct.ingv.it/Progetti/Iavcei/index.htm) for a detailed review of these models.
Empirical models, which may be regarded as the opposite to the previous line.No physical arguments are used.A mathematical function describing an appropriate shape is fit to the field data looking for minimum errors.These are not strictly models but geometrical constructions for display purposes.
Semi-empirical modeling, that may be developed fitting mathematical functions resembling characteristic solutions of the physical processes involved.Such functions have free parameters that may be related to the eruptive and meteorological observables.Fitting is searched adjusting these parameters according to the reported eruption and wind data, and the field thickness data.
These lines are not necessarily independent or exclusive.In fact, empirical factors are present in all of them, although in different degrees.Even sophisticated physical models involving a full system of fluid-dynamics differential equations are data-based and require assumed values of phenomenological parameters that usually are difficult to measure (like diffusion coefficients or viscosity within or near the volcanic cloud for instance), in order to produce satisfactory results.

Semi-empirical modeling
In the semi-empirical approach discussed here, a basic assumption is implicit: a model constructed using data of different eruptions is capable to recreate scenarios of eruptions not used in the calculation of the model parameters, and thus is capable to generate new scenarios.The validity and limits of this assumption is discussed below.
To construct the model, we looked first into the isopach maps of 14 well-documented eruptions, covering most of the eruption-magnitude range.Those eruptions are shown in Table 1.
The isopach data base of each eruption was digitized as 3-D matrices of the ash-thicknesses as a function of the polar coordinates of the sites setting the origin at the eruption source and placing the X-axis coincident with the wind vector direction.
Table 1.Eruptions used as a reference to study the spatial distribution of ashfall deposits as a function of eruptive intensity, duration and wind conditions.

Volcano and date of eruption
The values at the intersections of the isopachs with a set of radial lines drawn every 10 • , generated a matrix T ij of deposit thickness measured at the azimuthal angle θ i , and at a distance r ij from the source, where i = 1, 36 and j depends on the number of isopachs considered (González-Mellado, 2000 -Appendix A).

Radial dependence
The model was constructed considering separately the radial and the azimuthal components of the T ij tensor.The radial behavior of the deposit thickness as a function of the source distance for each azimuthal angle rendered two functions which fit all the data well; (1) where T is the tephra fall deposit thickness, r is the radial distance to the source and A, α and δ are adjustable parameters.
Although the exponential function has been more frequently cited in the literature (Pyle, 1989;Sparks et al., 1991;Fierstein and Nathenson, 1992;Bonadonna et al., 1998;Legros, 2000) to describe the thinning of the deposits, the power law seems to better describe the thickness reduction with distance of freshly deposited tephra, particularly in the near and intermediate fields.To illustrate this, we present in Fig. 3a the closest 120 km of the ashfall (dry and uncompacted) deposit thickness measured three days after the 28 March 1982 eruption of El Chichón volcano along a northeastward line (see Table 2).
The best power fit with the fresh deposit data was: On the other hand, the best exponential fit for radial distances greater than about 8 km from the source for the fresh deposit data was: Figure 3a shows that the power law fits better, particularly in the near field.
In contrast, the isopachs obtained along a direction of 60 • counterclockwise from the wind axis in a further field study of the deposits of the 1982 eruption of El Chichón performed in June 1982 and January 1983 (see Table 2) shows that a slightly better fit is obtained with the exponential function.
T e (r) = 32.79e−0.106 r C r = −0.997, (5) although the power function correlation is still quite good, Figure 3b shows the fits of the stale deposit data with these functions.The inspection of best-fitting power function yields an interesting result: we found remarkable that functions for fresh and older, staler, deposits, (3) and ( 6) respectively, remained almost unchanged and that the value of the coefficient of the older deposits was about half the value obtained from the fresh deposits.
At this stage, we speculated that the exponent of the power function, describing the rate at which the thickness of the deposits decays with distance did not change much with compaction while the change of the coefficient reflects the compaction of the deposits after several months comprising a rainy season.Some measurements of the fresh to compact density ratio of the El Chichón ash deposits (non-compacted density/compacted density ∼ 0.4 − 0.5) also support this.It is thus likely that fresh deposits are better described by the power function, and compactation and perhaps some reworking make older deposits to be better described by the exponential function.This idea gained force after analyzing isopachs of deposits produced by the eruption of Popocatépetl volcano on 21 December 1994 and on 30 June 1997 and in Chaitén volcano, on 3 May 2008 (see Table 2).Again, a better fit was obtained with the power function in several radial directions, particularly in the near field, as shown in Figs. 4 and 5.

Azimuthal dependence
To clarify the two-dimensional functional dependence of the ground deposits thickness with distance from the source, the published isopachs from different eruptions were analyzed at different azimuthal angles in the way described in the previous section.The deposit thickness calculated from the intersection of radio vectors traced every 10 • with the isopachs, i.e., the values of the array T ij described above, were used to test the goodness of fit of potential and exponential distributions for different azimuthal angles.
The results for the fits for different azimuthal directions for the 21 December 1994 eruption of Popocatépetl and for the 28 March 1982 eruption of El Chichón are shown in Fig. 6a and b, respectively.Similar plots for the correlation coefficient of the exponential and power functions at different azimuthal angles were obtained for all other eruptions listed in Table 1.In most cases, the values of the correlation coefficients for both functions were high, in the range 0.9 to near 1.0, except in the directions opposite to the wind.
In about half of the studied cases, the power function fit was slightly better than the exponential function.In the other half, the goodness of fit was otherwise.Most of the cases in which the best fit corresponded to the power function were obtained shortly after the eruptions.In older deposits, the number of cases in which the best fit was with the exponential function almost doubled the number of cases in which the best fit was a power law.
We therefore choose the power function expressed by Eq. ( 2) to describe the radial decay of the deposit thickness based on the following grounds: 1.The fit of power law is good in both, near and far fields.
2. The model is intended to simulate deposit thickness from ongoing or days-old eruptions.The power function performs better in such cases.
3. A power-law (like that in Eq. 2) may also describe the ash-particle concentration in a volcanic cloud as a function of the distance to the source.The exact solution of the stationary diffusion equation for a continuously emitting point source varies as r −1 (Csanady, 1973).This solution would hold for a horizontally spreading cloud of fine suspended particles.However, in a real cloud, larger and heavier particles would cause loss of material by sedimentation, increasing the absolute value of the exponent of the power law.

Proposed model
In a static atmosphere, the lines of equal concentration of ash cloud particles would be concentric circles.Under these conditions, on a flat ground, we would expect the isopachs to reflect this circular distribution, as evidenced by a few circular deposits: Fogo-A in Agua de Pau volcano, Azores (Bursik et al., 1992) and BF in Pululagua volcano, Ecuador (Papale and Rossi, 1993;Volentik et al., 2010).Similarly, if the ash cloud particle concentration radially decays according to a power law, we may expect that the thickness of the circular isopachs would also decay according to a power law, with different parameters though.
If the medium in which the diffusion process develops is moving, i.e., in an eruption injecting relatively fine tephra into the atmosphere in the presence of wind, the diffusion process may be perceived as in the static case, but referred to a coordinate system moving with the wind.This means that the geometry of the equal-concentration curves respect to the ground would be changed according to a coordinate transformation of the type x i = x i − U t, where U is the wind speed.In this case the lines of equal concentration of ash cloud particles would elongate along the wind direction taking a quasi-elliptic shape.We thus assume that the shape of the isopachs in the ground will deform in the same fashion.
Sustaining the validity of a uniform wind field assumption is not a simple matter.However, some considerations to this respect may be briefly addressed here.The "wind speed" U above does not really intends to represent an unrealistic constant and uniform wind field, but rather a mean value of the predominant wind near the eruption source.Random fluctuations around this value should average out over the duration of major eruptions, particularly in the "near" field, in which most of the eruptions analyzed here have produced essentially elliptic isopachs.Far field deposits, dominated by highly wind-sensitive fine grains may be strongly distorted by wind field fluctuations.Since the model proposed here is aimed to generate ashfall scenarios within a limited region around the erupting center, we assume U to be a constant.
Under such a consideration, the following semi-empirical function is proposed to describe the deposit-thickness field around an eruptive center: where α and γ are constants depending of the eruption parameters, and F (r,θ,U ) is a shape function derived from the above transformation controlling the ash-deposit thinning with distance r, angle between the wind direction and the vector from the ash emission center to the test point θ, and wind velocity U .Csanady (1973) found the solution for the diffusion equation for a gas continuously emitted by a point source in a fluid moving at constant speed U .The gas concentration χ at a distance r, and at an angle θ with the wind vector is: where D a is the molecular diffusivity and I is the ash emission rate at the source.Based on this, we constructed the function F (r,θ,U ) on semi-empirical grounds using a general form of the azimuthal component: where β is a parameter depending inversely from the diffusion coefficient, and f (r,θ) is a purely geometrical function governing the quasi-elliptic shape of the equal concentration curves in the ash cloud It is important to emphasize that we are assuming that the geometry of the isopachs are geometrically similar to the lines of equal concentration in the ash cloud.
The final form of the semi-empirical function describing the thickness distribution of fall deposits is then: where parameters α is dimensionless, and β and γ have units of [TL −2 ] and [L/L −α ], respectively.Preliminary values of the model parameters α, β and γ , were estimated from characteristic or typical values of the different eruption parameters and diffusion coefficients.Those values were then adjusted to fit the sets of thickness data T ij obtained from the 14 reported isopachs.The method of adjustment was a recursive search for minimum variances, calculated as the squared differences between the data base T ij and the model predictions.Table 3 shows the eruption parameters for each event, and the best estimates of the model parameters (González-Mellado, 2000 -Appendix B).
In this study, we assume that α, β, the model parameters acting as arguments of the radial and azimuthal distributions, are quantities depending on the observable eruption parameters governing the dispersion of volcanic ash, namely the eruptive column height H , and the duration τ of the eruptive phase producing that column.The dependence of the deposit thickness distribution with the wind velocity is explicit in Eq. ( 11), and thus assumes that the model parameters do not depend on U .Intuitively, one may also expect a dependence on the density of the falling material.However, for recently deposited ash it is very weak, since in most cases the density of the uncompressed fresh deposit was very close to 1100 kg/m 3 .
The next step of the model development was the search of the explicit form of the relations: The structure of the parameter γ is discussed separately.To determine the above functions, we made a simple analysis of the correlation between each model parameter and the eruption parameters.Analysis of the power law exponent α, which determines the rate at which the deposit thickness decays with distance, yields clear results.An almost linear dependence with the column height H results. Figure 7a shows the line α(H ) = 2.283−0.041Hobtained from the regression of the isopach and column height data of the twelve eruptions indicated in the figure caption.Inspection of that plot shows that two points, corresponding to the Vesuvius, Italy, 22 March 1944, andFuego, Guatemala, 14 September 1971 depart significantly from the line, reducing the value of the correlation coefficient to −0.885.This may be caused by the "twisting" of the farthermost isopachs observed in those eruptions, an effect probably derived from changing wind regimes in the far field.A new fitting ignoring Vesuvius and Fuego is shown in Fig. 7b, where the regression line now has a very satisfactory correlation coefficient −0.981.
In these equations H is measured in km.
On the other hand, attempts to correlate the power law decay parameter α with the duration of eruptions yields very small slopes and low correlation coefficients.Even smaller slopes and correlation coefficients result form the attempts to correlate α with the wind velocity U .We can therefore conclude that the potential decay parameter α depends only on the height of the eruptive column, which in turn is a nonlinear function of the mass rate and power output of the eruptions (Settle, 1978;Wilson et al., 1978;Mastin et al., 2009).
A similar analysis was performed searching for correlations for the "effective" diffusion coefficient D. We call D "effective" because it is an empirical parameter calculated from the relationships between the emission rate (represented www.nat-hazards-earth-syst-sci.net/10/2241/2010/Nat.Hazards Earth Syst.Sci., 10, 2241-2257, 2010  by the column height) and the isopach distributions.It thus plays the role, and has the units of diffusivity, but it is not necessarily a measure of the actual diffusivity coefficient.
Comparing the arguments of the exponential functions in the solution of Csanady (1973) given by Eq. ( 8) with the model given by Eq. ( 11), we obtain an inverse relationship between parameter β and effective diffusion coefficient D: If this was the case of a gas or a finely dispersed contaminant diffusing into the atmosphere, the parameter D would be just the diffusion coefficient of the gas.However, when an ash cloud spreads it is not evident that relationship (15) holds.
We thus assume its validity as an additional hypothesis of this model, and use it to evaluate D, using the values of β listed in Table 3.
The points in Fig. 8 show the values of D calculated from Eq. ( 15), and listed in Table 3, as a function of the column height H .A simple inspection of the distribution of the points indicates that two distinct regions showing opposite trends may be recognized.Points align with a negative slope below about 15 km.Above that, the points are more scattered, but with a clear positive trend.This behavior is consistent with the thermal structure of the atmosphere.In the troposphere, temperature decreases with height until an altitude that depends on the latitude.Above that, in the stratosphere the vertical temperature gradient inverts.The height of the border region, the tropopause, varies from about 8 or 9 km in the pole to about 17 or 18 km near the equator, depending on the season.Considering that the Einstein-Smoluchowski relation of kinetic theory, the actual diffusion coefficient varies linearly with the temperature, (D a = bkT , where b is the mobility, k is the Boltzmann's constant and T is the cloud absolute temperature), the trends in Fig. 8 indeed reflects the way in which the thermal structure of the atmosphere influences the particle diffusion within the volcanic cloud.
The volcanoes used for this study lie in latitudes lower than about 46 • , and the altitude of the slope change is thus representative of the troposphere-stratosphere boundary.In particular, the tropopause above the central part of Mexico averages an altitude of 16.5 km (Cortes-Luna, 1996).
From these considerations, we adjusted two lines to the points of Fig. 8 setting the intersection of the lines at a mean tropopause height H = 15.5 km.
The best fitting lines result as follows.If the column is tropospheric, D(H )=−4.189H+114.407,0<H <15.5 (C r =−0.992) , (16a) and for an eruptive column penetrating into the stratosphere, D(H )=52.822H−770.17, 15.5≤H ≤50 (C r =0.920) , (16b)   where H is measured in km and the units of D(H ) are km 2 /h.The value of D for the 19 April 1993 eruption of Lascar was omitted in this fit because the isopachs of these deposits lack sampling points in the direction perpendicular to the prevailing wind direction U .
To search for an explicit form of the coefficient γ in Eq. ( 11), we used dimensional analysis.First, we assumed that γ had a functional dependence on the main parameters of the eruption, namely, the eruption intensity I (mass eruption rate), the time τ in which such intensity holds, the density of the erupted material ρ, and the model parameters α and β.Equation (11) should thus satisfy the dimensional equation: where M, L and T are mass, distance and time unit, respectively.The deposit thickness T hick should have dimensions of distance [L] only, and therefore, the product of γ by the power law function must have the same dimensions (since the exponential function is dimensionless).A solution was found solving the dimensional equation where a, b, c, and d are the exponents of the variables such that both sides of the equation have the same dimensions.This generates the system of equations: Since this system has multiple solutions we choose the simplest one, a = 1.The solutions for the other variables thus are Thus, an expression for γ is: where A is a dimensionless constant.The intensity of eruptions and the height of the eruptive columns are related by the fourth-power law (Settle, 1978;Wilson et al., 1978;Mastin et al., 2009), where B is an empirical constant.Inserting (23) into Eq.( 22), we obtain, Replacing the values of the parameters of Table 3 in (24), the resulting set of values for the product AB is shown in Table 4.However, the best value of the coefficient of AB to fit the observed isopachs is 17.5% larger than the mean value obtained from Table 4.We thus have that the best value for AB is 1.175 • 12.15 = 14.28 (10 −11 /3.6) km −4 s −1 kg.A final expression for Eq. ( 11) is then, where α and D are given by ( 14), (17a) and (17b), respectively, and variables units are

SECCVO: a visual computer tool to simulate ash fall scenarios
Equation (25) yields a deposit thickness value for a pair of polar coordinates (r, θ) and a given set of parameters.For planners, Civil Protection authorities, and for teaching purposes it is more useful to display the expected deposit thickness for an area.This can be accomplished using a standard graphics package such as Matlab, Mathematica or a GIS, overlaying the ashfall model thickness values on a specified map.
A displaying package called SECCVO (Spanish acronym for Ash-fall Scenario Simulator) has been used for training Civil Protection personnel in volcanic hazard assessment.SECCVO has proved to be a useful application of model.A built-in options allows the reproduction of fallout deposits of important historical eruptions or to recreate possible scenarios based on user-provided eruption and meteorological data.Figure 9 depicts the expected ashfall deposit produced by a hypothetic eruption of Colima volcano similar to the 28 March 1982 eruption of El Chichón volcano, with prevailing winds directed to Colima City.Such a scenario estimates ash deposits about 4 cm thick over that city.
Figure 10 depicts the expected ashfall deposit produced by the same 28 March 1982 eruption of El Chichón volcano hypothetically occurring in Popocatépetl, with prevailing winds directed towards Mexico City.The model estimates 1.13 cm of ash deposited in the Benito Juarez International Airport of México City (AICM), 0.42 cm in Cuernavaca City, and 0.157 cm in Puebla City.This type of modeling helps planners and Civil Protection officers to visualize the impact of ashfall and to identify vulnerabilities (De la Cruz-Reyna, 2001).
Particularly, on 30 June 1997 at 19:26 LT, Popocatépetl volcano, produced a moderately large eruption causing some ashfall on the AICM, located 66 km to the NW of the volcano.Aircraft flying near the volcano reported that an ash column reached a maximum altitude of 13 000 m a.s.l.(about 8 km above the volcano summit) with predominant Northwesterly winds at that height (SCT-Mexico, 1997).Wind speed was estimated as 80 km h −1 by the National Meteorological Service.The Navigation Service in the Mexican Airspace (SENEAM, personal communication) reported the following day that the ashfall on the runways of AICM amounted to 70 gr m −2 , corresponding to a deposit thickness of about 0.006 cm assuming a density of 1100 kg m −3 the uncompacted ash.The airport was closed to all operations for a period of 8 h 30 min, beginning on 30 June 1997 until the ash cleanup tasks were finished (SCT-Mexico, 1997).Figure 11 shows a SECCVO ashfall scenario for the Popocatépetl eruption of 30 June 1997 with column height H = 8 km above the volcano crater, duration τ = 0.58 h and wind speed U = 80 km h −1 toward WNW (azimuthal angle ∼ 312 • N).The isopachs obtained by Martin Del Pozzo et al. (2008) for the same eruption have been added to the scenario.The airport is located between the 0.005 and 0.01 cm isopachs, in good agreement with the model result of 0.005 cm.2008).The AICM is located between the first and two outermost isopachs (0.005 and 0.01 cm).

Discussion and conclusions
The thickness of a non-compacted deposit is closely related to the concentration of ash in the eruptive cloud.While the concentration of particles in a cloud expanding by diffusion processes is governed by a r −1 law, the thickness of the freshly deposited ash layer follows a r −α power law.The effect of wind is introduced as a change of coordinates into a coordinate system moving with the wind, generating an elongation of the concentric circles of decreasing concentration produced by a windless diffusion process.These geometric factors are described by Eqs. ( 9) and (10).The coefficients and constants were estimated adjusting their values to fit available data from 14 well documented eruptions, shown in Table 1.
Estimating the ground distribution of deposited fallen-ash thickness, function (25) requires only basic eruption information: eruptive column height, duration of the eruptive phase that produced the column, wind velocity and direction.The diffusion-related parameters must be changed for volcanic cloud altitudes below or above about 15 ∼ 16 km according to Eqs. (17a) and (17b).
Ash cloud dispersion, ashfall deposition and other factors of the eruptive cloud dynamics depend on the grain-size distribution in the volcanic cloud and its interactions with the wind, as emphasized by Bursik et al. (1992), Sparks et al. (1997), andScollo et al. (2008) among others.The present model does not consider the grain-size distribution factors, since no information on them is explicitly incorporated in the inversion of the isopach data.The model just provides an estimate of the thickness-distance relationship that is simple to implement in any computer, and is intended to display ashfall scenarios that are sufficiently realistic to be considered in Civil Protection decision-making during situation of volcanic crisis.

Fig. 1 .
Fig. 1.Map of Popocatépetl volcano and its densely populated surroundings.Urban settlements within a radius of 100 km around Popocatépetl volcano (shown in yellow) may exceed 20 million.

Fig. 2 .
Fig. 2. Map of Colima volcano and surrounding area.Colima City and Ciudad Guzmán are the largest urban settlements, which along other smaller cities (shown as non-labeled yellow spots) may be seriously affected by ashfall.Guadalajara City (4 000 000), about 160 km North of Colima volcano may also be affected by lighter ashfalls.

Fig. 3 .
Fig. 3. Ashfall deposit thickness vs. distance data for the 28 March 1982 El Chichón eruption.(a) Best fittings along a radial line at 60 • (counterclockwise respect to the wind axis) of fresh ash deposits measured a few days after the 28 March 1982 eruption (diamonds) with the power function (Eq. 3) (solid line, C r =-0.995), and with the exponential function (Eq.4) (dotted line, C r = −0.845).(b) Thickness and best fits along a radial line at 60 • (counterclockwise respect to the wind axis), of the same deposits measured several months later (diamonds) with the exponential function (Eq.5) (dotted line, C r = −0.997)and with the power function (Eq.6) (solid line, C r = −0.950).All of the data in Table 2 are used for the inversions, but only the closest (120 km) data are shown in the figure for clarity.

Fig. 4 .
Fig. 4. Ashfall deposit thickness data for the 21 December 1994 Popocatépetl eruption.Field data (diamonds obtained from the isopachs of Martin Del Pozzo, 1995), and best fittings to the power function (solid line), and to the exponential function (dashed line).(a) Thickness values and best fits along a radial line at 290 • (counterclockwise) respect to the wind axis (solid line, C r = −0.983;dashed line C r = −0.973).(b) Same as in (a), but along a radial line at 320 • respect to the wind axis (solid line C r = −0.981;dashed line C r = −0.977).

Fig. 5 .
Fig. 5. (a) Ashfall deposit thickness data for the 30 June 1997 Popocatépetl eruption along a radial line at 0 • , the wind axis.Field data (diamonds obtained from the isopachs of Martin Del Pozzo et al., 2008), and best fittings to a power function (solid line, C r = −0.989)and to an exponential function (dashed line, C r = −0.960).(b) Ashfall deposit thickness data for the 3 May 2008 Chaitén eruption (Chile), along a radial line directed at 0 • (wind axis).Field data (diamonds obtained from the isopachs of Watt et al., 2009) and best fits to a power function (solid line, C r = −0.988)and to an exponential function (dashed line, C r = −0.974).

Fig. 6 .
Fig. 6.Correlation C r as a function of θ (measured counterclockwise respect to the wind axis) for the fits in different directions for (a) the 21 December 1994 eruption of Popocatépetl and (b) the 28 March 1982 eruption of El Chichón.The "x" are the correlations to the power function, and the "o" to the exponential function.

Fig. 9 .
Fig. 9. SECCVO ashfall scenario of an eruption similar to the 28 March 1982 eruption of El Chichón volcano occurring at Colima volcano, assuming a prevailing wind blowing towards Colima city (N 200 • SW).This scenario displays the topography and the main towns and villages.The simulator estimates an ash thickness of 4 cm in Colima City, 8.6 cm in Suchitlán, 3 cm in Alzada and 9.2 cm in San Jose del Carmen (Jalisco State).The program displays the ash thickness in any point selected with the cursor (Colima city in this example).

Fig. 10 .
Fig. 10.SECCVO ashfall scenario of an eruption similar to the 28 March 1982 eruption of El Chichón volcano occurring at Popocatépetl volcano, assuming a prevailing wind blowing towards Mexico city (N 315• NW).The scenario displays the Federal District of Mexico, and other main towns and villages.The simulator estimates an ash deposit 1.30 cm thick in "Central de Abasto" (main food distribution center of Mexico City), 0.97 cm in the National Autonomous University of Mexico main campus, 0.42 cm in Cuernavaca City, 10.55 cm in Amecameca, and 0.157 cm in Puebla City.The displayed Information window shows a deposit 1.13 cm thick at 66 km to the NW of the volcano, the position of the Mexico City International Airport marked with the cursor.

Fig. 11 .
Fig. 11.SECCVO ashfall scenario of the 30 June 1997 Popocatépetl eruption.The eruptive parameters were: column height H = 8 km, eruption duration τ = 0.58 h (32 min), wind velocity U = 80 km/h (blowing towards Mexico city: ∼312 • NW).The Information Window displays 0.005 cm ash-thickness at the International Airport of México City (where the cursor is) at 66.00 km from volcano.The isopachs (in green) are from Martin Del Pozzo et al. (2008).The AICM is located between the first and two outermost isopachs (0.005 and 0.01 cm).

Table 2 .
Distance (r) and ash thickness (T ) data for some recent eruptions.

Table 3 .
Eruptive parameters.Column height H , wind velocity U , duration of the high-intensity phase of the eruption τ , and the semiempiral model parameters α, β, and γ of the 14 well-documented eruptions.

Table 4 .
Value of (AB) −1 and AB for different eruptions.In all cases a density of 1100 kg m −3 has been used for the uncompacted ash deposit.