Forecasting volcanic ash deposition using HYSPLIT

Abstract

A major source of error in forecasting where airborne volcanic ash will travel and land is the wind pattern above and around the volcano. GNS Science, in conjunction with MetService, is seeking to move its routine ash forecasts from using the ASHFALL program, which cannot allow for horizontal variations in the wind pattern, to HYSPLIT, which uses a full 4-D atmospheric model. This has required some extensions to the standard version of the HYSPLIT program, both to get appropriate source terms and to handle the fall velocities of ash particles larger than 100 microns.

Application of the modified HYSPLIT to ash from the Te Maari eruption of 6 August 2012 from Tongariro volcano gives results similar to the observed ash distribution. However, it was also apparent that the high precision of these results could be misleading in actual forecasting situations, and there needs to be ways in which the likely errors in atmospheric model winds can be incorporated into ash models, to show all the areas in which there is a significant likelihood of deposited ash with each particular volcanic eruption model.

Introduction

One aspect of GeoNet, the geological hazard monitoring system for New Zealand run by GNS Science, is forecasting the expected thickness distribution of volcanic ash if any volcano was to erupt. As well as the four volcanoes that have erupted in the last 50 years (Ruapehu, Ngauruhoe, Tongariro and White Island), possible ash distributions are calculated daily for six other volcanoes seen as significant risks, as shown in Table 1. There are six cases calculated for each volcano, with two column heights and three eruption volumes. If any volcano is erupting or has erupted recently, the forecasts for that volcano are distributed to relevant organisations and the public, whereas during normal times the information is kept within GNS Science to avoid causing confusion.

Table 1 The standard eruption model input parameters used in New Zealand by GNS Science, with three volumes and two column heights for each volcano

These forecasts are currently produced by the program ASHFALL, originally described in Hurst (1994), which is now somewhat out-dated. ASHFALL follows the ash particles from an initial vertical eruption column as they are moved horizontally by the wind, and at the same time fall vertically with a velocity which depends on their size and density. It was designed to have a minimal processing requirement, based on the processing power of mid-1990s Personal Computers. To obtain reduced processing times, ASHFALL ignored vertical diffusion of the ash, although its effect was partly accounted for by using a large value for the horizontal diffusion coefficient. This horizontal diffusion coefficient is effectively a “fudge factor”, which spreads out the expected ash distribution to incorporate uncertainties in the wind pattern and other parameters. ASHFALL uses MetService (New Zealand Meteorological Service) forecasts of the wind at various levels above each volcano. However, it does not allow for any horizontal variations in the wind pattern, i.e. ASHFALL can only use 1-D spatial wind fields, although it does allow for wind changes with time.

Hurst and Turner (1999) found that the forecast wind field was the major source of error in the volcanic ash distribution from Ruapehu eruptions. During investigations as to how the wind information used by ASHFALL might be improved, we looked at the atmospheric dispersion model HYSPLIT (Stein et al., 2015), already in use by MetService for the movement of fine particles travelling in the atmosphere, such as volcanic ash, and found it has the capability for estimating the distribution of ash particles on the ground. As well as capturing the influence of realistic atmospheric motion in ash transport, HYSPLIT can also model in-cloud (rainout) and below-cloud (washout) wet deposition processes.

HYSPLIT is a hybrid Lagrangian dispersion model, developed by NOAA/ARL which is used by MetService, in its role as a Volcanic Ash Advisory Centre (VAAC), to model airborne volcanic ash, with meteorological data provided by external and in-house NWP (Numerical Weather Prediction) models, which operate with 3 spatial dimensions and time. The spatial resolution of these models is 4 km, and the temporal resolution one hour. HYSPLIT is operationally used by several of the 9 Volcanic Ash Advisory Centres for aviation forecasting. A by-product of the HYSPLIT volcanic ash dispersion simulations is the ash deposition at the ground surface. The ground surface elevation is modelled at the same resolution as the atmospheric model. Aside from the dispersion of volcanic ash, HYSPLIT is used in several other atmospheric transport applications, including the dispersion of hazardous materials (e.g. nuclear material after the Fukushima reactor accident), air quality modelling (e.g. ozone, visibility, haze, and dioxin), dust storms, smoke, and the transport of biological material (e.g. pollen and mould spores) (Stein et al., 2015). Rather than using a horizontal dispersion coefficient, the dispersion was calculated from the friction velocity, height and boundary layer height following Kantha and Clayson (2000), where there are different equations for stable/neutral and unstable layers for both the surface and boundary layer. Above the boundary layer a mixing coefficient obtained from mixing length theory is used to calculate the velocity variances (Draxler and Hess, 1997).

