DEVELOPMENT OF A VOLCANIC HAZARD MODEL FOR NEW ZEALAND: FIRST APPROACHES FROM THE METHODS OF PROBABILISTIC SEISMIC HAZARD ANALYSIS

We commence development of a volcanic hazard model for New Zealand by applying the wellestablished methods of probabilistic seismic hazard analysis to volcanoes. As part of this work we use seismologically-based methods to develop eruption volume frequency distributions for the Okataina and Taupo volcanoes of the central Taupo Volcanic Zone, New Zealand. Our procedure is to use the geologic and historical record of large eruptions (erupted magma volumes ~ 0.01 cubic km for Taupo and ~ 0.5 cubic km for Okataina) to construct eruption volume-frequency distributions for the two volcanoes. The two volcanoes show log-log distributions of decreasing frequency as a function of eruption volume, analogous to the shape of earthquake magnitude-frequency distributions constructed from seismicity catalogues. On the basis of these eruption volume-frequency distributions we estimate the maximum eruption volumes that Taupo and Okataina are capable of producing at probability levels of relevance to engineers and planners. We find that a maximum eruption volume of 0.1 cubic km is expected from Tau po with a I 0% probability in 50 years, while Okataina may not produce a large eruption at this probability level. However, at the more conservative 2% probability in 50 years, both volcanoes are expected to produce large eruptions (0.5 cubic km for Okataina and 1 cubic km for Taupo). Our study therefore shows significant differences in eruption probabilities for volcanoes in the same physiographic region, and therefore highlights the importance of establishing unique eruption databases for all volcanoes in a hazard analysis.


