Natural Hazards and Earth System Sciences Volcanic ash forecast – application to the May 2008 Chaitén eruption

We model the transport and subsequent deposition of ash from Chait́ en volcano, Chile, during the first week of May 2008. The simulation couples the Weather Research and Forecasting (WRF) meteorological model with the FALL3D dispersion model. We only use semi-quantitative volcanological inputs based on the first eruption reports. We consider two different run types based on forecasted and hindcasted meteorological conditions. The first simulation type can be regarded as a syn-eruptive operational forecast for the 2– 8 May period. We predict the evolution of the ash cloud position, the concentration of ash on air, the expected deposit thickness, and the ash accumulation rates at different localities. The comparison of model results with observed cloud arrival times and satellite images shows the goodness of the combined WRF+FALL3D forecast system and points out the feasibility of combining these two models for shortterm forecast of volcanic clouds and ash fallout.


Introduction
Fallout of tephra following explosive volcanic eruptions threats the communities located around active volcanoes and beyond.The deposition of lapilli and ash has many negative effects on human settlements and infrastructures at short to medium distances from the source.On the other hand, large quantities of fine ash and aerosols injected into the atmosphere can remain air-borne from days to months and be transported by high level winds thousands of kilometers downwind.Such fine particles affect the air quality and threat the aerial navigation causing drifting of aircraft routes to prevent degraded engine performance, loss of visibility, and possible failure of navigational instruments.Modeling and fore-Correspondence to: A. Folch (arnau.folch@bsc.es)cast of volcanic ash clouds is, consequently, a relevant issue for civil protection agencies and aerial aviation organizations (mainly concerned with the atmospheric ash concentration and ground ash deposition at airports).
In the light of the above, the International Civil Aviation Organization (ICAO) promoted the creation of the Volcanic Ash Advisory Centers (VAACs) to keep pilots and international aviation industry informed of the location and movement of eventual volcanic ash clouds.Whenever there is an on-going eruption, the regional VAAC makes use of direct observations by pilots, satellite imagery, and modeling, to issue warning bulletins and to forecast the short-term evolution of the cloud.Dispersion models used by VAACs run at a regional scale and use both Lagrangian, like HYSPLIT (Draxler and Hess, 1998), NAME (Ryall and Mayron, 1998) or PUFF (Searcy et al., 1998), and Eulerian approaches, like CANERM (D'Amours, 1998) or MEDIA (Sandu et al., 2003).One advantage of Lagrangian particle tracking models is their capacity to predict an approximate ash cloud trajectory with low computational effort, a key aspect during emergency situations.However, with few exceptions (e.g.Tanaka and Yamamoto, 2002), such models do not provide quantitative predictions for neither atmospheric ash concentration or ground accumulation, a drawback that makes them less attractive to civil protection authorities.In this paper we use the recent Chaitén eruption, Chile, to test the operational capacity of the FALL3D Eulerian model (Costa et al., 2006;Folch et al., 2008) to forecast ash concentration on air and expected deposit load and thickness.We run a one-week simulation coupling the Weather Research and Forecasting (WRF) meteorological model (Michalakes et al., 2005) and the parallel version of the FALL3D code.The model runs as if responding to the eruption, i.e. we only use as volcanological inputs the semi-quantitative data available at that time.It is important to point out that our goal is not to perform a detailed study of the Chaitén eruption but to use it as a blind test to confront short-term model predictions and semi- quantitative syn-eruptive observations.Actually, we remark that, at the time of writing this article, the Chaitén eruption was still going on and no quantitative field data was available yet.

The May 2008 Chaitén eruption
Chaitén volcano (Fig. 1), a caldera volcanic complex filled with an Holocene rhyolitic lava dome (Naranjo and Stern, 2004), reawakened unexpectedly on 2 May 2008 after a long period of quiescence.During the last days of April, the Chilean "Servicio Nacional de Geología y Minería" (SER-NAGEOMIN) (www.sernageomin.cl)and the Chilean National Emergency Office (ONEMI) (www.onemi.cl)raised the level of alarm after recording anomalous seismic activity in the area.The eruptive process started on 2 May early morning, and was characterized by a vigorous plume that rose between 10.7 and 16.8 km into the atmosphere according to the Smithsonian's Global Volcanism Program report (www.volcano.si.edu).As often occurs in South Andean eruptions, the dominant regional winds directed the ash cloud over the Andes and caused fallout in the Argentinean Patagonia (a part, of course, of producing an intense precipitation at the volcano proximal areas in Chile).Here we give a brief chronology of events occurred during the first week of eruptive activity in order to compare them with the predictions of the model (see the ONEMI webpage and the local press agencies for details).
-2 May (Friday -4 May (Sunday).The ash cloud directs towards SE.Copious fallout at the Chilean town of Futaleufú (located at ∼75 km from the volcano), where a deposition from 10 to 30 cm of ash is reported.
-5 May (Monday).The Chilean authorities evacuate Futaleufú.In Argentina, precipitation of small quantities of fine ash at Trelew, Rawson and Puerto Madryn.Shut down of several Argentinean regional airports due to the lack of visibility.
-6 May (Tuesday).Increase of the eruptive activity.ONEMI reports two big energetic bursts (one at 09:30 and another at 16:00 LT) which temporally enlarge the eruptive column up to 30 km in height (see Fig. 2).Puerto Madryn and Trelew airports do not operate.Fallout at El Bolsón (135 km from the vent), where particle sizes up to 0.5 mm are reported.During the afternoon, the ash cloud starts to drift towards NE across the Argentinean Río Negro province to cover, at around midnight, the South of neighbor province of Neuquén.
-7 May (Wednesday).Light fallout reported at Bahía Blanca and Mar del Plata, located at more than 1000 km away from the volcano.The eruptive column stabilizes at 12-13 km height.
-8 May (Thursday).The presence of ash at medium to high atmospheric levels (from 3 to 10 km) is reported at Buenos Aires.Several national and international flights are cancelled in both airports of the Argentinean capital.
3 Modeling strategy

Meteorological modeling
We used the Weather Research and Forecasting (WRF) model (Michalakes et al., 2005) to simulate the meteorological conditions from 2 May 2008 at 00:00 UTC to 9 May 2008 at 00:00 UTC.Note that, in order to warm up the meteorological model, there is a time-lag of 10 h between the beginning of the meteorological simulation and the beginning of the eruption, assumed to start at 10:00 UTC (06:00 LT).WRF is a fully compressible, Eulerian non-hydrostatic mesoscale model that solves the equations of atmospheric motion.The model was configured to integrate the primitive equations using the ARW dynamics solver (Skamarock and Klemp, 2008).The following physical parameterization was considered: the single-moment 3-class microphysics scheme (Hong et al., 2004), the Kain-Fritsch cumulus scheme (Kain and Fritsch, 1990), the Noah Land-Surface model (Skamarock et al., 2005), the Yonsei University PBL scheme (Noh et al., 2003), the RRTM long wave radiative model (Mlawer et al., 1997), and the short-wave radiative based on Dudhia (1989).A high-resolution domain which covers most of Chile and Argentina was defined (Fig. 1).The domain was centered over Argentina and discretized using 200×200 grid points at 12 km of horizontal resolution with 38 σ -vertical levels (11 to characterize the boundary layer).The model pressure top was set to 50 hPa (∼20 km).We considered two types of runs based on forecasted and hindcasted meteorological variables.The forecast run is configured using initial and 6hourly boundary conditions for the WRF mesoscale model from the high-resolution (0.5 • ×0.5 • ) Global Forecast System (GFS) results.The GFS is the global forecast model of the National Centers for Environmental Prediction (NCEP) that runs 4 times per day and produces forecasts for the next 180 h.GFS uses meteorological analysis as initial conditions for its forecasts.Meteorological analysis are a grid-base description of the state of the atmosphere at a given time obtained from several observation sources.The hindcast run is configured using as initial and boundary conditions the 6-hourly meteorological analysis of the GFS system, which should produce more accurate results than the forecast run.In the first case the FALL3D model runs exactly as if operational and the simulated period is split in 3-days intervals (that is, we run the 2-5 May interval using the 72 h forecast available on 2 May, then restart the run for the 5-8 May period with the data available on 5 May, and so on).The second case serves to describe the dispersion episode and as a reference case to validate the accuracy of the forecast run.

Ash transport and deposition modeling
The ash cloud forecast was done using FALL3D, a 3-D timedependent advection-diffusion-sedimentation model for the atmospheric transport and dry deposition of particles.We generalized the original model to include curvilinear coordinates.It allows to consider large-extension domains in which the Earth's curvature is not negligible (actually, it was a major limitation of the previous version of the code).The model inputs are meteorological data, topography, grain-size distribution, shape and density of ash particles, and Mass Eruption Rate (MER).
The computational domain spans in longitude from 58 E to 74 E (17 meridians), in latitude from 34 S to 48 S (15 parallels), and has 17 km in altitude.The spatial discretization contains 251×251×18 points, corresponding to a grid spacing of 1 km in the vertical and, depending on the latitude, between ∼5 and ∼7 km in the horizontal.The simulation covers 158 h, from 2 May at 10:00 UTC to 9 May at 00:00 UTC.The WRF meteorological variables were hourly interpolated from the WRF mesh to the FALL3D mesh.
We consider a Gaussian granulometric distribution with 10 particle classes ranging in size from =−4 (16 mm) to =5 (31 µm), peaked at m =0 (1 mm), and with a deviation σ =2.5.These values are representative of sub-plinian eruptions.For example, Bonadonna and Houghton (2005) find σ =2.4 and m =−1 for the 1996 Ruapehu eruption.However, we point out that this "analog" granulometry can actually differ from the real one, for which there were no measurements at the time of running the simulations.We also assume that particles have a sphericity of 0.9 and a sizedependent density which varies linearly from 1600 kg m −3 for =−4 to 2600 kg m −3 for =5.The source term is given by the 1-D radially-averaged Buoyant Plume Theory (BPT) (Bursik, 2001).In order to reproduce the observed column height (∼14 km on average), the BPT requires a MER of 3×10 7 kg s −1 (we assume an outflow velocity of 120 m s −1 and a mixture temperature of 850 • C, a characteristic value for a rhyolite).For simplicity, we consider the MER to be constant during the whole simulation (note that, at the time of forecasting, one has no information about the future evolution of the plume height).Note also that a constant MER does not imply a constant column height because the time-varying wind intensities produce different plume bent-over (see Fig. 3).This MER gives a total erupted mass of 1.7×10 13 kg for the 2-8 May period.The horizontal diffusion is internally computed by FALL3D using the approach of the RAMS model (Pielke et al., 1992).The diffusion coefficient depends on the grid scale and, for the grid size considered here, has values in the range of ∼30×10 3 m 2 s −1 .Such values are substantially larger than the 5-10×10 3 m 2 s −1 typically found when solving an inverse problem (inverse problems are normally solved in smaller domains, from short to medium distances, and hence typically involve finer grid resolutions, in the range of 1-5 km).

Meteorological results
The synoptic simulation of the episode under study is summarized in Fig. 4. The meteorological situation on 2 May 2008 is characterized by a region of high pressures over the Atlantic ocean extending towards Argentina and a cut-off low south-west of southern Chile.The Chaitén region is affected by a ridge in the upper troposphere with moderate to weak SW synoptic winds (Fig. 4a) and important wind shear in the upper troposphere (strong vertical change in direction and speed).The synoptic evolution from 2 May to 4 May is characterized by the extension of high pressures over central Patagonia leading to a large region with weaker winds (below 5 m/s) located slightly north of the Chaitén latitude.This configuration leads to strong advective NW motions above latitude 45 S and from moderate to calm winds at lower latitudes (Fig. 4b).The subsequent displacement of high pressures towards the Atlantic ocean allows the penetration of moderate W flows that dominate during the next two days (Fig. 4c).From 7 May on, the high pressure system located over the Pacific Ocean (latitude 36-39 S) advances towards the Chilean coast causing a SW veering of synoptic winds.Low to moderate winds (10-20 m/s) blow over North Patagonia and Pampa (Fig. 4d).

Ash forecast
The simulation predicts an initial ash cloud movement towards NE during the first hours of eruptive activity.However, a rapid SE direction drift directs the cloud over Argentina to reach the Atlantic coast at the latitude of Comodoro Rivadavia during the first hours of 3 May.During 4 May morning and afternoon the ash cloud keeps quite steady in the ESE direction (passing over Futaleufú) and then shifts eastwards during the evening.Thus, during the morning of 5 May, the direction is exactly E affecting Esquel and, to a lesser extent, the North of the Chubut province coast (Trelew, Mardyn and Rawson).During the evening the ash cloud directs again towards SE for few hours only because a rapid change in wind direction during 6 May drifts the plume NE.Due to the low intensity winds (see Fig. 4d) the diffusion of ash dominates over its advection and the air-laden ash spreads fan-like across Río Nergo and Neuquén provinces during 7 May, affecting lightly Viedma and Bahía Blanca.At about midnight of 8 May the cloud direction is already NNE but, again, it changes to ENE along the day.Moderate to low concentrations are predicted over a vast territory, including the Buenos Aires province.
The evolution of the ash cloud during the 2-8 May week is illustrated in Fig. 5, which shows the predicted cloud column mass at different days.Cloud column mass is computed by integrating concentration along the vertical from the surface to the top of the domain, and is a useful quantity (usually expressed in Tn/km 2 ) to compare with satellite imagery.Predictions for cloud trajectory and arrival times are in good agreement with observations (for comparison see the chronology of events in Sect.2). Figure 6 compares visible Moderate Resolution Imaging Spectroradiometer (MODIS) satellite images and predicted cloud column mass.It should be noted that these two types of images are not directly comparable because the MODIS ash detection threshold and the reflectivity coefficients of volcanic ash are not well constrained.However, the figure illustrates the capability of the model to predict the variation of the cloud position with time.
Figures 7 and 8 show, respectively, the predicted ground deposit load at different days and the expected deposit thickness at the end of the simulation (9 May at 00:00 UTC).The deposit thickness is estimated dividing the deposit load by the class-dependent ash density and assuming a package factor of 0.75.The model predicts accumulations of ∼10 cm at Futaleufú, 1-3 cm at Esquel and El Bolsón, few mm at Bariloche, and traces at most Atlantic coast sites.Such predicted thickness compare semi-quantitatively well with preliminary values reported by eyewitness.Figure 9      A deposit packing factor of 0.75 is assumed.
"proximal" and distal localities.Thus, for example, the simulation predicts a continuous fallout at Futaleufú during 2-5 May, with a higher accumulation rate during 4 May or, in the case of Comodoro Rivadavia, light fallout from 3 May at 07:00 UTC (03:00 LT) to 4 May at 15:00 UTC (11:00 LT).FALL3D calculates also concentrations at different vertical levels, a relevant quantity for aerial navigation.Figure 10 shows the evolution of concentration at flight level 300 (FL300), which corresponds to a standard pressure level of 30 000 feet (∼9 km) and is the normal flight level of aircrafts.Unfortunately, the threshold for hazardous concentration of ash is not well established.Notwithstanding this, we adopt the Montreal VAAC criteria and plot three different concentration contours of 10 −3 , 10 −4 , and 10 −5 gr m −3 (Witham et al., 2007).The 10 −4 gr m −3 contour is an estimate of the visual ash cloud.Assuming these values, the Esquel airport lays within the critical zone during the whole simulated period.Critical limits are also predicted at Comodoro Rivadavia during 4 May and at Trelew and Puerto Madryn airports during 5-6 May.The limit zone shifts towards NE during 7 May, affecting Bariloche and Esquel.Finally, during 8 May, the "conservative" isoline threshold contains also Santa Rosa and, to a lesser extent, Buenos Aires.This sequence of events agrees with what actually occurred during the 2-8 May week, as illustrated in Fig. 11.

Comparison between forecast and hindcast configuration results
Dispersion models need meteorological forecasts to run in operational mode.Global meteorological models like GFS give short and mid-term predictions (up to 10-15 days) and, consequently, ash-clouds and fallout from long-lasting eruptions could theoretically be forecasted several days in advance.However, such mid-term forecasting is likely to fail because meteorological forecasts loose accuracy in time (usually after the third day) and the evolution of the eruptive parameters is very often unpredictable.A more reasonable strategy for long-lasting events is to constrain the simulation to the next 48-72 h and then restart the run once the updated   forecasts and volcanological inputs are available.Here we have first simulated the first week of the Chaitén eruption using 3-days forecast intervals and then, at the end of the simulated period, we have compared the results with a single one-week run using the NCEP data analysis as described in Sect.3.1.Figure 12 compares the results from both simulations at 5 and 8 May (00:00 UTC), i.e. at the end of the 3-days interval, when maximum differences are expected.The figure illustrates how such differences are not very significant and evidences, at least for this particular event, the goodness of the 3-days forecast strategy.

Summary and conclusions
We simulated the transport and deposition of volcanic ash during the first week of May 2008 Chaitén volcano activity.The simulation spans from 2 May at 10:00 UTC to 9 May at 00:00 UTC.The purpose was to test the performance of the combined WRF+FALL3D models in forecasting volcanic clouds and ash fallout.The first model simulates meteorological conditions whereas the second deals with advection, diffusion and sedimentation of air-borne particles (ash in our case).As input parameters for the transport model we only used semi-quantitative information available shortly after the eruption onset.-In order to reproduce a column height of 14 km (average from observed values), the BPT predicts a MER of 3×10 7 kg s −1 .It would correspond to a VEI 4-5 eruption (Newhall and Self, 1982).
-Assuming a constant MER, the total mass erupted during the 2-8 May period is 1.7×10 13 kg, corresponding to ∼6 km 3 of volume (DRE).Note that it does not correspond to the in-land predicted deposit volume because part of the ash falls over the Atlantic ocean.We point out that this value is probably an overestimation due to the assumption of constant MER.
-The forecast predicts a cloud position in good agreement with subsequent observations: an initial rapid NE to SE cloud direction drifting during 3 May, a stabilized ESE direction along 4 May that becomes E on 5 May morning and, finally, a progressive NE drifting during the next days which are characterized by a diffusiondominant ash spreading due to the low wind intensities.
-The predicted deposit for the 2-8 May week is rather complex but shows two preferential directions, one along SE (along the Chaitén-Futaleufú-Comodoro Rivadavia alignment) and another along NE (along the Chaitén-El Bolsón-Neuquén alignment).Predicted thickness range from ∼5-15 cm at "proximal" localities (Futaleufú or Esquel) to few millimeters or less at more distant sites (like Neuquén or the Atlantic coast).
-The movement of the predicted FL300 critical concentration contour presents a notable coincidence with the actual status of the regional airports.
-The differences between simulations using 72 h meteorological forecasts, corresponding to a test of syneruptive model performance, and simulations using NCEP data analysis are low.
An important aspect during emergency situations is to dispose of aid-supporting predictions.A poor-accuracy but fast model can be, in this context, more useful for emergency managers than a high accuracy but time-consuming complex model.The total CPU time required for the WRF+FALL3D simulation was 3.5 h using an IBM Power PC 970MP processors cluster (128 CPU's for WRF and 42 for FALL3D).It corresponds to approximately 30 min of cluster CPU time per simulated day and for this large extension domain.We argue that such CPU times are acceptable even during emergencies and conclude that the WRF+FALL3D combination can predict efficiently, at both local and regional scales, cloud position, concentration on air, deposit load (thickness), and ash accumulation rates.
A major limitation of nearly real-time ash cloud forecasting comes from the large degree of uncertainty associated to the inputs of the models (for a detailed study on this subject see Scollo et al., 2008).Consequently, it is important to stress that forecasted quantities must be understood as a first order semi-quantitative approach.The results presented here can differ from reality for two main reasons.Firstly, we assumed a constant MER during the whole simulated interval, whereas the real Chaitén eruptive column showed multiple pulsating phases and two strong bursts on 6 May.Secondly, we considered a granulometry based on analogous sub-plinan eruptions.Clearly this is a rather subjective "input" but note that this is the "kind of data" that exists at the time of syn-eruptive forecasting, when no field measurements are available yet by obvious reasons.However, notwithstanding the remarks above, we point out the goodness of the combined WRF+FALL3D simulation when compared with MODIS satellite images and with semi-quantitative syn-eruptive reports on cloud arrival times and ash ground accumulations.
The Chaitén eruption was still going on at the time of writing this article.Certainly, more detailed studies including field data will be necessary to quantify this event properly.Nevertheless, preliminary results suggest that this eruption could be comparable to the 1991 Hudson event, which erupted about 7 km 3 of material (Scasso et al., 1994;Bonadonna and Houghton, 2005) and caused a severe economic recession in the Santa Cruz province.

Fig. 3 .
Fig. 3. Time evolution of the plume height according to the BPT for a MER of 3×10 7 kg s −1 .Variations in height are caused by different wind intensities.Note that the column heights do not correspond with the real values actually observed after the forecast.

Fig. 6 .
Fig. 6.Left: photo-like true color images from the MODIS on the Terra NASA satellite.The corresponding UTC time instants are, from top to bottom, 3 May at 14:35, 5 May at 14:25, 6 May at 19:15, and 8 May at 15:00.Right: predicted cloud column mass in Tn/km 2 at similar time instants.The inner frames indicate the extent of the MODIS counterpart image.

Fig. 11 .
Fig. 11.Time evolution of predicted concentrations at flight levels FL100 (blue) and FL300 (red) at three different localities: Esquel, Comodoro Rivadavia, Trelew, and Buenos Aires.The shaded zone indicates the actual periods in which the airports of these localities shut down or suffered flight grounding due to poor visibility.

Fig. 12 .
Fig. 12.Comparison of FALL3D results using meteorological forecasts (left) and hindcasts (right).Column mass is shown for 5 May and 8 May at 00:00 UTC, when maximum divergence is expected.The bottom plots show the deposit load predicted in both cases.