We initially compared the ash deposition results from HYSPLIT with ASHFALL for similar eruptions and wind patterns. This showed that alterations to the standard fall velocity model of HYSPLIT were required to deal with ash particles larger than about 100 microns, which make up the bulk of ash deposits near a volcano. This is not a consideration in the aviation applications of HYSPLIT, as these particles are the ones that rapidly fall out of ash clouds.

Ash sizes and fall velocities

HYSPLIT was originally designed to use atmospheric models to calculate the tracks of pollutants in the atmosphere (Stein et al., 2015). In regard to volcanic ash, it is used for tracking the small ash particles with long residence times that are hazardous to aircraft generally at considerable distances from the volcano. Accordingly by default it only considers ash particles less than 100 microns in diameter (e.g. in VAFTAD, Heffter and Stunder (1993)). The terminal velocity for these small particles can be calculated by Stokes Law, which applies in the laminar flow region, with Reynolds Number less than 0.4.

The ash size distributions used for calculating ash deposit thicknesses in the work presented here were based on historic and well-preserved ash deposits, after making allowance for the fraction of very fine particles that do not fall out within the area of identifiable ash deposits. In New Zealand, we use three size distributions, based on the rock type of the volcano (Fig. 1 and Table 1). For rhyolitic volcanoes the distribution is based on the Hatepe Ash (Walker, 1981), for andesitic volcanoes it is based on the AD 79 eruption of Vesuvius (Macedonio et al., 1988), while for the basaltic volcanoes it is based on the 1973 eruption of Heimaey (Self et al., 1974). These ash size distributions are similar to those used for hazard analysis in Hurst and Smith (2010). Size-dependent densities are used to convert size to mass for each of the three ash types. The distribution shown as VAAC in Fig. 1 is derived from Heffter and Stunder (1993). The normal practice at MetService, who operate the Wellington VAAC, when using HYSPLIT for Aviation Hazard purposes, is to assume that 5% of any total eruption volume is comprised of fine ash, with this VAAC size distribution, and to ignore the 95% of coarser ash as it is falling out quickly.

Fig. 1

The three ash size distributions used by ASHFALL, for Andesitic, basaltic and Rhyolitic volcanoes, compared to the ash size distribution used for VAAC modelling of ash as an aviation hazard. On the same particle diameter scale are also shown particle fall speeds for the three equations of Bonadonna et al. (1998) and the Ganser equation

Aggregation of ash particles can have significant effects on ash distribution (Carey and Sigurdsson, 1982, Mastin et al., 2016), however this is not considered in either HYSPLIT or ASHFALL, and the focus of this work is the incorporation of realistic wind fields.

The larger (>100 micron) ash particles fall at higher terminal velocities, which produce turbulent air flow, with Reynolds Numbers greater than 500, and their actual terminal velocities are much less than that given by Stokes Law (Rose, 1993). Volcanologists typically use empirical methods or a range of equations for particles of different sizes, based on the Reynolds Number (Bonadonna et al., 1998). Figure 1 also shows how the three equations used for different ash sizes by Bonadonna et al. (1998) can all be matched by an equation derived by Ganser (1993) covering all particle sizes. This Ganser equation has recently been incorporated in HYSPLIT (Dare, 2015) and we have used this in our ash modelling.

It was also necessary to provide a more accurate initial ash concentration as a function of height, with a concentration based on a typical umbrella-shaped eruption column, rather than uniform across all levels. The distribution used was that suggested by Suzuki (1983), in which the ash concentration at height h is proportional to (1-h/H) e-^A(1-h/H), where H is the maximum column height, and A is the Suzuki constant. We used a constant of 4, commonly used for sub-plinian and plinian eruption columns (Pfeiffer et al, 2005), for which the distribution is plotted in Fig. 2.