INTRODUCTION
In recent years there have been major developments in the assessment of geological hazards.One of the biggest advances has been in the development of probabilistic hazard analysis, which produces estimates of the probability of occurrence of a hazard at a site according to the location, severity and frequency of occurrence of hazardous events around the site.For example, probabilistic seismic hazard analysis (PSHA: Cornell, 1968) is now routinely used the world over (e.g.Frankel et al., 1995;Stirling et al., 1998 and2002).PSHA provides estimates of the probability of occurrence of a suite of ground motion levels at a site from data describing the location, magnitude and frequency of occurrence of earthquakes in the surrounding region, and from relationships that provide estimates of ground motion as a function of earthquake-to-site distance.The general 4-step procedure of PSHA is shown in Figure 1.While the approach of incorporating occurrence rate information into a hazard analysis makes it more complex than the traditional approach of defining deterministic or scenario-based hazard (i.e.hazard based on the closest and severest of hazardous events, irrespective of the frequency of occurrence of that event), the resulting hazard curve is usually more useful, since it provides information on the level of hazard expected in any given return period.Therefore, a single probabilistic hazard analysis can be useful for different end users who each need to know the expected level of hazard for a particular return period.For example, town planners usually need to know about the levels of earthquake shaking expected with a return period of about 500 years, whereas nuclear power plant operators need to know about the very rare and higher levels of shaking expected with a return period of 10,000 years.
In contrast, development of probabilistic volcanic hazard models is less well advanced.This may be due to the fact that probabilistic hazard analysis was largely developed in the arena of earthquake hazard, and application of the methodology to other hazards has been a more recent development.It could also be due to the inherent 'messiness' of volcanic phenomena, in that hazards posed by eruptions vary greatly with the style of activity (e.g.fall from a high plume versus generation of pyroclastic density currents, versus generation of lava flows or domes) and such styles can change rapidly, and affect markedly different areas (particularly due to wind changes affecting dispersal of a high eruption plume) during the course of a single event.We are working towards development of a comprehensive probabilistic volcanic hazard (PVH) model for New Zealand (in particular reflecting the most frequent and widespread hazard, that of tephra fall) comparable in scope to that already available for seismic events (e.g.Stirling et al., 2002).In this paper we adopt the 4-step procedure of PSHA, and undertake our analysis to address "STEP 2" of a 4-step Probabilistic Volcanic Hazard Analysis (PVHA; Fig. 2), equivalent to its seismic counterpart (Fig. 1).In other parts of the world, individual elements of the PVHA methodology shown in Figure 2 have been applied to forecasting the size and timing of future events at individual volcanoes (e.g.Bacon, 1982;Wadge, 1982), to assess the probability of a hazardous volcanic event at a given site from one or more 267 sources (e.g.Ho, 1991;Connor et al., 2000Connor et al., , 2001)), and to derive probabilistic forecasts of tephra fall from individual volcanoes (e.g.Rhoades et al., 2000).The elements or steps of PVHA are as follows (numbered according to the steps shown in Fig. 2

Log Peak Acceleration
Figure 1.The four-step procedure of probabilistic seismic hazard analysis, adapted from Cornell (1968).
(1).The volcanic sources have to be identified with respect to any given site.This in general is accomplished for New Zealand, as all active or potentially active areas of volcanism have been identified (e.g.Smith, 1986;Latter, 1986).
(2).Eruption volume -frequency relationships.This is the step we address in this paper, and is analogous to developing earthquake magnitude-frequency distributions in PSHA ("STEP 2" of Fig. 1).To date eruption volume -frequency distributions have only been established in New Zealand for some volcanic systems, and only for the relatively recent time periods that can be dated using radiocarbon techniques (e.g.Froggatt & Lowe, 1990).Records of earlier activity (pre-ca.65,000 years BP) are undoubtedly incomplete for the smaller eruptions, due both to erosion or burial of older deposits and a lack of detailed documentation.Thus the record from only the early large events (largely-silicic caldera-forming events in TVZ) has been dated (Houghton et al., 1995), and even then there is some uncertainty as to whether activity representative of the post-65,000 year record in TVZ was present as a background to the large events, or if the style(s) of activity has fundamentally changed.The data used in this step are taken from the post 26,500 year record for Taupo and Okataina, as available information indicate that the eruptive systems of both volcanoes were significantly different in terms of chemistry and eruptive patterns (Nairn, 1989;Jurado-Chichay and Walker, 2000).Though limiting our attention to Taupo and Okataina, our study also develops a methodology for addressing "STEP 2" at other volcanoes.(3) In PSHA, this is the most complex step, because of the differing controls of attenuation from source earthquakes and local ground response that will affect the shaking acceleration.With PVHA, this step is exceptionally complex (e.g.Blong, 1984), for two reasons.First, the styles of different kinds of activity vary widely and nonsystematically.For example, fall activity will deposit tephra in thicknesses that will diminish gradually away from source, but in direction(s) controlled by the wind; consequently the tephra fall thickness is a complex function of distance and azimuth from source.Significant work has now been done to understand and model tephra attenuation from New Zealand volcanoes (Rhoades et al., 2000).Second, certain kinds of volcanic activity, such as lava flows, and pyroclastic flows and lahars, will produce locally extreme effects, regardless of distance from source, but their distribution will be controlled to a greater or lesser extent by the local topographic relief (e.g.Wadge et al., 1994;Iverson et al., 1998).We do not portray these localised hazards in the PVHA diagram in Figure 2.

STEP2
(4) The final step is to construct a hazard curve, which in our example of PVHA is a graph showing the frequency or probability of exceedance for a suite of tephra fall thicknesses.Again, the hazard curve shown in "STEP 4" of Figure 2 does not represent the more localised hazards such as pyroclastic flows and lahars.The hazard curve shown can be used to estimate the maximum thickness of tephra fall that can be expected for any user-specified return period (e.g.c. 500 year return period).
In this paper we focus our attention on the Taupo and Okataina volcanoes in the central Taupo Volcanic Zone (TVZ) of New Zealand (Fig. 3).These volcanoes are responsible for the largest and most destructive eruptions in historic times (The 1886 AD Tarawera eruption from Okataina) and the last 25,000 years (Taupo 1,800 yrs BP) in New Zealand.Existing approaches to hazard assessment for these volcanoes have used a deterministic approach (e.g.Nairn et al., 1998), but the variability in eruption sizes and styles means that choices have to be made of eruption scenarios that may not be realistic for future events.While the main thrusts of this paper are in developing the framework of PVHA and addressing STEP 2 of PVHA for Taupo and Okataina volcanoes, we also use the eruption volume-frequency distributions to make first-order estimates of the maximum volume of eruption expected from the volcanoes with a 10% and 2% probability of exceedance in 50 years.In PSHA these are probabilities and time periods of interest to engineers and planners (e.g.Frankel et al. 1995;Stirling et al. 1998Stirling et al. , 2002)).While not constituting STEP 4 of Figure 2, these probabilistic estimates give an indication of the eruptions that should be of concern to engineers and planners located near the two volcanoes.

