Models of volcanic ash (tephra) fallout are increasingly used to assess volcanic hazards in advance of eruptions and in near–real time. These models often approximate the wind field using simplistic assumptions of the atmosphere that do not account for four-dimensional variations in wind velocity. The fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5) is used to improve forecasts of tephra dispersion. MM5 is a 3D model that can predict circulations in areas with sparse meteorological observations and complex terrain, such as volcanic plateaus. MM5 is applied to the 1995 eruption of Cerro Negro, Nicaragua. Validation of MM5 is achieved by comparing the simulated winds with rawinsonde observations. Estimates of diffusivity and particle settling velocities are used in conjunction with MM5-generated wind fields to forecast the major axes of the tephra dispersion. The predicted axes of dispersion derived from the MM5 winds approximate very closely the observed bilobate tephra accumulation and tephra plumes observed in satellite images. MM5 winds provide far more accurate spatial and temporal forecasts than do the wind assumptions that had been used previously to assess Cerro Negro tephra hazard.
Tephra, commonly known as volcanic ash, is defined as airborne volcanic ejecta of any size (Tilling et al. 1987). It is a frequent product of volcanic eruptions. Volcanoes create numerous natural hazards that pose a risk to society, but wind-driven tephra threatens far larger areas than any other volcanic hazard. Tephra deposited on roofs can easily cause building collapse (Wolfe 1992). Tephra can severely damage airplanes, automobiles, and electrical equipment (Casadevall 1994; Scott et al. 1995; Newhall and Punongbayan 1996). Tephra often causes water contamination, and airborne tephra can cause adverse effects on respiration and irritation of eyes and nose (Latter 1989; Baxter 1990). Deadly mudslides known as lahars, made of tephra and water, often pick up more debris as they travel, and can move over 10 m s−1 (Iverson 1997). Tephra that reaches the stratosphere affects climate—for example, lower surface temperatures following the 1991 eruption of Mount Pinatubo (Robock 2002).
Because tephra poses numerous hazards, volcanologists have modeled its dispersion for just over two decades. Although tephra models have increased in complexity over time, their meteorological component has been conventionally simple, often out of necessity. As demand grows for real-time prediction and scenario-based hazard assessment of tephra deposition, improvements to the atmospheric component of tephra modeling are essential. There is an important need for volcanic eruption hazard mapping at many volcanoes around the world under different wind regimes. Many volcanoes are located in countries that do not have the resources required to undertake detailed hazard mapping. Atmospheric models in combination with tephra dispersion models can provide very valuable information to satisfy this requirement.
This study incorporates winds from a high-resolution atmospheric model to predict the major axes of tephra dispersion and deposition. Winds produced by the fifth-generation Pennsylvania State University–National Center for Atmospheric Research (NCAR) Mesoscale Model (MM5) are applied to the eruption of Cerro Negro (Fig. 1), Nicaragua, that occurred on 19 November–6 December 1995.
a. Tephra modeling
To understand tephra plumes better, many volcanologists have numerically simulated volcanic eruptions (Suzuki 1983; Carey and Sigurdsson 1982; Connor et al. 2001; Bonadonna and Phillips 2003). Physical assumptions are made about the eruption, eruption column, and tephra plume. Some eruption parameters are estimated. For example, satellite-estimated heights of eruption plumes can be prone to error (Tupper et al. 2004). Geologists often attempt to reconstruct eruption dynamics from ancient deposits. This activity is essential in volcanic hazard assessment. However, it is not always clear which aspects of the deposit reflect volcanological conditions, such as eruption column height, and which aspects of the deposit can be explained by meteorological conditions.
Modeling tephra movement in the atmosphere is crucial for modeling tephra accumulation. Once tephra particles in the eruption column have achieved their maximum height, their movement depends upon their characteristics (density, shape, terminal fall velocity) and, very important, wind. For most eruptions, tephra plumes are advected by winds and become more diffuse with time and distance downwind (Bursik 1998). In general, ballistic ejecta reach the ground very quickly, whereas the finest tephra settles last (Fig. 2). Vertical wind shear can create irregular deposition patterns. Although many tephra deposits exhibit axisymmetrical patterns, isopach maps of tephra deposits often reveal nonlinear axes, multiple deposition axes, and multiple thickness maxima (Ernst et al. 1994). Despite observations of complex tephra deposits, winds, more often than not, are assumed to be uniform (Hill et al. 1998). This assumption may be made for simplicity’s sake, or because wind data are not easily accessible and, in many cases, are incomplete. Wind data are extrapolated sometimes from a single observation from a weather station tens of kilometers away (Barberi et al. 1990; Macedonio et al. 1994; Hill et al. 1998). This is a serious problem because the simplistic wind assumption neglects a major impact on deposition location.
Much of the research on tephra dispersion involving atmospheric models has focused on aviation interests. The Volcanic Ash Forecast and Transport and Dispersion model (VAFTAD, Heffter and Stunder 1993); the Canadian Emergency Response Model (CANERM), a 3D Eulerian model (Simpson et al. 2002); and “PUFF” (Searcy et al. 1998) are designed for real-time tracking of volcanic ash clouds. These three models are intended for tracking middle-to-upper-tropospheric airborne ash. They do not predict low-tropospheric dispersal.
Atmospheric dispersion models applied to surface tephra accumulation have not been as widely investigated as those used for aviation. Suzuki (1983) was the first to model tephra transport and deposition numerically. He used a diffusion–advection equation for horizontal particle motion in a temporally and spatially uniform wind field. The Suzuki model has been adjusted to represent increasingly intricate volcanic eruptions (Glaze and Self 1991; Armienti et al. 1988; Connor et al. 2001; Bonadonna et al. 2002; Bonadonna and Phillips 2003) and modified for parallel computation in the “TEPHRA” model (Bonadonna et al. 2005). Although TEPHRA is sophisticated in its representation of volcanic eruptions, it uses a single wind velocity profile as its input and assumes uniformed winds over each x–y domain for each altitude.
In a study similar to ours, Turner and Hurst (2001) used atmospheric and dispersion models to forecast tephra deposition from Mount Ruapehu, New Zealand. They compared “ASHFALL” (Hurst and Turner 1999), a simplified Suzuki model, with a combination of the Regional Atmospheric Modeling System (RAMS) and the Hybrid Particle and Concentration Transport model (HYPACT). Although HYPACT had merit in tracking tephra plumes, it was not suitable for quantifying tephra deposition depth because of shortcomings related to the distribution of particle fall velocities. However, the RAMS–HYPACT suite produced more accurate forecasts than did ASHFALL and demonstrated some benefits of incorporating high-resolution, spatially varying winds into a volcanic ash dispersion forecast system. The Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Draxler and Hess 1997) is another program used for modeling the dispersion of volcanic ash in the atmosphere. The Active-Tracer High-Resolution Atmospheric Model (ATHAM), a nonhydrostatic plume model, simulates the dispersion of gas and particles (Herzog et al. 2003).
b. Mesoscale atmospheric modeling and dispersion modeling
Mesoscale models are of particular relevance to tephra modeling because they address the same spatial scales (on the order of 10–100 km) as a typical tephra plume (Brasseur et al. 1999). MM5, like other mesoscale atmospheric models, simulates an atmosphere that evolves based on the physical laws of motion, and the laws of conservation of energy, mass, and momentum. It is recognized that high-resolution grids (kilometer scale or smaller) are necessary to depict wind patterns in regions of complex terrain (Cox et al. 1998; Salvador et al. 1999; Grell et al. 2000). Various implementations of MM5 have competently simulated orographically driven circulation such as gap winds, downslope winds, lee waves, and topographically forced precipitation (White et al. 1999; Zängl 2002; Rife et al. 2004). Those models have been effective at resolving circulations in areas with sparse data and/or complex terrain—hence, the suitability for simulations over Nicaragua. For example, MM5 was applied to tropical, mountainous Colombia where it replicated diurnal shifts in precipitation, as well as other localized phenomena (Warner et al. 2003). During the Mesoscale Alpine Program, a version of MM5 performed well in complex terrain (Ferretti et al. 2003). A 1998 sand “event” in China was simulated using winds generated using MM5 (In and Park 2002). The winds were fed into an aerosol transport model that predicted well the arrival of suspended sand in Korea. MM5 winds were used in a transport–chemistry–deposition model (Kitada and Regmi 2003) for the mountainous region of Kathmandu, Nepal. It revealed large diurnal changes in wind patterns, as well as seasonal peaks in air pollution.
Because versions of MM5 have proven adept when applied to various dispersion problems, it should prove useful for modeling tephra transport and deposition. Therefore, the aim is 1) to evaluate the accuracy of MM5 forecast winds for the study domain, 2) to predict the major axes of tephra dispersion, and 3) to verify forecasts of tephra deposition using observed tephra deposit and satellite images. In a hazard response environment when the balance between accuracy and timeliness of information is critical, questions can be raised about how complex of a wind field is necessary to predict tephra dispersal. Is the simple assumption of uniform winds sufficient for assessing the tephra hazard for similar eruptions? Do readily available global analyses produce an accurate forecast? Global reanalyses from the National Centers for Environmental Prediction (NCEP) and NCAR were used for probabilistic assessment of tephra dispersal by Bonadonna et al. (2005) and wind fields used in PUFF are derived from NCEP twice-daily forecasts (Searcy et al. 1998). To assess the adequacy of forecast wind fields for this type of eruption, MM5 simulations are performed with two minimum grid sizes. The MM5-driven trajectories are compared with those derived from the uniform wind and the NCEP–NCAR reanalysis.
2. Volcanological setting
Cerro Negro volcano is located in the central Marrabios Range, which runs parallel to Nicaragua’s Pacific Ocean coast. The chain runs northwest to southeast throughout Central America, marking a convergent plate boundary where the Cocos plate is subducted under the Caribbean plate. Cerro Negro is 728 m in elevation and is located at 12.5°N, 86.7°W. It sits 20 km east-northeast of León, Nicaragua’s educational and cultural capital and its second-largest city (estimated 200 000 residents).
a. 1995 eruption parameters and tephra deposit pattern
A basaltic cinder cone, Cerro Negro has had 23 eruptions since its formation in April of 1850. According to the Global Volcanism Network (1995), the 1995 eruption occurred between 19 November and 6 December. On 25 November, a possible low-level ash cloud moving toward the west-southwest was seen on visible satellite images. Between 27 and 28 November, the ash plume shifted toward the northwest. The period of the most prolific tephra output was toward the west and west-southwest from 29 November to 2 December. The total volume of tephra erupted was estimated at a dense rock equivalent of 1.3 × 106 m3 (Hill et al. 1998).
The eruption produced a bent-over plume, indicating that the vertical velocity of the eruption column was the same order of magnitude as the crosswinds it encountered. The maximum column height was estimated at 2400 m (Hill et al. 1998). The volcanic explosivity index (VEI; Newhall and Self 1982) for the 1995 eruption was 2. For comparison, the 1991 eruption of Mount Pinatubo had a VEI of 6 and the 1992 eruption of Cerro Negro had a VEI of 3 (Connor et al. 2001). Although the 1995 eruption did not cause any human fatalities, it forced the evacuation of 12 000 people, of whom 6000 were from 15 rural villages, and caused U.S. $700,000 of damage to agriculture and infrastructure (Connor et al. 2001).
The tephra deposit from the 1995 eruption was sampled immediately after the eruption (Hill et al. 1998). Thicknesses of 30 cm extend 2 km west of the vent, and the deposit was about 0.5 cm thick as far away as León (Fig. 1b). The most striking feature of the deposit is that it has two lobes—one running west-southwest to León and the other west-northwest toward the town of Telica.
3. Model descriptions
MM5 is a nonhydrostatic limited-area model that uses a terrain-following vertical coordinate system, a one- or two-way nested grid system, and full physics (Grell et al.1995; Dudhia et al. 2002). Because this study modeled winds that are in the wake of Cerro Negro, a 1-km grid is utilized. To determine whether a lower-resolution grid would be sufficient for predicting the tephra dispersal, MM5 is also run with 9-km minimum spatial grid. Simulations are performed with three nested domains and five nested domains (Fig. 3, Table 1). All grids use 23 terrain-following or sigma vertical levels extending from approximately 15 m to about 17 km AGL. The smallest vertical interval (500 m) is in the boundary layer. MM5 was initialized with data from the NCEP–NCAR reanalysis project (Kalnay et al. 1996). These retrospective global data are on 2.5° grids at 17 different pressure levels, from 1000 to 10 hPa. This simulation uses the Blackadar planetary boundary layer scheme for all domains, and the Grell cumulus scheme is used for the three outer domains. Four-dimensional data assimilation (FDDA) is used to “nudge” forecast values toward global analyses, every 12 h, to reduce the likelihood of output that substantially diverges from realistic values over time. FDDA was recommended for longer simulations by Grell et al. (2000).
Computing costs are minimized by utilizing the MM5 multiprocessor mode on a Linux cluster. An eight-processor mode was found to be the most stable configuration on the cluster and each minute of computing time accounted for 4 min and 40 s of forecast time.
To account for spinup, as the mesoscale model is given time to create dynamically and physically consistent meteorological fields, the model was started on 25 November, although the continuous phase of the eruption did not begin until 27 November. This configuration of the model enables assessment of the effects of small-scale topographic flows on ash dispersal and deposition and provides guidance to emergency response for mountainous, rural villages.
b. Particle fall model
In tracing the movement of tephra particles, some simplifying assumptions are made. First, tephra particles from the 1995 eruption are all assumed to be released into the atmosphere at the exact same location, an initial height of 2400 m above sea level, directly above the center of the volcanic vent. This model is a dramatic simplification but was deliberately chosen because it allows us to isolate the effect of the wind field on tephra fallout from the effects of other processes, such as convection in the ascending volcanic plume, that also influence tephra transport.
Although friction in the vertical direction is used to determine the terminal fall velocity of each particle, horizontal motion of tephra particles is determined solely from MM5-predicted winds. Horizontal friction and turbulent eddies are neglected. These assumptions are suitable for simple eruptions in which plumes are not very diffused. The impact of different wind fields can then be assessed from the major axes of dispersion rather than dispersed particle fields.
The 3D gridded wind fields enclose the 1995 tephra deposit out to the 0.1-cm isopach (Fig. 1b). The u and υ wind fields are stored within a 40 × 32 (horizontal dimensions) × 34 (vertical dimension) grid. The respective position of each released particle is advanced based upon the velocity field in the given grid cell. The particle falls vertically according to its terminal fall velocity until its z coordinate reaches zero, which is sea level.
As a particle is “released” into the grid, it is advected in 1-s increments by the wind vectors it encounters. A positive value for u means winds are blowing from west to east, and a positive value for υ means winds are blowing from south to north. Tephra particles are generally oblong and cannot be regarded as spheres in modeling (Armienti et al. 1988). The settling velocity of a particle is tied mainly to its dimensions and density. Particles between 1 μm and 1 cm can reach their settling velocity in 5 m or less, and so they can be assumed to have this fall velocity upon release (Armienti et al. 1988). Particles fall vertically as determined by Suzuki’s (1983) equation for settling velocity:
ρt is tephra density (kg m−3), ρa is air density (kg m−3), η is air viscosity (N s m−2), g is gravitational acceleration (m s−2), pf is particle shape factor, pa, pb, and pc are relative dimensions of the three axes of a particle, with a being the major axis, and ϕ is particle diameter.
Most of the variables in Eq. (2) are held constant at assumed values. The tephra density ρt is set at 1100 kg m−3. Air density ρa is assumed to be 1.229 kg m−3, whereas the true air density is inversely related to altitude. Air viscosity η is set at 1.73 × 10−5 N s m−2, and g is 9.8 m s−2, a common assumption. The particle shape factor pf is a dimensionless parameter held at 0.5. This value would occur for tephra particles for which the major axis is 2 times as long as either of the minor axes of the particle. Particle size is the one “variable” in the Suzuki equation that was varied. Diameter ϕ was tested at 0.26, 0.5, 1, 2, and 5 mm. The 0.26-mm particle size was chosen because it produced a simple terminal fall velocity of 1.0 m s−1, which made error checking easier. Because air viscosity and air density are held constant, this particle fall model should only be applied to the lower atmosphere.
4. Observations and model data
The data and method are illustrated in Fig. 4. The NCEP–NCAR data were obtained from NCAR Scientific Computing Division (SCD) archives in gridded binary (GRIB) format, which enabled direct input into the MM5 modeling system. Surface observations and soundings, also from the NCAR SCD archive, were objectively analyzed in the mesoscale domain.
Sounding data from 19 observation stations (Fig. 3) were employed for model validation. The station sites are variable in terms of the local terrain and vegetation characteristics. The number of observations available at each station varied. Nine stations had all observations (0000 and 1200 UTC) for the 9-day period.
Wind components (u and υ) and geopotential height were extracted from the MM5 output for a 40 km (x direction) × 32 km (y direction) × 20 vertical pressure level grid centered on Cerro Negro, which encompassed the grid for the particle fall model. The pressure levels ranged from 1000 to 100 hPa, at 50-hPa intervals. Then, the u and υ winds were vertically interpolated to 500-m incremental levels from 0 to 16.5 km. The horizontal grid spacing was approximately 1 km, and its extent was based on the geometry of the 1995 tephra deposit.
Multispectral satellite images from the Geostationary Operational Environmental Satellite (GOES) were used to track the volcanic ash plume and to verify the simulations. Visible images were used to track the plumes during the day and split-window differencing of three infrared bands adapted from Ellrod et al. (2003) were used to track the plume at nighttime. These methods, while not uniformly effective at classifying ash clouds (Oppenheimer 1998; Tupper et al. 2004), were sufficient because of the absence of major weather systems during the simulation. Only on 27 November was there strong convection associated with an upper-level trough and frontal passage.
Predicted tephra deposit axes using MM5 winds were compared with those produced using NCEP–NCAR reanalyses only and uniform winds across the domain. Using the geopotential heights, the reanalyses winds were first interpolated vertically to the grid of the particle fall model. The wind data at each altitude were then interpolated horizontally, from the two closest global analysis grid points to each of the 1280 x–y grid points in the particle fall grid, using an inverse-distance-weighted mean.
5. Mesoscale model validation
MM5 forecasts were compared with atmospheric soundings to ascertain the applicability of this model for this region and period. Model-produced soundings were generated by horizontally interpolating MM5 output to the latitude–longitude coordinates of each of the meteorological stations and vertically interpolating from terrain-following sigma levels to pressure levels. The variables tested statistically were twice-daily winds (u and υ components) and geopotential heights h for the 850-, 500-, and 200-hPa levels. These heights are representative of the atmospheric boundary layer, the middle-tropospheric steering level, and the upper troposphere. To quantify model error over time and space, root-mean-square (RMS) errors were computed at each of the 19 meteorological stations, at three levels, every 12 h from 0000 UTC 27 November to 0000 UTC 4 December.
In general, the simulations are not expected to be perfect, but rather to be reasonable. The observational data being used to initialize the model are scarce; all soundings are more than 380 km from the volcano and none exist for Nicaragua. It is impossible to directly validate MM5 in the Cerro Negro vicinity.
Two separate sets of RMS errors were calculated for model validation: station averages for the entire simulation period and hour-by-hour time series. The difference between the model-predicted value and the observation given by a sounding was defined as the model error for a particular variable at each time, pressure level, and station. In the station-by-station method, the squared errors of each variable (u winds at 850 hPa, height at 200 hPa, etc.) were summed over the 9-day model run and divided by the number of entries, and the square root of this quotient is the RMS error for the variable. The RMS errors for each station are summarized in Fig. 5, which shows that there are no spatial trends in the average errors across the domain. At 850 hPa (Fig. 5b) the component directional errors were highest over Costa Rica and east of Nicaragua. For most of the domain, component errors were 3 m s−1 or less. The model consistency was evaluated by comparing the errors hour by hour. This required summing the squares of the errors of all stations then taking an average, for each sounding time.
For zonal (u) winds, the RMS error at 850 hPa was the lowest for 13 of the 19 observation times (Fig. 6a). The model performed better for observations that were recorded at 0000 UTC (hours 0, 24, . . .), as opposed to 1200 (hours 12, 36, . . .). Differences may be a result of MM5 not adequately representing changes in the boundary layer that are due to the diurnal cycle.
For the meridional (υ) winds, MM5 performed best at the 500-hPa level, with an RMS error mostly between 1 and 3 m s−1 (Fig. 6b). The error for 850 hPa was only slightly higher, and the 200-hPa level had more variable RMS errors, mostly ranging from 3 to 6 m s−1. The errors were relatively consistent over time.
a. MM5 validation summary
Model validation has demonstrated that MM5 produced a reasonable simulation of wind fields for the period of the 1995 Cerro Negro eruption. Model error did not grow significantly with time. This result is important, considering the error trends exhibited in the study by White et al. (1999), in which RMS error increased in almost every subsequent model time step. Average geopotential-height RMS errors for this study were considerably lower on the whole than those reported by White et al. (1999). Wind-component error values are consistent with other simulations in mountainous terrain, such as Ferretti et al. (2003). Averages across the horizontal domain did not reveal any spatial trends. Further details of the model validation are available in Byrne (2005).
6. Model prediction—Tephra
a. Prediction of major axis of dispersion
The axis trajectories predicted by the high-resolution MM5 wind field correspond very well to deposit measurements. Figure 7 shows trajectories generated for the five chosen particle sizes at 0000 and 1200 UTC.
For the 1-km MM5 winds, only 2 of the 60 particle paths terminate outside the 0.1-cm isopach, and the path locations appear to mimic the thickness contours well on the whole. These trajectories can account for deposition along both major axes of the deposit. The trajectories demonstrate that wind speed and direction changed significantly during the modeled period. The more distant paths that reach León represent stronger winds (nearing 9 m s−1) encountered by the smallest particles (0.26 mm in diameter). The major axis of dispersion shifts between the azimuths of 240° and 280° over the course of the 9-day model run, a range of 40°. For the most part, the winds encountered by the majority of erupted particles are just north of easterly (70° azimuth, toward 250°), as evidenced by the deposit.
It appears that the 9-km MM5 winds (Fig. 7b) shifted the axes in a more southwesterly direction relative to the 1-km output (Fig. 7a). This result could be due to the decreased influence of localized topographic effects.
The particle trajectories produced by the reanalysis winds (Fig. 7c) seem to follow two major axes that are 15° out of phase with the deposit axes. The reanalysis winds have a stronger northerly component, which results in tracks that push too far to the south. This result is evident in particular at 0000 UTC 27 November, 29 November, and 1 December. Those winds also produce fewer northward-moving particles. In addition, the trajectory paths on 27, 28, and 30 November and on 1 December are shorter than the MM5-driven trajectories.
The uniform wind (Fig. 7d) blowing from 70° at 9 m s−1 (Hill et al. 1998) matched the large southern lobe of the deposit very well but completely failed to account for the deposition in the northern lobe. Of the three wind fields, MM5 winds had the highest spatial variability closer to the ground and were more variable in the first few days of the simulation, 27–29 November.
A vertically exaggerated 3D view of particle trajectories, driven by the 1-km MM5 winds, exhibits how the northwestern particle tracks can account for the creation of the smaller northwestern lobe of the 1995 deposit (Fig. 8). These northwestern tracks were generated from wind fields on 27, 28, and 30 November, all at 1200 UTC. According to the MM5 simulations, there is a diurnal wind shift during the first few days of the constant-tephra-production phase. Morning winds (1200 UTC) have a greater tendency to have a southerly component, blowing tephra to the northwest, whereas evening winds (0000 UTC) tend to have a northerly component, distributing tephra toward the southwest.
b. Wind field comparison—Satellite imagery
Satellite images have proven to be one of the most useful tools for detecting eruptions and tracking volcanic ash plumes (Oppenheimer 1998; Tupper et al. 2004). Although in meteorology winds are conventionally described by the direction from which they are blowing, it is more intuitive to describe the tephra plume in terms of the direction in which it is moving. Because the plume in these satellite images (e.g., Fig. 9) consists of primarily finer-grained tephra than what was selected in the particle fall model, a smaller particle size of 0.01-mm diameter was tested to simulate particles that stay in the atmosphere longer. These particles remained in the 2000 ± 250-m wind level for the duration of their paths in the wind grid, over 40 km from the vent. The 0.01-mm-diameter particles were tested with the MM5 winds, NCEP–NCAR reanalysis winds, and uniform winds. The most important shifts in the plume direction and corresponding predicted trajectories are presented in Table 2.
At 1315 UTC 28 November, the plume is blowing to the north-northwest at an azimuth of about 320° (Fig. 9a). MM5 (1 km) simulates a track at 315° at 1200 UTC 28 November. The 9-km MM5 winds place the plume track slightly more to the west at 310°, and the NCEP–NCAR winds predict the plume axis at 285°.
The plume shifted to the west-southwest on 29 November (Fig. 9b), is visible over the Pacific Ocean, and appears to be spreading laterally at a constant rate. The MM5 1-km trajectory is toward 280°, and the 9-km-grid forecast is toward 270°. In this case, the NCEP–NCAR forecast of 240° is closest to the satellite observation.
At 2215 UTC 30 November (Fig. 9c), the plume is blowing toward approximately 250° and is clearly visible to the Pacific Ocean, at which point it appears to bend toward the west-northwest. At 0000 UTC 1 December, MM5 simulates the plume to move toward 250° and 245° for the 1- and 9-km grids, respectively. NCEP–NCAR reanalysis winds bring the plume farther south at 235°. On 1745 UTC 1 December (Fig. 9d), the plume fans between 265° and 250° and is visible over the ocean. At 0000 UTC 2 December, MM5 1-km winds move toward the west and then bend toward the west-southwest at 250°. The 9-km gridded MM5 winds take the plume toward 260°, and the NCEP–NCAR winds move the plume toward the 240° azimuth.
In summary, the 1-km MM5 winds exhibited the highest spatial variation and fit the tephra deposition pattern better than did the three other wind models. MM5 wind-predicted trajectories from both resolutions were found to be accurate within 10° in azimuth relative to satellite-observed plume direction, with the exception of the evening of 29 November. NCEP–NCAR reanalysis winds agreed with the deposit somewhat less, forecasting particle tracks too far to the south. Errors between NCEP–NCAR reanalysis winds and satellite observations were significantly higher than those of MM5, on the whole. The uniform winds agreed only with the southwestern lobe of the tephra deposit.
This study has shown that a mesoscale atmospheric model such as MM5 can provide a relatively accurate and high-resolution depiction of the winds that a volcanic plume will encounter after an eruption. This result supports the findings of Turner and Hurst (2001), a study that used the RAMS over New Zealand. MM5 has great utility in regions where wind data are not readily available and where mountainous terrain causes localized wind effects for which a nearby sounding could not substitute. Moreover, when strong eruptions are able to reach the upper troposphere and stratosphere, tephra can be advected over large areas where winds vary in magnitude and direction, rendering the assumption of uniform wind woefully inadequate.
MM5 could be extremely useful for short-term forecasting of tephra dispersion for real-time eruptive events. Given adequate input data for initial and boundary conditions, MM5 is a dependable forecasting tool. Upper-air sounding data taken near an erupting volcano are certainly useful, but the data cannot be extrapolated very far into the future. A forecast from MM5, on the other hand, could provide a more comprehensive prediction of winds.
The application of atmospheric dispersion models to ground accumulation of tephra has not been widely used operationally for emergency-response procedures. The results presented in Fig. 7 are important in illustrating the impact of grid resolution on model-based dispersion of tephra, a result that has implications in emergency-response situations in which the value of increased accuracy versus expedited information must be assessed. Turner and Hurst (2001) compared 2.5- and 40-km-grid simulations and found an approximate difference of 10° in the horizontal orientation between the two sets of simulated plumes. They were unable to verify the quality of the high-resolution RAMS simulations but suggested that for operational purposes the 40-km-grid winds would be adequate. They found that about 35° in azimuth was a useful range for ash-fall prediction if accurate meteorological input is assumed while allowing for uncertainties in the eruption timing and strength. Both sets of our MM5-driven trajectories were accurate to within 10° azimuth of satellite plumes. Notwithstanding the adequacy of the low-resolution guidance, we agree with the recommendation of Turner and Hurst (2001) that, if computationally practical, the increased accuracy to be gained from the higher-resolution, varying winds makes their use imperative in any tephra-dispersal forecast system.
Ernst et al. (1994) stated that diffusion–advection models that predict tephra deposition are applicable for eruptions with strong plumes but not for weaker plumes that can bifurcate. A plume with two vortices and constant buoyancy should split into two distinct plumes that continue to separate as they move downwind. While tracking the central axis of dispersion would be meaningless for a plume with two distinct axes, satellite imagery from the 1995 Cerro Negro eruption essentially dismisses this possibility. Imagery indicates that the weak 1995 Cerro Negro plume did diffuse significantly but did not bifurcate. Therefore, shifts in wind direction primarily account for the deposition pattern.
a. Chronology of deposition
The MM5 winds suggest that the smaller, northwestern lobe of the 1995 Cerro Negro tephra deposit was created primarily before and at the beginning of the major phase of the eruption, on 27 and 28 November (Fig. 7). Satellite imagery (Fig. 9) indicates that winds had a southerly component during this time.
The larger southwestern lobe of the deposit was then created after 28 November. In several satellite images, tephra is seen reaching all the way to the Pacific Ocean (e.g., Figs. 9b,c). These particles are finer than most of the particles modeled in this project and did not result in a measurable deposit.
b. Model limitations
Although the nested-grid setup maximizes the efficiency of MM5, spatial-resolution and temporal-resolution problems still persist. The innermost model domain (1 km) is not small enough to accurately describe Cerro Negro itself (1-km diameter).
A major drawback in the particle fall model and, indeed, in many volcanic plume models is that they neglect the effects of particle aggregation. Although there was little evidence of particle aggregation during the 1995 eruption of Cerro Negro, this phenomenon is believed to take place in many tephra-producing eruptions (Brasseur et al. 1999). With particle aggregation, clumps of fine tephra are removed from the atmosphere more quickly than they otherwise would be, because clusters have a larger diameter and thus a higher fall velocity. Veitch and Woods (2001) detailed examples of eruptions in which particle aggregation significantly changes the pattern of tephra deposit.
The volcano is an additional source of heat, moisture, and momentum that has not been considered in the model specification. If nearby surface and sounding data were available during the eruption, those data could have been assimilated into the model and would have influenced the forecast.
Although MM5 produces higher-resolution winds than most wind fields used in volcanology, it cannot recreate atmospheric phenomena smaller than a few kilometers in width (Dudhia et al. 2002). For this reason, small-scale eddies and turbulence within the plume cannot be simulated explicitly. Because the particle fall model only tracks particles along the major axis of dispersion, however, the diffusive effects of turbulence are not of primary concern.
The fall velocity of each particle is calculated solely as a function of its size, because all particles are assumed to have the same shape and density. Changes in air density and viscosity are not considered. For this particular eruption, variability of the air column is not as crucial, because the particles are only falling through the bottom 2400 m of the atmosphere. For more powerful eruptions, especially those whose columns reach the tropopause, changes in the air density become much more important to consider.
To simplify the particle fall program, the tephra particles in the model are all assumed to fall to sea level. True elevations in the region of the deposit range from 50 m toward the Pacific up to the base of the volcano at 450 m.
c. Sources of error and recommended improvement
Several sources of error can complicate a model run—for example, grid resolution, physics schemes, initial conditions, and boundary conditions. Despite these possible sources of error, model validation showed that this configuration of MM5 produced sufficiently reasonable forecasts. Note that this configuration of the model was constrained by the assimilation of the global reanalyses. Nevertheless, the predicted trajectories demonstrated greater accuracy at the highest resolution, where local and mesoscale effects such as topography have more influence. With the lower-resolution grid, the influence of the global analyses is evident in the shift toward more northeasterly winds and southerly trajectories. In an operational scenario, a mesoscale model like MM5 would be nested within a global forecast model, and the uncertainties in the forecast tephra dispersion will increase accordingly.
Although MM5 is a very complex atmospheric model, the particle fall program used to track tephra for this study was simplistic. A number of logical updates to the model would make it more robust and realistic. One relatively simple addition would be to incorporate varying release heights for the tephra particles, based on ground observations or satellite imagery. Another modification to the particle fall model would allow it to take into account the effects of wind during buoyant column rise. This effect is more applicable to bent-over plumes, for which the vertical velocity of the eruption column is close to the magnitude of the crosswind. Another progression of the model would be to add a diffusion term so that spatial distributions of tephra on the surface can be forecast. Unlike the MM5 model, which uses topography in creating a wind field for a particular region, the particle fall model does not have a built-in digital elevation model. Another advance, then, would be to assimilate topography into the fall model so that a particle would cease moving when it reached the true surface elevation at a given x, y location.
8. Concluding remarks
A mesoscale atmospheric model, MM5, was applied to the 1995 eruption of Cerro Negro and reproduced a reasonably realistic four-dimensional wind field, confirmed by tephra deposition patterns, satellite imagery, and meteorological station soundings. The major axis of dispersion for this eruption was mapped as a series of particle trajectories for the period of highest tephra output and showed very good agreement with the observed tephra deposit. It was also found that the axis of dispersion shifted significantly in space and time over the course of the eruption. Ninety-seven percent of the axis tracks fall within the observed pattern of tephra deposition for the 1995 eruption. This result is very promising, considering that MM5 had no input station observations within 380 km of the volcano and did not utilize any data extrapolated from the tephra deposit.
Although the eruption parameters and atmospheric conditions seemed to be favorable for plume bifurcation (Ernst et al. 1994), MM5 winds and satellite imagery show that the bilobate tephra deposit from the 1995 eruption can be accounted for by variation of winds. In addition, evidence shows that the northern lobe of the tephra deposit was produced, in large part, at the outset of the continuous phase of the eruption on 27–28 November.
Particle trajectories revealed spatial and temporal differences among the dispersion patterns produced by MM5 winds, NCEP–NCAR reanalysis winds, and uniform winds. Reanalysis winds, scaled to the small domain, forecast the plume to be too far to the south most of the time. The uniform winds were adequate in producing the larger, southern lobe of the deposit but could not account for the northern lobe. The uniform winds match satellite observations about one-half of the time. The significant shifts in wind direction, verified by satellite imagery, show that the assumption of uniform winds, even for a relatively weak eruption such as this, is faulty. Model winds vary over x–y space noticeably but vary by much greater differences for small changes in altitude. The differences among the forecast patterns illustrate the impact of grid resolution on the model-based dispersion of tephra. When timeliness rather than greater accuracy is critical, simpler model configurations would be optimal. The 9-km MM5 wind field would be adequate for this kind of eruption because it produced trajectories that matched the tephra deposit pattern fairly closely.
In addition to providing improved guidance for surface accumulation hazards, this paper presents a low-level case that is important to major airlines for approach and departures from an airport in the path of the ash but is most important to regional carriers who fly low-level routes. This investigation also provides another link between volcanology and meteorology, using methods that already have proven useful in their respective fields. More such collaboration is needed to understand tephra dispersal better and to improve volcanic hazard mitigation.
Thanks are given to Cindy Bruyere, Dr. Wei Wang, Brian Smith, Dustin Binge, and Laura Connor for assistance with running the models. Suggestions by Drs. Jason Knievel and Thomas Warner and by three anonymous reviewers greatly improved the manuscript. USF Research Computing Core, NCAR SCD, and NCAR/MMM provided computing facilities. Satellite images were provided by the Space Science and Engineering Center of the University of Wisconsin, sounding data were provided by the University of Wyoming Web site, and NCEP–NCAR reanalysis data were provided by the NOAA–CIRES Climate Diagnostics Center, Boulder, Colorado, from their Web site (http://www.cdc.noaa.gov/). Part of this work was performed while the first author was a visitor at NCAR/MMM. This research was supported by NSF Grant EAR-0130602.
* The National Center for Atmospheric Research is sponsored by the National Science Foundation
Corresponding author address: Arlene Laing, NCAR/MMM, P.O. Box 3000, Boulder, CO 80307. Email: email@example.com