Fig. 2

Suzuki distribution with constant of 4, representing how ash is concentrated near the top of an eruption plume as it slows down due to its decreasing temperature and a lower density atmosphere

Comparisons of HYSPLIT and ASHFALL showed broadly similar results on days when there were uniform wind patterns, but had large differences on days with more complex flows, as in the example of Fig. 3, for 5 August 2016. In this case, with a complex low pressure system over the central North Island, although there was a consistent north-westerly wind above Ruapehu, once ash emitted from Ruapehu had gone some distance to the east it would have encountered a completely different wind pattern, with much of it blown back towards the west.

Fig. 3

Example of ash distribution forecasts for a hypothetical short-duration Ruapehu eruption at 0500UT on 8 August 2016, with a 12.5 km column height and a total volume of 0.1 km³. Colours and grey 0.01 mm contour are from HYSPLIT, black 0.01, 0.1, 1 and 10 mm contours from ASHFALL. The ASHFALL distribution is entirely based on the expected wind at Ruapehu, whereas HYSPLIT forecasts that a significant amount of ash will be deposited to the south and west under the influence of a complex wind pattern

Accordingly, we have started the implementation of routine ash forecasts from HYSPLIT, using the same ten volcanoes, three sizes of eruptions, and two eruption column heights as we use for ASHFALL.

Ash distribution comparisons for Te Maari eruption of 6 August 2012

We did not compare a HYSPLIT model with the actual volcanic ash thicknesses from the 1995 and 1996 eruptions of Ruapehu, used to evaluate the performance of ASHFALL predictions by Hurst and Turner (1999), because of the difficulty in obtaining an accurate 4 dimensional wind data for that time. However, there was a small ash eruption from the Upper Te Maari crater of Mt Tongariro on 6 August 2012, and NWP data for this date was archived by MetService, making it comparatively easy to calculate a hindcast for this event with HYSPLIT and ASHFALL. This NWP data was from the operational 8 km resolution WRF-ARW model, using European Centre for Medium Range Weather Forecasting (ECMWF) Integrated Forecasting System (IFS) data for initial and boundary conditions.

This eruption, which followed a period of low-level seismicity (Hurst et al., 2014), was the first eruption of Mt Tongariro for over 75 years (Scott and Potter, 2014). It occurred at night, commencing at 11:52 pm NZST (11:52 UTC), so there were no visual observations, but weather radar observations showed an approximate ash column height of 7800 m above sea level (the altitude of Upper Te Maari is 1485 m), and a duration of less than 15 min (Crouch et al., 2014). Ash samples were collected along roads in the vicinity of the volcano, and downwind to the east, by researchers from Massey University (Pardo et al., 2014). The ash cloud was observed in both visible and infra-red satellite images taken during the early hours of 7 August (NZST), as shown in Fig. 4.

Fig. 4

A moonlight image of part of the North Island, New Zealand from the NASA-NOAA Suomi satellite using the VIIRS (Visual Infrared Imaging Radiometer Suite) showing the ash cloud travelling east from the Te Maari eruption (white triangle) taken at 0055 NZST on 7^th August 2012. NOAA photo

Because relevant atmospheric model data had been archived, it was comparatively easy to run HYSPLIT with the same information for the Te Maari eruption as would be available for an eruption today. The results have been compared with the main plume that travelled eastwards. Some ash travelled north from Te Maari, but it was clear (Lube et al., 2014; Pardo et al., 2014; Turner et al., 2014) that these ash deposits were produced by low-angle directed blasts.