SETTING AND ERUPTIVE HISTORY OF TAUPO AND OKATAINA VOLCANOES
Taupo and Okatmna volcanoes are located in the central TVZ in the North Island of New Zealand (Fig. 3, Houghton et al., 1995).Within the central TVZ, volcanism has been overwhelmingly rhyolitic and concentrated at a number of caldera centres, of which Taupo and Okataina are the highly active representatives at the present day.
The eruptive history of the central TVZ is established on two contrasting time scales.Over a long-term time scale, a framework is avmlable that deals with the largest (calderaforming) eruptions, with volumes between about 30 and 1,000 cubic km, and a mean recurrence period of about 50,000 years (Houghton et al., 1995).Because of decreasing quality of exposure with age, the minimum eruption size captured by the primary subaerial record (Houghton et al., 1995) may be as large as 100 cubic km.Further information is potentially available from primary fall or secondary reworked tephras in sedimentary sequences onshore and offshore (e.g.Shane, pers comm.), but methods have yet to be established to identify sources and/or eruption volumes.Over a shorter-term time scale, activities over the past c.
As mentioned in the Introduction, the inbuilt complexity of volcanic data presents a potential problem for the application of the four-step procedure of PSHA to PVHA (Figs. 1 & 2).This is because a given volume of magma can be erupted in one or more of several different forms, such as lava flows, tephra fall deposits, pyroclastic flow deposits and subsequent lahars, each of which has a characteristic dispersal and impact.Thus the form which a future eruption will take, its products and hazards, for a given volume of magma cannot be reliably assessed, except in generalised terms.For instance, at Taupo, all eruptions greater than 2 cubic km in the last 25,000 years have developed pyroclastic flow deposits, though of greatly varying extents.Two contrasting treatments are therefore equally appropriate for quantifying the size of volcanic eruptions.One is the total volume of erupted magma, since this is a direct function of the productivity and timing of the magmatic system.The other is the volume of pyroclastic deposits (and subsequent reworked materials such as lahars), as these are the manifestations of the volcanism that represent the greatest hazard.Since our primary concern is to understand the relationship between eruption volume and frequency of occurrence ("STEP 2" of Fig. 2), our preference is to adopt the simplest approach and measure eruption volume by total volume.Volume!: Volume ofpyroclastic material demonstrated to exist.This is the measure of the fall (and some flow) deposits that produce the greatest impact.Volume2:The total magma volume of eruption.This is the eruption measure used for developing our eruption volume-frequency distributions.Note that for Taupo these are simply taken to first order to be the same as for the demonstrable volumes of tephra (see Wilson, 1993), whereas at Okataina values represent tephra reduced to magma volume plus the volumes of coeruptive lavas (Nairn, 1989) Repose Periods: The time intervals between the eruptions.Repose Period Before: The time period from the previous eruption to this eruption.Since this contains one more repose period than "Repose Period After", we use this for constructing our eruption volume-frequency distributions.Repose Period After: The time period from this eruption to the next eruption.Age: 40 Ar/ 39 Ar ages from Houghton et al. (1995).Data Sources for younger eruption record: Nairn (1989); Wilson (1993).* Denotes eruptions that are highly uncertain in terms of eruption volume and age, and we are certain that this part of the record of the Table is incomplete.

ERUPTION VOLUME • FREQUENCY DISTRIBUTIONS
Probabilistic volcanic hazard analysis requires information on the size and occurrence rate of eruptions ("STEP 2" of Fig. 2).We therefore use the Taupo and Okataina data (Table 1) to construct eruption volume-frequency distributions (the number of eruptions per year greater than or equal to a given eruption volume).In general terms the distributions are constructed by sorting the eruptions in Table 1 by volume, counting up the number of eruptions greater than or equal to each eruption volume (the cumulative number of eruptions), and then dividing these numbers by the total timespan of the eruption record to get the cumulative annual rate of occurrence for a particular eruption volume (rate at which an eruption volume is equalled or exceeded).This is the same method used to construct an earthquake magnitude-frequency distribution, in which the number of earthquakes in a catalogue of a particular magnitude or greater are counted up to give the cumulative number of earthquakes, and then the cumulative number is divided by the timespan represented by the catalogue to get the cumulative annual rate.
The calculation of cumulative annual rates for eruption rates for eruption volume-frequency distributions first requires a complete volcanic record for a given timespan, down to some minimum eruption volume.We therefore have to determine whether Table 1 has recorded all eruptions greater than the minimum eruption size listed for each volcano in Table 1, which is 0.01 cubic km for Taupo since 22,600 years B.P., and 0.5 cubic km for Okataina since 25,000 years B.P. We doubt that the record is complete for both volcanoes over these time periods, since the geologic record is unlikely to have preserved all of the smaller eruptions that have occurred in the timespan represented by Table 1.Note that we exclude the older large Okataina eruptions (i.e. the eruptions with an * at the base of Table 1) from our assessment of eruption completeness since we are certain that this part of the record is incomplete.Natural processes of erosion are likely to have removed the tephra deposits of some eruptions, and some deposits may have been obliterated by much larger, more recent eruptions.
To assess the level of completeness of the volcanic record of Table l we apply well-established methods for assessing the completeness of earthquake catalogues.Earthquake catalogues typically can be broken into "sub-catalogues", each having a complete record down to some minimum magnitude.For example the New Zealand catalogue has recorded all earthquakes of magnitude 4 and above since about 1964, magnitude 5 and above since about 1940 and magnitude 6.5 and above since about 1840 (e.g.Stirling et al. 2002).The reason for the progressive improvement in detection threshold is the progressive increase in human population across New Zealand (more written records) and the advent of a computerised seismic network in 1964.These three "sub-catalogues" are collectively used in seismic hazard assessments in New Zealand because they maximise the use of complete earthquake records for calculating occurrence rates.For example, calculating the rate of occurrence of magnitude 6.5 and above earthquakes that have occurred since 1840 will be a more reliable rate than the magnitude 6.5 rate derived from the much shorter post-1963 record.In light of the above we identify complete "subcatalogues" in Table 1 from graphs of the cumulative number of eruptions versus time (Fig. 4).With the assumption of time-invariant eruption rates (i.e. a steady occurrence of eruptions as a function of time), these plots should show a roughly constant slope for the eruption volumes and timeperiods associated with a "sub-catalogue".The graphs in Figure 4 show a sharp increase in the number of recorded eruptions at about 7,500 years B.P. for Taupo, and a steady rate of eruption occurrence since about 25,000 year B.P. for Okataina.We attribute the Taupo observations to an incomplete record for the earlier, smaller eruptions, rather than a marked acceleration in the rate of volcanism (the latter explanation is possible, but we consider it a less likely explanation than the former).
These "sub-catalogues" are defined by plotting various subsets of the data as shown in Figure 5.The plotted lines that show roughly constant slopes (i.e.constant rates of occurrence of eruptions as a function of time) are assumed to indicate a complete "sub-catalogue".In Figure 5 we show the two "sub-catalogues" of the Taupo data that we consider complete.Since the entire Okataina dataset (excluding unreliable pre c. 25,000 year data) plots a roughly straight line in Figure 4 then we assume that the eruption record is complete to the smallest eruptions given in Table 1 for Okataina (i.e. down to 0.5 cubic km).