The initial result from the HYSPLIT model with a point ash source had an extremely narrow ash plume heading east. With such a narrow plume we need to take account of the horizontal extent of the umbrella portion of the initial eruption column, as the width of this will significantly affect the near-source ash distribution. This has not been necessary with ASHFALL, because of the large horizontal diffusion coefficient it used. The GNS webcam south of Ruapehu, with an infra-red sensitive camera which observed the top portion of the Te Maari eruption column (Crouch et al., 2014), suggested a diameter of the initial eruption column of approximately 1.5 km, so we modelled the initial ash distribution as a vertical cylinder with this diameter, with uniform ash concentration horizontally, and a Suzuki distribution vertically. This technique was used by Crawford et al. (2016), but on a much larger scale, when modelling the 2008 eruption of Kasatochi. It was also apparent from the radar images in Crouch et al., (2014) that the eruption column was not vertical, but directed northwards, with its top region between one and two kilometres north of the source vent. We did not try to model a sloping column, but since most of the ash is in the top region of the column, we simply shifted the source one kilometre north of the Upper Te Maari crater.

Figure 5 shows a comparison of the actual ash distribution from (Pardo et al., 2014) along the approximately north-south State Highway 1 and nearby roads with the output from both ASHFALL and HYSPLIT for a north-south profile at 175.81° E (about 12 km east of Te Maari), using a point source for ASHFALL and a 1.5 km diameter cylindrical source for HYSPLIT. Both models used a source centred 1 km north of the Upper Te Maari vent, and had their maxima close to the actual maximum ash thickness. The width of the actual ash distribution was intermediate between the two models, but generally closer to the very narrow distribution predicted by HYSPLIT. This is confirmed by the narrow initial plume seen heading east from Te Maari in Fig. 4, which can also be seen in Fig. 6.

Fig. 5

North-south ash distribution profile from the Te Maari eruption of Tongariro at approximately 12 km downwind at 175.81° E, showing modelled thicknesses from ASHFALL and HYSPLIT, and actual ash thicknesses from Pardo et al. (2014), with maximum thicknesses normalised to 100

Fig. 6

Ash thickness observations from Te Maari eruption (black triangle), compared to forecasts from ASHFALL and HYSPLIT. ASHFALL contours are 10 μm, 100 μm and 1 mm. Ash observations from Pardo et al. (2014), with blue for trace amounts, and other colours approximately equivalent to thicknesses of 0.25 mm (dark green), 0.5 mm (light green), 1 mm (yellow) and 2 mm (orange)

The ash distribution predicted by ASHFALL for the Te Maari eruption has already been discussed by (Turner et al., 2014), which also compared both the ground and satellite measurements of where ash particles travelled with the results of the UK Met Office aerial dispersion model NAME-III (Jones, 2004), which used the NWP model NZLAM-12. The NAME and HYSPLIT models are very similar; both being Lagrangian dispersion models employed for a range of atmospheric pollution problems. NAME-III and HYSPLIT gave fairly similar ash distributions from this eruption, and in particular both had narrow ash plumes compared to ASHFALL. We will not discuss the differences between the ash distributions from these modelling programs further, other than to note that they probably are as much due to using different atmospheric models as to differences between the programs.

The region further downwind from this eruption was mountainous and virtually unpopulated until at least 80 km from Te Maari, and by then only trace amounts of volcanic ash were present, so all that could be said was that it was present or not present. Figure 6 shows which of the small settlements from 80 to 200 km reported either the presence or absence of ash. We have compared these reports to the outer limits of the ASHFALL and HYSPLIT predictions, effectively for an ash thickness of a few microns.

In Fig. 6 it can be seen that the narrow range of ash directions predicted by HYSPLIT does not cover all the places where ash was observed, suggesting that the wind was directed slightly more to the north than predicted by the atmospheric model. The available ground truth data in the far downwind area is too limited to let us distinguish between the narrow ash distribution of the HYSPLIT model and the much wider distribution from ASHFALL. However, the narrowness of the ash distribution in the near-source region does suggest that we have a narrow plume, but directed slightly north of where the atmospheric model used for HYSPLIT predicted it to be.

Discussion

The ash predictions produced by HYSPLIT from a 4-D atmospheric model use much more information about the atmospheric wind pattern, and so are potentially much more accurate than the ASHFALL model that is currently used which assumes a horizontally uniform wind profile. For the 6 August 2012 eruption from the Upper Te Maari crater of Tongariro, the HYSPLIT forecast using an archived atmospheric model agreed well with the actual ash distribution recorded about 12 km downwind. There was not such good agreement between HYSPLIT and the reports of where small amounts of ash were present in areas over 80 km downwind. As found by Hurst and Turner (1999) for the 1995 and 1996 eruptions of Ruapehu, the main reason for the difference between forecast and actual ash distributions was the difference between the forecast and actual wind pattern.