Graphs of the cumulative number of eruptions versus time for the two subsets of eruption data for Taupo (Table 1 ).
See the text for further explanation.
The next step is to construct eruption volume-frequency distributions (Fig. 6) by combining the eruption rates calculated for the two Taupo subsets, and using the entire c. 25,000 year record for Okataina.Note that for illustrative purposes only we also show the pre-25,000 year Okataina data on Figure 6.This is simply to illustrate how the 25,000 year record might extrapolate out to larger, less frequent events.
The Taupo and Okataina eruption volume-frequency distributions are plotted in log-log scale, consistent with the construction of earthquake frequency distributions.Earthquake frequency distributions typically show the log of the number of events per year for a given magnitude, magnitude being a log-scale measure of the size of an earthquake.The two eruption volume-frequency distributions show the same general log-log trend of decreasing cumulative rate as a function of eruption volume.This is consistent with the shape of eruption size-frequency distributions developed for other volcanoes (e.g.Pyle, 1998).
The Taupo and Okataina eruption volume-frequency distributions appear quite different in the total range of eruption sizes, and in the general slope and shape of the distributions.The Taupo distribution shows a lower limit of eruption size that is ten times lower than the Okataina distribution on Figure 6.Furthermore, the Taupo distribution has a slope (b) of about 0.5 in an equation of the form logN/yr = a-b(logV) (N/yr is the cumulative annual rate of eruptions greater than or equal to eruption volume V), whereas the Okataina distribution has a slope (b) of about 0.9.Since the distributions are log-log, these different slopes translate to large differences in the rates of eruptions at the two volcanoes.Our particular choice of equation is analogous to that of the seismologically-based Gutenberg-Richter relationship logN/yr = a-bM.The Gutenberg-Richter relationship is used to describe the relationship between magnitude (M) and the number of earthquakes of magnitude M or larger (N/yr) for regions, countries, and even the globe (Gutenberg and Richter, 1944).The Okataina distribution also shows a "roll off' or shallowing of the slope of the distribution at volumes less than about 5 cubic km, which in seismology is usually attributed to incompleteness of an earthquake record for the smallest earthquakes.However, in this case there are three reasons why we infer this roll off to be real.
First, the proximal tephra exposures around Okataina are exceptionally well exposed, and intensive studies (Nairn, 1981(Nairn, , 1989) ) have failed, in contrast to Taupo (Wilson,I 993 ), to reveal any trace of smaller eruptions than those recorded in Table l.Second, a detailed record from a peat bog down-wind of Okataina (Lowe et al., 1999) records all previously documented eruptions 2: 0.2 km 3 from Taupo and Okataina over the past -18,000 years, and implies that no eruption of that magnitude or greater would have escaped detection.

ERUPTION PROBABILITIES
Assuming that the shapes of the Taupo and Okataina eruption volume -frequency distributions (Fig. 6) represent the real size distribution of eruptions at the two volcanoes, we can use the distributions to estimate the largest eruptions that can be expected with probabilities of interest to engineers and planners.We therefore estimate the maximum eruption volumes that can be expected with 10% and 2% probabilities of exceedance in 50 years.We make these estimates by assuming a Poisson model (e.g.Ho, 1991), in which the time since the last eruption is not considered, and the production of eruptions through time is assumed to be constant.Use of the Poisson model is in keeping with the assumptions embodied in our assessment of eruption record completeness in Figures 4 and 5.The following equation is used to calculate the maximum eruption volumes expected from by the two volcanoes according to a Poisson model: (1) In which P is the probability of an eruption of a certain size or greater occurring in time T, and r is the cumulative annual rate of that eruption size.If P is equal to 0.1 (10% probability) and T is equal to 50 years, r works out to be about 0.0021 events per year (return period of 475 years) in equation 1.If we examine Figure 6 we find that the eruption volumes that correspond to rates of 0.0021 are about 0.1 cubic km for Taupo and 0 for Okataina.With P equal to 0.02 (2% probability) and T equal to 50 years, the value of r is about 0.0004 (return period of 2,475 years).This rate corresponds to an eruption volume of about 1 cubic km for Taupo, and 0.5 to 5 cubic km for Okataina.At face value then, Figure 6 shows that Taupo volcano is expected to produce hazardous eruptions while Okataina may lie dormant at the 10% probability level in 50 years, whereas both volcanoes are expected to produce hazardous eruptions at the more conservative 2% probability level in 50 years.
We emphasise that our analysis has been kept at a relatively simple level.For instance, we have not examined the uncertainties in the calculated rates of eruption volumes in Figure 6, and addressing these uncertainties could influence the comparisons of the expected eruption volumes from the two volcanoes.Treatment of uncertainties is outside the scope of our national-scale probabilistic volcanic hazard analysis for New Zealand.Also, we have made a simple assumption that eruption production at the two volcanoes is a time-independent process, since the data do not show any clear indications of acceleration or deceleration in eruption rates with time.Accumulation of new data may test this assumption in the future.Lastly, we have simplistically characterised eruption volume as total volume of erupted material, but in the context of PVHA in New Zealand we will need to meld our approach to addressing "STEP 2" of PVHA (Fig. 2) with the results of recent studies that have addressed "STEP 3" (Rhoades et al. 2000).Rhoades et al. have restricted their attention to the pyroclastic portion of erupted material.Ongoing efforts will be in addressing the above for the Taupo and Okataina volcanoes, and in developing eruption volume-frequency distributions for the other active volcanoes in the TVZ and elsewhere.