HYSPLIT should provide a much better match to the ash distribution of an actual eruption than ASHFALL, because of the more detailed atmospheric model, and because it more completely accounts for ash removal processes (e.g. wet deposition). However, when using these programs for forecasting purposes, there are two factors that need to be considered. Firstly, HYSPLIT outputs have more detail, compared to the simple ellipses produced by ASHFALL, so they give a higher impression of accuracy. Secondly, without the non-physically large horizontal diffusion coefficient used in ASHFALL, HYSPLIT sometimes forecasts narrow ash distributions, as was seen for the Te Maari eruption, for which a very small change in wind direction will displace the forecast ashfall area to the extent that it no longer overlaps with the actual ashfall area some distance from the vent. This difficulty in demonstrating value in more detailed forecasts is something already encountered in NWP for normal weather forecasts, where for instance high resolution forecasts of severe convection are seen subjectively as being more accurate, but traditional verification metrics strongly penalise displacement errors, which high resolution models by nature are more prone to (e.g. Mass et al 2002). Human forecasters (e.g. Roebber et. al. 2004) and/or the use of model ensembles (e.g. Peralta et. al., 2012) or other post-processing techniques (e.g. Theis et. al., 2005)) are needed to account for uncertainty in this high resolution model output, and appropriately convey risk to end users.

To provide more useful automated quantitative ashfall forecasts from a HYSPLIT/WRF system that account for uncertainty in the NWP winds and eruption parameters, we could follow the example of high resolution NWP and use an ensemble of model runs based on multiple NWP forecasts and a range of eruption parameters.

For airborne ash the Wellington VAAC does this to some extent by running a small HYSPLIT ensemble, with a range of eruption parameters, and with meteorological data from 3 NWP models. Dare et al., (2016) examined the impact of atmospheric uncertainty with a 24 member NWP model ensemble, again for airborne ash. Large ensemble systems are not yet used operationally by VAACs (WMO 2012, 2015). We believe similar ensemble based approaches that also include eruption source uncertainty are needed for ash fall forecasts.

An immediately achievable way of including uncertainty into HYSPLIT predictions would be to follow the Wellington VAAC approach and make the predictions based on multiple NWP models for each hypothetical eruption. Further investigations would be required to determine how best to present this output; options include the most conservative approach of outputting the union of ensemble ash fall areas, the contour of a specific percentile of ash depth (e.g. Draxler et al 2003), or the probability of exceedance for a specific ash depth.

By moving to the HYSPLIT/WRF we have eliminated the gross forecast errors that ASHFALL will produce in a complex wind field. The approach that the existing ASHFALL model takes to address meteorological uncertainty (an artificially large diffusion coefficient) apportions too much of this uncertainty close to the vent, and even upwind, whereas NWP wind errors should grow downwind from the vent.

It should also be noted that as well as the uncertainties of the wind pattern for ash forecasting, there will be large uncertainties in the eruption parameters. We cannot predict what eruption (i.e. magnitude, plume height) will next take place from any particular volcano, and when an eruption starts we often have very limited information.

Conclusions

We have successfully set up the use of HYSPLIT to give forecasts of volcanic ash deposition, were a New Zealand volcano to erupt. These forecasts use a much more detailed atmospheric model than the forecasts produced by ASHFALL, in particular the wind pattern varies in 3 spatial dimensions, rather than one. This gives potentially more accurate forecasts, particularly in situations with complex wind patterns.

Comparison of HYSPLIT hindcasts of the ash deposition from the 8 August 2012 Te Maari eruption from Mt Tongariro showed that HYSPLIT would have accurately forecast the narrow ash plume observed about 12 km downwind of the vent, but there was apparently a discrepancy between the forecast and actual wind direction further downwind, that significantly altered where ash was deposited. This emphasizes that we need to give an indication of what areas are at risk of volcanic ash, rather than just giving a single best estimate of the ash distribution.