DISCUSSION AND CONCLUSIONS
We have commenced development of a volcanic hazard model for New Zealand by applying the well-established methods of probabilistic seismic hazard analysis to volcanoes.As part of this work we have used seismologically-based methods to develop eruption volumefrequency distributions for the Okataina and Taupo volcanoes.We have then used these distributions to estimate the maximum eruption volumes that the volcanoes will produce with 10% and 2% probability in 50 years.From our results it appears that the two volcanoes show log-log distributions of decreasing eruption frequency as a function of eruption volume.Taupo is the only volcano of the two that is expected to produce a hazardous eruption with a 10% probability in 50 years, but the two volcanoes are expected to produce hazardous eruptions at the more conservative 2% probability in 50 years.Our study also shows that the methods of probabilistic hazard analysis can be readily applied to New Zealand volcanoes, and also shows the importance of establishing unique eruption databases for all volcanoes, despite them being close together within physiographic regions like the TVZ.
We have not entered into the well-studied realm of timepredictable or volume-predictable eruption statistics in our analysis (e.g.Bacon, 1982;Froggatt, 1982;Wadge, 1982).Embodied in these time or volume-predictable methods is the assumption that there is a predictive relationship between eruption volume and repose period.The typical approach of these studies is to plot cumulative eruption volume versus cumulative time, and then use the resulting curve to forecast the timing and size of future eruptions.We have not followed this approach here because: (a) it is only our intention to consider Poissonian or time-independent eruption probabilities for our New Zealand volcanic hazard model, since this is consistent with our characterisation of other hazards at a national scale (e.g. the national seismic hazard model; Stirling et al. 2002); and (b) our own analysis of the available data show no such predictive relationships to exist.To illustrate, we plot eruption records from the activity at Coso, California (used because it is one of the best eruption records available for examining eruption timing and volume; Bacon, 1982), together with the data used here from Okataina and Taupo (Figure 7), to demonstrate that there is a seemingly random relationship between repose periods either before of after given eruptions and the size of that eruption.o 1000 2000 3000 4000 5000 6000 repose period before eruption, years repose period after eruption, years Plots of eruption volume versus the repose period before (A, C, E) and after (B, D, F) eruptions for the rhyolite eruptions at Coso (A, B: data from Bacon, 1982), Taupo (C, D:from Wilson, 1993), andOkataina (E, F: data from Nairn, 1989).

Figure 2 .
Figure 2. Our four-step procedure of probabilistic volcanic hazard analysis, adaptedfrom Figure I.

Figure 3 .
Figure 3. Map showing the location of the Taupo Volcanic Zone, and the Taupo and Okataina volcanoes (outlined).

Figure 4 .
Figure 4. Graphs of the cumulative number of eruptions versus time for Taupo (A) and Okataina (B) volcanoes.The data arefrom Table1.
Figure 5. Graphs of the cumulative number of eruptions versus time for the two subsets of eruption data for Taupo (Table 1 ).See the text for further explanation.

Figure 6 .
Figure 6.Eruption volume-frequency distributions for Taupo and Okataina volcanoes.Data points surrounded by parentheses are considered to be unreliable in terms of volume and frequency of occurrence (i.e.not an output from our statistical analysis of the eruption record), and are shown for illustrative purposes only.See the text for further explanation.