xml version 1.0 encoding UTF8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd
leader nam 2200421Ka 4500
controlfield tag 006 m d s
007 cr mnu
008 041209s2004 flua sbm s0000 eng d
datafield ind1 8 ind2 024
subfield code a E14SFE0000485
035
(OCoLC)57709174
9
AJU6859
b SE
040
FHM
c FHM
090
QE26.2 (ONLINE)
1 100
Weller, Jennifer N.
0 245
Bayesian inference in forecasting volcanic hazards
h [electronic resource] :
an example from Armenia /
by Jennifer N. Weller.
260
[Tampa, Fla.] :
University of South Florida,
2004.
502
Thesis (M.S.)University of South Florida, 2004.
504
Includes bibliographical references.
516
Text (Electronic thesis) in PDF format.
538
System requirements: World Wide Web browser and PDF reader.
Mode of access: World Wide Web.
500
Title from PDF of title page.
Document formatted into pages; contains 60 pages.
520
ABSTRACT: Scientists worldwide are increasingly faced with the need to assess geologic hazards for very infrequent events that have high consequence, for instance, in siting nuclear facilities for volcanic hazards. One of the methods currently being developed for such assessments is the Bayesian method. This paper outlines the Bayesian technique by focusing on the volcanic hazard assessment for the Armenia Nuclear Power Plant, (ANPP), which is located in a Quaternary volcanic field. The Bayesian method presented in this paper relies on the development of a probabilistic model based on the spatial distribution of past volcanic events and a geologic process model. To develop the probabilistic model a bivariate Gaussian kernel function is used to forecast probabilities based on estimates of &lambda t, temporal recurrence rate and &lambda s, spatial density.Shortcomings often cited in such purely probabilistic assessments are that it takes into account only known features and the event, new volcano formation, is rare and there is no opportunity for repeated experiments or uniform observations, the hallmarks of classical probability. One approach to improving such probabilistic models is to incorporate related geological data that reflect controls on vent distribution and would improve the accuracy of subsequent models. Geophysical data indicate that volcanism in Armenia is closely linked to crustal movement along major right lateral strikeslip fault systems that generates transtension across region. The surface expression of this transtension is pullapart basins, filled with thick deposits of sediment, and antithetic normal faults. Volcanism in Armenia is concentrated in these deep sedimentary basins as is reflected in regional gravity data surveys.This means that low gravity anomalies are likely good indicators of future volcanic activity and therefore would improve probabilistic hazard models. Therefore, gravity data are transformed into a likelihood function and combined with the original probability model in quantitative fashion using Bayesian statistics. The result is a model that is based on the distribution of past events but modified to include pertinent geologic information. Using the Bayesian approach in this example increases the uncertainty, or range in probability, which reflects how well we actually know our probability estimate. Therefore, we feel it is appropriate to consider a range in probabilities for volcanic disruption of the ANPP, 14 x 10 per year (t=1 yr). We note that these values exceed the current International Atomic Energy Agency standard, 1 x 10 per year by at least one order of magnitude.
590
Adviser: Connor, Charles B.
653
Baye's Theorem.
ANPP.
kernel function.
spatial density.
gravity.
690
Dissertations, Academic
z USF
x Geology
Masters.
773
t USF Electronic Theses and Dissertations.
949
FTS
SFERS
ETD
QE26.2 (ONLINE)
nkt 1/5/04
4 856
u http://digital.lib.usf.edu/?e14.485
PAGE 1
Bayesian Inference In Forecasting Volcanic Hazards: An Example From Armenia by Jennifer N. Weller A thesis submitted in partial fulfillment of the requirement for the degree of Master of Science Department of Geology College of Arts and Sciences University of South Florida Major Professor: Charles B. Connor, Ph.D. Sarah Kruse, Ph.D. Robert Watts, Ph.D. Date of Approval: November 9, 2004 Keywords: baye's theorem, anpp, kernel function, spatial denisty, gravity Copyright 2004, Jennifer N. Weller ii
PAGE 2
Dedication To Chuck: Undertake something that is difficult; it will do you good. Unless you try to do something beyond what you have already mastered, you will never grow. Osborn ii
PAGE 3
Acknowledgements I would like to thank the Armenian National Academy of Sciences for excellent field and collaborative support. This work is supported by a NATO research grand and the USF foundation. Maps were created using GMT (Wessel and Smith). In addition, I would like to thank my advisor, Chuck Connor, for his continued time, help and support in all aspects of this thesis and for constantly reminding me that I was capable of working through it. Thanks also to my committee members, Sarah and Rob, for their editorial comments and time. Furthermore, I would like to thank the following people for their immeasurable support, without which this thesis would not have been possible: Tim, Patty, Don, and Jamie. iii
PAGE 4
Table of Contents List of Figures i Abstract ii Chapter One: Introduction 1 Chapter Two: Methods 7 Defining the Database 7 Developing the Probability Model 9 Developing the Gravity Model 16 Bayesian Model 19 Chapter Three: Discussion 23 Chapter Four: Conclusions 28 References 30 Appendices 34 Appendix A: Armenia Volcano Location Data 35 Appendix B: Computer Codes 39 ii i
PAGE 5
ii List of Figures Figure 1. Location Map 3 Figure 2. Photograph of ANPP 4 Figure 3. Gravity Map 6 Figure 4. Bayesian Flowchart 8 Figure 5. Data Model 15 Figure 6. Recurrence Rate Map 16 Figure 7. Smoothing Parameter vs. Probability 17 Figure 8. KolmogorovSmirnov Plot 18 Figure 9. Bayesian Map 23 Figure 10. Stylized Models 25
PAGE 6
iii Bayesian inference in forecasting volcanic hazards : An example from Armenia Jennifer N. Weller ABSTRACT Scientists worldwide are increasingly f aced with the need to assess geologic hazards for very infrequent events that have high consequence, for instance, in siting nuclear facilities for volcanic hazards. One of the methods currently being developed for such assessments is the Bayesian method. Th is paper outlines the Bayesian technique by focusing on the volcanic hazard assessment for the Armenia Nuclear Power Plant, (ANPP), which is located in a Quaternary vol canic field. The Bayesian method presented in this paper relies on the development of a probabilistic model based on the spatial distribution of past volcanic events and a geologic process model. To develop the probabilistic model a bivari ate Gaussian kernel function is used to forecast probabilities based on estimates of t temporal recu rrence rate, and s spatial density. Shortcomings often cited in such purely probabilistic assessments are that it takes into account only known features and the event, new volcano formation, is rare and there is no opportunity for re peated experiments or uniform observations, the hallmarks of classical probability. One approach to improving such probabilistic models is to
PAGE 7
iv incorporate related geological data that reflect controls on vent distribution and would improve the accuracy of subsequent models. Geophysical data indicate that volcanism in Armenia is closely linked to crustal movement along major right lateral strikeslip fault systems that generates transtension across region. The surface expression of this transtension is pullapar t basins, filled with thick deposits of sediment, and antithetic normal faults. Volcanism in Armenia is concentrated in these deep sedimentary basins as is reflected in regional gravity data surveys. This means that low gravity anom alies are likely good indicators of future volcanic activity and therefore would improve probabilistic hazard models. Therefore, gravity data are transformed into a likeli hood function and combined with the original probability model in quantitative fashion using Bayesian statistics. The result is a model that is based on the distribution of past events but modified to include pertinent geologic information. Using the Bayesian approach in this example increases the uncertainty, or range in probability, which reflects how well we actually know our pr obability estimate. Therefore, we feel it is appr opriate to consider a range in probabilities for volcanic disruption of the ANPP, 14 x 10 6 per year ( t =1 yr). We note that these values exceed the current International Atomic Energy Agency standard, 1 x 10 7 per year by at least one order of magnitude.
PAGE 8
1 Chapter 1 Introduction Scientists worldwide are increasingly faced with the need to assess hazards associated with pointlike f eatures such as volcanoes and earthquake epicenters on various temporal and spatial scales. Commonality among these phenomena exists because the analysis of their distribution a nd geologic setting can be used to estimate hazards quantitatively. Often, these geologic hazard assessments must evaluate the likelihood of very infrequent events that have high consequences (Haneberg 2000). For example, in the last two decades longter m probabilistic volcanic hazard assessment has increasingly been used in siting nuclear facilities worldwide (Crowe et al. 1982; Stamatakos and Ferrill 1996; Connor et al. 2000; McBirney et al. 2003; McBirney and Godoy 2003; Martin et al. 2004). Often, the central issue in these assessments is the likelihood of a new volcano forming by eruptions in close proximity to the facility. At such facilities, haza rds with probabilitie s on the order of 10 6 10 8 per year are often considered high (Connor et al. 1995, Martin et al. 2004) because overall the risks associated with such facilities must be very low. Geological hazard assessments for pointlike features should present robust estimates of hazard rates, based on the freque ncy of past events and insights about the geological processes that control such events One challenge associated with longterm
PAGE 9
2 probabilistic assessment of future volcanism is th at models of volcanic processes, such as the generation and ascent of magma, are inhe rently uncertain. One approach to making hazard assessments based on such models more robust is to modify probabilistic analyses by incorporating additional datasets th rough Bayesian inference (Von Mises 1957; Connor et al. 2000; Martin et al. 2004). Esse ntially, Bayesian inference allows us to combine two or more states of informati on (e.g., geophysical) to forecast the probability of volcanic events, such as formation of a new volcano, based on our understanding of volcanic systems, rather than solely based on the limited, and often in complete, record of volcanic events. If we consider the freque ncy of volcanic events to be a physical property of a magmatic system, we are faced with the conclusion that the limiting value of the frequency of volcanic events is unknow n. The event, formation of a new volcano, is rare and there is simply no opportunity for repeated experiments or uniform observations, the hallmarks of classical proba bility. Consequently, we are forced to update hazard forecasts using disparate obser vations of geologic and/or geophysical data that we believe impacts hazard forecasts. Bayesian inference provides a practical approach to incorporating such information. In this paper, we analyze volcanological and geophysical data from Armenia with the goal of calculating the hazard associated with the disruption of the Armenian Nuclear Power Plant (ANPP) (Karakhanian et al. 2003), outline the technique, and illustrate the problems inherent to such analyses. We do this through the construction of an improved model that focuses on the proba bility of renewed volcanism that would impact the ANPP, by combining the probabilistic and geophysical models using Bayesian inference. The
PAGE 10
3 ANPP ANPP is located in the northwestern part of the Ararat Depression (a sedimentary basin between Mt. Aragats in the north and Mt. Ararat in the south) in close proximity to the town of Metzamor, and 28 km west of Yerevan, the capital of Armenia (Figure 1). The ANPP is a Chernobylstyle reactor that sits at the base of the southern foothills of Mt. Aragats, the largest composite volcano in Armenia. Mount Ararat, another large composite volcano in Turkey, is 55 km south of the ANPP. The ANPP is located on the Shamiram Volcanic Plateau and is only 1.3 6 km south of 38 small cinder cones arranged in four local clusters (Figure 2) (Karakhanian et al. 2003). An additional source of volcanic hazard for the ANPP and the capital city of Yerevan are the volcanoes of the Ghegam Ridge located 52 km to the east of the site and just west of Figure 1. Location map of Armenia showing 554 Quaternary volcanoes used in study and faults. The rate of convergence just north of the ANPP is 1819 mm/yr based on REVEL 2000 models. The ANPP is shown just south of a cluster of 38 cinder cones.
PAGE 11
4n Lake Sevan (Figure 1). Some of these volcanoes have been dated as Holocene and their Late Pleistocene valley flow terminates 25 km east of the plant site. The most recent volcanic eruptions othe Ghegam Ridge have been dated between 4500 to 4400 + 70 yr BP (Karakhanian et al. 2003). Figure 2. Photograph of the ANPP clearly seen in the background are a cluster of cinder cones and to the far left the base of Mt. Aragats. Armenia is an appropriate choice for this type of analysis due to both its volcanic and tectonic setting. In the Quaternary (1.6 million years to the present), 554 basaltic to andesitic cinder cones (Savov et al. 2003) developed in response to mostly monogenetic activity. Monogenetic activity is characterized by the formation of a new volcano, such as a cinder cone or lava dome, and duration of volcanic activity at monogenetic volcanoes is thought to be typically less than 100 years (Connor and Conway 2000). After cessation of eruptive activity at any individual monogenetic volcano, renewed volcanism in the area builds a new monogenetic volcano. Thus, for this type of volcanism, the number of volcanoes reflects the number of volcanic events for which probabilistic forecasts are made. Because of the nature of this volcanic activity, the
PAGE 12
5 volcanoes themselves can be considered point like features (Figur e 1) and the hazard assessment reduces to the problem of estimating the density distribution of these features. Examples of other hazard assessments for m onogenetic volcanic activity include nuclear power plants and storage facili ties (Connor et al. 1995; Karakhanian et al. 2003; Martin et al. 2004) and urban centers such as Auck land, New Zealand (Magill et al. 2004) and Mexico City (Bloomfield 1975; Martin del Pozzo 1982). This distributed, monogenetic volcanism re sults from the complex tectonic history of the region that lies within a broad zone of deformation th at forms part of the AlpineHimalayan collision belt. Overall, volcanism describes an arc across Armenia (Figure 1) that is subparallel to this collision belt. Pullapart basi ns can be delineated by mapping anomalies in the Earths gravity field, cause d by density variations between the sediments filling the basins and the surrounding crust (T suboi 1979). Presently, as a part of the Alpine fold belt, the uplift occurring across Armenia is a result of the northward motion of the Arabian plate with respect to Eurasia (Philip et al. 2001). Th e rate of convergence of these two plates is 1819 mm/yr based on REVEL 2000 models (Dixon and Mao 2002). Volcanism across the region is linked to subduction and subsequent collision, and may result from slab steepening and breako ff which provides a viable mechanism for magma generation (Keskin 2003). In any cas e, volcanism is closely linked to NS compression and EW extension (Philip et al. 2001). The main geologic structures produced in this tectonic sett ing are northwest trending righ tlateral strikeslip faults. These faults produce areas of transtension th at create pullapart basins within which volcanism is localized. Contrasts in crustal structure reflected in the distribution of
PAGE 13
6nd en n ne n in analysis of volcano distribution and preparation of volcanic hazard forecasts. gravity anomalies have been shown to influence both the generation and ascent of magma (Tamura et al. 2002) and monogenetic volcanism has been shown to preferentially occur in pullapart basins and similar areas undergoing crustal extension (Connor aConway 2000). This correlation betwevolcano distributioand gravity anomalies occurs in Armenia as well (Figure 3), and may result from magmageneration by decompression melting of mantle previously enriched by subduction zoprocesses (Pearce 1990; Keskin 2003; Savov et al. 2003). Such a positive correlatiobetween volcanism and gravity anomalies make it appropriate to consider gravity data Figure 4. Gravity distribution across Armenia, coordinates are in UTM. Volcanoes typically occur in regions of broad gravity lows which correspond to deep sedimentary basins formed from extension and pullapart mechanisms. The above observations point to a strong correlation between patterns of volcanism and structure.
PAGE 14
7 Chapter Two Methods In this paper we outline a method that al that permits systematic analysis of haza f interest. The simplest way to visualize the process is in chart form (Figure 4). The flowchart provides a stepbystep summary of the Bayesian process from the in itial step of defining databases to the final step of evaluating hazard at th e location of interest. Each element of the flowchart relies on information about the region to co nstruct a reasonable hazard model. create them. In this paper we use two data sets. The first dataset consists of the geographic coordinates of Armenian volcanoes. These coordina tes were provided by the Armenia National Academy of Sciences and list all known Quaternary volcanoes in Armenia. The coordinates are provided in the universal transverse mercator (UTM) grid, WGS84 datum. The second dataset comprises gravit y data (Ohanissyan, 1985), also provided by the Armenia National Academy of Sciences. Th ese gravity data consist of measurements of the relative change in the ear ths gravity field, corrected fo r topographic effects. They lows geologic data to be cast in a way rd at a site o Defining the Database Volcanic hazard assessments are only as good as the data used to
PAGE 15
8 Figure 4: Flowchart of Bayesian method.
PAGE 16
9follow somewhat different paths in the development of a probability model and a geophysical model until they are transformed so they can be combined using Bayesian methods (Figure 3). We will describe both these paths beginning with the development of the probability model. Developing the Probability Model The first and perhaps most important step in developing a probability model for volcanic hazard is to explicitly define events in the region of interest (Connor and Hill, 2000; Martin et al. 2004). Given the complexity of volcanic processes, what event is specifically forecast by the probability model? This is important for a number of reasons, not least of which is that the entire development of the probability analysis depends on the event being defined as a point process. Second, having defined the event as a point process, it is important to unambiguously define what constitutes an event and what does not. For instance, an event could be the epicenter of an earthquake, or the location of a sinkhole, or any other process that can be isolated to a precise location spatially. In the consist of geographic location, in the same UTM projection used in the volcano dataset, and relative change in gravity (mGal). Following this compilation, the two datasets hazard assessments presented in this study, the definition of an event is limited to the formation of a new volcano, and is estimated directly from the distribution of existing Quaternary volcanoes and gravity anomalies. This assessment specifically does not looat the spatial or temporal distribution of volcanic eruptions at existing volcanoes or at tvolcanic hazards associated with eruptions, such as lava flows, lahars, or pyroclastic flows. This is because the hazard assessment is intended to assess the probability of k he
PAGE 17
eruptions from new volcanic vents in the area around the ANPP that would have a potential deleterious impact on operation of the nuclear facility, and containment of radionuclides (Crowe et al. 1982; McBirney and Godoy 2003). 10The volcano dataset used in this study includes both monogenetic volcanoes (one eruptive sequence only) and polygenetic volcanoes (volcanoes that have more than one eruptive sequence). Volcano types include cinder cones, domes, composite volcanoes, maars, and calderas. Individual vents are treated equally in the analysis, regardless of volcano type, because we are modeling the probability of formation of a new volcanic vent. Estimates of the probability of dike injection, without volcano formation, are not considered in this hazard assessment but may be important in other types of hazard assessments, such as for high level waste repositories (Woods et al., 1999; Connor and Conway 2000). Formally, each volcanic event in this assessment is considered to be independent of the other volcanic events in the data set. e l. in ets Alternative methods for defining volcanic events can be used in probabilistic volcanic hazard assessments. For example, polygenetic volcanoes may be weighted morheavily depending on the number of past eruptions (Martin et al. 2004). Also, closely spaced and similarly aged vents can be grouped together as a single event (Connor et a2000). Furthermore, only cones younger than a specific age, for example those formed inthe last 100,000 years, might be included in the analysis as volcanic events or these could be weighted more heavily. Because vents may be grouped into single volcanic eventsvarying ways depending on their timing, distribution, and episodes of activity, data s
PAGE 18
11igure by alternative definitions does not yet exist for Armenia. t P (1) er r time period for the expected operation of the ANPP. Once these s can be defined in different ways which would change the probability of an event occurring. In this case, we use a simple definition in which each mapped volcano (F1a) is a single event, partly because the additional data required Once the event is explicitly defined, we move on to the mathematical development of the probability model (Figure 3). Ultimately, the probability of an evenoccurring at a specific location is given by where A is the effective area of interest (e.g., the area within which volcanism wouldimpact the ANPP, should it occur). In this study we define A as 50 km AtsteN1]1[ 2 around the powplant. The time period of interest is given by the parameter t, (e.g., the duration of expected operation of the ANPP). We estimate a 100 yea ite specific parameters are specified, we are left with the problem of how to estimate s spatial density (events per kilometer) and t temporal recurrence rate (events per year). For Armenian volcanism, we use a simple estimate for temporal recurrence rate, ytttN (2) 01
PAGE 19
12 n, N is ion years and ty to 0, or present. This gives t a value of 3.5 x 10 f volcanic mating t have been employed (Ho 1991a,b; Condit and Connor 1996). kernel estimation technique that has been used in several studies to estimate spatial density of volcanism, including in the Pi nacate Volcanic Field, Mexico (Lutz and Tohoku region, Japan (Martin et al. 2004). In this technique, spatial variation in sfunction of distance to n earestneighbor volcanoes and a smoothing parameter, h. The kernel spreads probability away from the event (Diggle 1985). Differe nt kernel functions can be used including the Cauchy kernel (Martin et al and Gutmann 1994), and the Gaussian kernel (C onnor and Hill 1995). It is widely agreed that the shape of kernel function chosen in th is type of analysis generally has a trivial because the ages of individual volcanoes ar e not well constrained. In this equatio the total number of volcanoes, t 0 is the age of oldest event, and t y is the age of the youngest event. We know that all the volcanoes used in this study are all of Quaternary age so we set t 0 equal to 1.6 mill 4 years. That is, we expect one event every 2900 years. Wher e ages o events are better constrained, more detailed analyses have revealed time trends in volcanic activity and other methods for esti Compared with t estimation of spatial density, s is more difficult. However pioneering work by Diggle (1985) and Silverma n (1986) has led to the development of a Gutmann 1994), the Yucca Mountain region, USA (Connor and Hill 1995), and the is a function is a probability density function that is sy mmetric about the origin and 2004), the Epanechnikov kernel (Lutz
PAGE 20
13oblem clustering of vents. The probability surfaces generated by this model are continuous, as opposed to consisting of abrupt changes in probability that must be introduced in spatially homogeneous models (Connor and Hill 1995). Continuous probability surfaces ution. The bivariate Gaussian kernel is given by: impact on probability calculations compared to other parameters (Connor and Hill 1995; Lutz and Gutmann 1994). The Gaussian function was chosen for the Armenian prbecause volcanoes are treated as discrete events in time and space and the Gaussian model responds well to the patterns generally recognized in volcano distributions, such asallow for relatively easy comparison to other empirical data sets, e.g. gravity data that shed light on volcano distrib 2211),(hdieyx 22sNh (3) where d i is the distance from the point x,y, where s is estimated, to the i th volcano location, h is a smoothing parameter, and N is the number of volcanoes (points) thaused to estimate t are l of probability density function. Probability estimates made using equations 1 and 3 depend e s (x,y). Due to the fact that N occurs in the denominator, the integra s (x,y) across the map will be unity. Therefore the spatial density, s is a bivariate on the value chosen for h. Using a bivariate Gaussian kernel, events will have a high estimated probability in proximity to existing volcanoes if the value chosen for h is small, but low estimated probability away from the volcano. On the other hand, a large value of h will yield a more uniform estimate of probability distribution across the region. In thGaussian kernel, the smoothing factor is equivalent to the standard deviation of a
PAGE 21
14 function chosen. This model gives an estimate of spatial density based on positions of volcanic vents for Armenia (Figure 5). The kernel model is useful becau can be made allowing ease of comparison with other geologic informa tion; (2) there is no a at symmetric, bivariate Gaussian di stribution. Therefore, the kernel function depends on the assumption that the smoothing parameter is estimated in a geologically and/or statistically significant way (Connor and Hill 1995). In this study a wide range of smoothing constants are used (Figure 4); in addition, the range of reasonable smoothing constants is further constrained by use of th e ClarkEvans (CE) spatial cluster analysis (Clark and Evans 1954). This analysis s hows that volcanic events in Armenia are clustered across a variety of scales with >99% confidence (Blyth and Ripley 1980; Cressie 1991). Applying these tools, we use h=3000 m. As mentioned above, an additional assumption in the kernel estimation technique is the shape of kernel se (1) probability maps need to define zones of volcanic activity, as is required in a homogenous approach (Margulies et al. 1992); and (3) uncertainty in the distribution of individual events is easy to assess. The final step in developing the probability model is to prepare a contour map of spatial density (Figure 6). The map is made using a UTM projection for equal are measurements and the units of probability are typically conve rted into their logarithms because volcano vent density commonly varies by orders of magnitude across a region Also, testing the suitability of the smoothing pa rameter chosen is appropriate at this point (Figure 6). The maximum probability of an event occurring is 3.4 x 10 4 where h has a value of 4000 m. More importantly, this plot demonstrates that probability is sensitive to the value of h chosen and varies rapidly from a maximum probability of 3.4 x 10 4
PAGE 22
15 Figure 5. Modeling the data. Complimentary cumulative distribution function showing distance of nearestneighbor volcano versus its cumulative probability of occurrence for the original data set and various values of smoothing parameter h. Notice that while h=1500 appears to fit th e data well overall, it is a very poor fit in the tail of the dist ribution. Conversely, h=3000 fits the data we ll in the tail of the distribution.
PAGE 23
16 F igure 6. Map of spatial recurrence rate, coordinates are in UTM. Quaternary ent density is greatest in two volcano clusters in the Ararat Depression, ncluding the region of the ANPP. Vent density varies by approximately three rders of magn vioitude across the area of interest.
PAGE 24
Figure 7. Plot showing values of h ranging from 1,000 m to 40,000 m versus the probability of an event occurring at the ANPP. The greatest hazard exists when h = 4,000 m. 4000 m to 7.5 x 105 at approximately 25,000 m. This sensitivity of the probability to the smoothing parameter is important to consider because we dont have a robust method for estimating h f volnoes. In oparatively small, volcanism might yet occur in areas that have no record of previous activity. Second, the kernel estimation technique does not attempt to account for additional geologic information that might influence our assessment of the distribution of future volcanoes. Here, we take gravity precisely. Using the above defined parameters, the acceptable range in probability of an event occurring that would disrupt the ANPP and using equation 1 is 1 to 3.4 x 104 for 1800 m < h < 17000 m (Figure 7). Development of the Gravity Model Although a probabilistic volcanic hazard assessment could be based solely on the m ethods outlined above, there are clear shortcomings. First, the distribution ocanoes in the region may incompletely sample the possible distribution of volcather words, because our sample of events is com 17
PAGE 25
Cu m u lati ve Fra c ti o n Gravity distribution near volcanoes Gravidistriovera ty bution ll data and cast gravity anomalies as another probability density function in an effort to he ultimate goal being toavity data in a probability However, before the likelihood function candeveloprelationshipbetween volcano distribution and gravity must be shown to exist. Gravity (mGals) 18 Figure 8. KStest comparison cumulative fraction plot. the gravity distribution near volcanoes shows that volcanoes based on the maximum difference in the two distributions, l Comparison of gravity distribution over the area of study and tend to occur in areas of anomalously low gravity. A KS test, indicates that the tendency for volcanoes to occur in gravity ows i s stat i st i ca ll y s i g nifi ca n t ( > 99% co nfi de n ce). develop a more geologically realistic probability estimate, t modify s with gravity information. However, to use the grmodel, the gravity observations must be transformed into a likelihood function. This likelihood function, like s, is a probability density function be ed a A powerful tool for ch fit s assessing the relationship between two distributions suas volcano distribution and gravity anomalies is the KolmogorovSmirnov (KStest) goodnessoftest (Chakravarti and Roy 1967). In this case, gravity at each specific volcano location i
PAGE 26
compared to gravity data collected across the entire region. Both distributions are plotin cumulative form (Figure 8) and the maximum difference between the two distributions, the KolmogorovSmirnov statistic, (KSstatistic), is measured. For Armenian gravity data, the KSstatistic is 0.202 and the two distributions are different value is ~ 0.0). This indicates volcanoes are clustered in low gravity regions, and thght be modified elihood function. ction and is 04). The first ted (Pis correlation suggests that gravity anomalies may be an indicator of the distribution of future volcanic activity. In other words, we need to consider how s miby the presence or absence of gravity anomalies. Now that a relatiied between volcano location and gravity g the KS test, the next step is to transform the gr into a lik e is currently no standard method for developing funfore somewhat subjective as it is not purely based on statistics and relies to a certain extent on the expertise of those conducting the analysis (Martin et al. 20step in developing the likelihood function is to make some general observations about the relationship between gravity anomalies and volcano distribution. One othere is an inverse relationsh distributions that is, volcano density is highrvati onlyo thhat the p greatshould be correspondingly small. Further observation shows that most volcanoes (about onship has been identif usinTherthere avity data the likelihood bservation is that on, volcanoese conclusion ter than 10 mGal 90%) occur where gravity values are less than 15 mGal. This suggests that the gravity anomalies indicate a threshold in vent distribution. That is, volcanoes are equally likely ip between the two in areas with low gravity anomalies. Expounding on this obse occur in areas where gravity is less than 10 mGal. This leads trobability of volcanoes occurring in areas with gravity values 19
PAGE 27
200 mGal. In he if (15 mGal < G(x,y) < 10 mGal) then W(x,y) = 0.01 assigned to gravity values based on the observetion: to occur anywhere gravity anomalies are less than some threshold value (such as 15 mGal) and not likely to occur at all in areas with gravity values greater than 1this case, we develop a step function to transformation to map gravity values into tlikelihood function using Boolean logic. For example: if (G(x,y) > 10 mGal) then W(x,y) = 0.001 if (G(x,y) < 15 mGal) then W(x,y) = 0.1 where G(x,y) is the measured value of gravity at map location x,y, usually after interpolation onto a grid, and W(x,y) is the weight d distribution of volcanoes with respect to gravity anomalies (Figure 8). Adding additional ifthen statements smoothes the mapped transformation. Using the weights derived from the above process, gravity values are transformed into a likelihood func XYyxWyxGyxL),(),(),( (4) where ),(yxL yxWyxG),(),( is the likelihood function, which integrates to unity over the region of interest XY, and is the set of weighted gravity values. Clearly, the way in which ),(yxL is calculated is a reflection of the geologic interpretation of the data and the experimenters belief about geologic processes governing magma generation, ascent, aeruption. Finally, having developed and confirmed the validity of the gravity model, it is then normalized, recast as a probability density function, and set to the same grid spacing as the probability model. As with the probability model, the final step in developing the gravity model is to prepare a contour map of normalized gravity. The map is made using nd
PAGE 28
21ine the probability PDF (a priori function) of spatial dento a joint probability density function, or a po The product (or s a UTM projection for equal area measurements and the units of gravity are in mGal (Figure 2). Bayesian Model The final step in the Bayesian process is to comb sity and gravity PDF (likelihood function) in steriori PDF for weighted spatial density. intersection) of these two states of information is given by RssdxdyyxLyxyxLyxyx),(),(),(),(),( (5)where x and y represent grid point locations in the region of interest A, is the gravity data, likelihood function. The posteriori PDF is normalized to unientire volcanic field so that cumulative probability does not change, but the shape of the A l egligible or zero in regions where there are no or extremely sparse volcanic events, s (x,y) is the spatial density, which is modified by L(x,y), the gravity PDF, or ty by integrating over the 2D surface distribution will be modified by the gravity PDF (Figure 9) (Martin et al., 2004). imitation with the standard Bayes rule above in equation (5) is the inability toweight the respective PDFs. Further, since it is conditional probability, probabilities will always be n irrespective of how irrefutably geophysical or other information show the potential of new volcanic events to form. We therefore modify equation (5) by combining PDFs through addition (or union) and assign weights:
PAGE 29
22 s syxLayxayx),()1(),(),( RsdxdyyxLayxa),()1(),( (6) ighting function of factor. Relative weight is assigned to the probability and gras es both PDFs equal weight and produces the same result as achieved using equation 4 (Figure 8). Tarise from using this method. First, the aforementioned subjectivity in assigning weights values under this threshold. For values over this threshold gravity values are given i sian where a, is the we vity PDFs subjectively. Giving a large weight to the probability model suggests that future volcanism can be predicted almost exclusively by the distribution of volcanic events. A large weight applied to the gravity model suggests that future volcanism ibetter predicted by gravity anomalies. Setting a = 0.50, giv wo problems to the PDFs and second, a step function must be introduced to determine how much weight should be given to individual gravity measurements. This process is also somewhat subjective but can be given credibility by looking at the graph of volcano fraction versus gravity (Figure 8). The observation that volcanoes cluster in regions under 20 mGals lends credence to weighting the step function heavily in areas with gravity ncreasingly less weight in the step function. When the parameter a is given equal weight in equation 6, we achieve the same result as using the Bayes Theorem equation (5). Now, the new value s (x,y) is used in equation 1, rather than the unmodified value s (x,y). For the ANPP the probability of volcanic disruption is 1 x 10 4 using Bayeinference. Interestingly, this is lower than the probability calculated using the unmodified value s (x,y) of 4 x 10 4
PAGE 30
23 Figure 9. Bayesian output. Map showing results from equation 5 and equation 6 (parameter a is 0.50) where the spatial recurr ence rate PDF and the gravity likelihood function are combined.
PAGE 31
24 Chapter Three Discussion Von Mises (1957), in his treatise on probability and statistics, shows that it is possible to develop an accurate model of spatial density, given enough events. That is, provided the experiment can be performed many times. In geologic hazard assessments we have only one spatial experiment to work with, in this case the observed distribution of volcanoes, and models developed to estimate hazards solely based on these data are inherently uncertain. The risk in using these models is that the distribution of past events may poorly reflect the distribution of potential future events. Any additional information that sheds light on the spatial distribution is therefore worth analyzing. A model based solely on spatial distribution of past events cluster probability around known events and does not predict a future event at positions far from the cluster (Figure 10a). However, a geologic model, such as one based on gravity anomalies, allows probability of future events to be modified to account for our understanding of the geoloot formed in the Quaternary where gravity values are higher than 10 mGal, but do occur through regions with gravity values below 15 mGal. Taking this into account modifies the probability model (Figure 10 b,c). By combining the statistical model with gravity gic setting of volcanism. For example, we know that volcanoes in Armenia have n
PAGE 32
25 x x yy Spatial Density Model: 22121),(hdsieyx 2Nh Gravity Model: XYyxWyxGyxWyxG),(),(),(),( Bayesian Model: yxL),( s RssdxdyyxLayxayxLayxa),()1(),(),()1(),( x y yx),( Figure 10 a,b,c. Stylized diagram showing (a) spatial density model; (b) gravity model; anto a high probability of a future volcanic event, yellow areas to moderate probability of an evassociated wofo occur at location x but would be expected at location y. In the gravity model (b) both x and y fin either loca available for assessing the hazard at both points. d (c) Bayesian model that combines the spatial and gravity models. Red areas correspond ent, and blue areas to regions of low probability of an event. Notice the probability ith location x and location y. In the spatial density model (a) x falls in a region low probability and y a region of high probability, i.e. an event would not be expected tall in a regions of low gravity (or high probability) and an event would be equally likely ation. In the Bayesian model (c) an event is still likely at position y but there is moderate probability of an event at location x. This model reflects the most information
PAGE 33
data through Bayesian inference we arrive at a different, and, if our geologic models are good, an improved volcanic hazard forecast. ear limitations to a purely statistical approach to hazard assessment. One being that this approach takes into account only known features. For instance, in Armenia volcanoes form in association with deep pullapart basins that are a result of transtension and are associated with high rates of sedimentation. The majority of volcanoes in Armenia are monogenetic and hence do not have a great deal of topographic relief. It is conceivable that some Quaternary vents are buried beneath sediment in these basins, and are not accounted for in the estimate of spatial density, but should be. Furthermore, the distribution of volcanoes in the region may incompletely samistribution of volcanoes, i.e. because our sample of events is statistically small, volcanism could occur in areas that have no previous record of activity. Conversely, abrupt changes in geologic structure are common, especially at active plate margins, such as in Armenia. Estimates of spatial density fromhe nt distribution may inadvertently indicate that areas of very different geology are equally likely to host future volcanic activity. supports the idea that In this context, there are several advantages to the Bayesian approach. First, this approach steers us toward a geological basis for making probabilistic hazard assessmentsIn this case we use gravity anomalies as the geologic base from which we modify our probability model. We use gravity because (1) other studies have shown that distributed, monogenetic volcanism often correlates with gravity anomalies (e.g., Connor et al, 2000);(2) the tectonic setting of Armenia and resulting geologic features poi t There are several cl ple the possible d 26
PAGE 34
27y. However, gravity is not unique other pertineer o n this correlation should occur in Armenia; and (3) the KS test statistically shows there is correlation between volcanic events and gravit nt geologic information could be used either alone or in combination with othdata to enhance probabilistic hazard assessments. These include seismic tomographs (Martin et al., 2004), geochemical data (Condit and Connor, 1996), fault and other structural data, and magnetic data. Second, this approach raises the issue that there is no standard method fordevelopment of the likelihood function. Transforming geophysical observations intPDFs and deciding how best to combine them is a subjective process that relies on the experimenters interpretation of geological data and on general observations about the relationship between datasets. Third, uncertainty in probability estimates may actually increase using the Bayesian approach. For instance, the range of probabilities may increase because of uncertainty in how to create the likelihood function, and in how much weight to assigthe likelihood function. Since the posteriori function for ),(yxs is a PDF, increasing probability in regions of low gravity w ith no volcanoes decreases probability in the immed e iate vicinity of volcano clusters. The amount of change depends on the weightingof geologic models. Rather than being a negative outcome, increase in uncertainty more accurately reflects our understanding of the geology. In the case of the ANPP, thprobability of volcanic disruption based solely on kernel estimation techniques is 3.4 x 10 4 for volcanic events impacting the site (A=50 km 2 t = 100 yr, h=3000 m). When the kernel estimation technique is modified with gravity using Bayesian inference, the
PAGE 35
28 es ing ( t =1 yr). These values, whilst apparently low on uman cy probability of volcanism decreases to 1 x 10 4 If other models of volcano formation are used to assess volcanic hazards to the AN PP, the probability may change, depending on both the geologic model and probability model chosen. Perhaps a more useful estimate for this site is simply stating that a range of models based on spatial density of volcano and additional geologic information yield probabilities of volcanic eruptions impact the ANPP site of 14 x 10 6 per year h timescales, do in fact exceed the curr ent International Atomic Energy Agen standard, 1 x 10 7 per year (McBirney and Godoy, 2003), by at least one order of magnitude.
PAGE 36
29 an anism Regionally, Quaternary volcanis m is concentrated in pullapart basins associated with low gravity anomalies, as a result of crustal exte nsion and decompression melting of an enriched mantle. Gaussian kernel estimates of spatial density are greatest in two clusters in the Ararat Depression, including the cluster in proximity to the ANPP. Furthermore, vent density varies by three or ders of magnitude across the region. These models lead to estimates of probability of volcanic disruption to the ANPP of between 1 x 104 and 3 x 104 in 100 years, which exceeds the cu rrent IAEA standard of 1 x 105 in 100 years. This variation indi cates that modification of m odels through the incorporation of additional geologic information, may in crease the range in probability estimates. Ultimately, this paper provides a pathway towards the incorporation of geologic process models in volcanic hazard assessments by allowing these models to be combined with more traditional probabilistic mode ls through Bayesian inference. Chapter 4 Conclusions The approach presented here points to the applicability of assessing pointlike geologic features, in this case volcanoes, through spatial statis tics, specifically, Bayesi statistics. This study shows that geologic stru cture controls the di stribution of volc in Armenia.
PAGE 37
30 Bloomfield, K. 1975, a lateQuaternar volcanic field in central Mexico. Geologische Rundshau v. 64, pp. 476497. Byth, K., and B.D. Ripley. 1980, On sampli ng spatial patterns by distance methods, References y monogenetic Biometrics v. 36, pp. 279284. Chakravarti, Laha and Roy. 1967, Handbook of Methods of Applied Statistics, Volume I relationships in populations. Ecology, v. 35, pp. 445453. Condit, C.D. and C.B. Connor. 1996, Recurrence rates of volcanism in basaltic volcanic Society of America, Bulletin v. 108, pp. 12251241. Connor, C.B. and F.M. Conway. 2000, Basa ltic Volcanic Fields. Chapter in: John Wiley and Sons, pp. 392394. Clark, P.J. and F.C. Evans. 1954, Distance to nearestneighbor as a measure of spatialfields: an example from the Spring erville volcanic field, Az, USA. Geological Encyclopedia of Volcanology Academic Press, 331343. Connor, C.B.; Stamatakos, J.A.; Ferrill, D.A. ; Hill, B.E.; Ofoegbu, G.I.; Conway, M.F.; e basaltic volcanism: Application to a volcanic hazards assessment at Yucca probability of basaltic volcanism: App lication to the Yucca Mountain region. l structural controls on vent distribution, Springerville volcanic field, Arizona, Cressie, Noel. 1991, Statistics for Spatial Data Sagar, B.; and J. Trapp. 2000, Geologic f actors controlling patterns of smallvolum Mountain, Nevada. Journal of Geophysical Research v. 105, n. 1, pp. 417432. Connor, C.B. and B.E. Hill. 1995, Three nonhomogeneous Poisson models for the Journal of Geophysical Research v. 100 (B6), n. 10, pp. 107125. Connor, C.B., Condit, C.D., Crumpler, L.S., a nd J,C, Aubele 1992, Evidence of regiona Journal of Geophysical Research v. 97, pp. 349359. John Wiley and Sons, Inc., New York, 900.
PAGE 38
31 rowe, B.M.; Johnson, M.E.; and Beckman, R. J. 1982; Calculation of the probability of volcanic disruption of a highle vel radioactive waste reposi tory in southern Nevada, USA. Radioactive Waste Management v. 3, pp. 167190. iggle, P.J. 1985, A kernel method for smoothing point process data. Applied Statistics v. 34, pp. 138147. Dixon, T. H. and A. Mao. 2002, REVEL: a m odel for recent plate velocities from space aneberg, W. 2000, Deterministic and probabi listic approaches to geologic hazard o, ChihHsiang. 1991a, Time trend analys is of basaltic volcanism for the Yucca Ho, CH. 1991b, Eruptive probability calcu lation for the Yucca Mountain site, USA: Karakhanian, A. Djrbashian, R., Trifonov, V., Philip, H., Arakelian, S., Avagian, A., a su sionrelated volcanism in Eastern Anatolia, Turkey Geophysical Research Letters v. 30, n. 24, pp.?. Lutz, co. ion. 92ring New York. g; asaki Takahashi. 2004, Modeli ng longterm volcanic hazards through Bayesian inference: an example from the Tohoku arc, Japan. Journal of Geophysical Research v. 109, B10208, pp. 120. C D geodesy. Journal of Geophysical Research v. 107, pp.132. H assessment. Environmental and Engineering Geoscience v. 6, n. 3, pp. 209226 H Mountain site. Journal of Volcanology and Geothermal Research v. 46, pp. 6172. statistical estimation of recurrence rates. Bulletin of Volcanology v. 54, pp. 5056. Baghdassaryan, H., Davtian, V., and Y. Ghoukassyan. 2003, Volcanic hazards in the region of the Armenian Nuclear Power Plant. Journal of Volcanology and Geothermal Research v. 126, pp. 3162. Keskin, M. 2003, Magma generation by sl ab steepening and breakoff beneath bductionaccretion complex: an alternativ e model for colli T.M., and J.T. Gutmann. 1994, An impr oved method of determining alignments of pointlike features and its implications for the Pinacate volcanic field, Mexi Journal of Geophysical Research in press. Magill, C.; Blong, R.; and J. McAnency. 2004, Volcanic loss for the Auckland reg Abstract: Workshop on Statistics in Vol canology, Bristol, England, March 23, 2004. Margulies, T., Lancaster, L .,Eisenberg, N., and L. Abramson. 1992, Probabilistic analysis of magma scenarios for assessing geologic waste repository performance, Rep WA/SAF11. American Society of Mechanical Enginee Martin, Andrew J., Umeda, Koji; Connor, Char les B.; Weller, Jennifer N.; Zhao, Dapen and M
PAGE 39
32 Bulle McBirney, A., Serva, L.,Guerra, M., and C.B. Connor. 2003, Volcanic and seismic Geoth pp. 1130. hazar al Research 126, pp. 19. Ohan Martin, del Pozzo, A.L. 1982, Monogenetic volcani sm in Sierra Chichinautzin, Mexico. tin Volcanologique v. 45, pp. 924. hazards at a proposed nuclear power site in central Java Journal of Volcanology and ermal Research 126 McBirney, A. and A. Godoy. 2003, Notes on the IAEA guidelines for assessing vocanic ds at nuclear facilities Journal of Volcanology and Geotherm issyan, Sh. S., 1985. Structure of the Earth Crust in Armenia Geophysical fields and structure of the Earth crust in the TransCaucasus Nauka, Moscow, pp. 4363 [in Pearce, J.A., Bender, J.F., DeLong, S.E., Kidd, W.S.F., Low, P.G., Guner, Y., Saroglu, Research v. 44, pp. 189229. Philip stron nik fault pp. 205232. c Tectonic Interaction workshop in Iceland. Sheri nic hazards at volcan Popocatepetl, Mexico, EOS, American Geophysical Union 82:185189. Silve Russian] F., Yilmaz, Y., Moorbath, S., and J.G. Mitchell. 1990, Genesis of collision volcanism in Eastern Anatolia, Turkey. Journal of Volcanology and Geothermal H., Ritz, J.F., S. Rebai. 2001, Estimati ng slip rates and recurrence intervals for g earthquakes along an intracontinental fault: example of the PembakSevanSu (Armenia). Tectonophysics v. 343, n. 34, Savov, I.P., Connor, C., DAntonio, M., a nd Y. Jrbashian. 2003, Volcanotectoni interactions in NE Armenia results from petrology, geochemi stry, and Quaternary volcano and fault distribu tion. Presented at the RIDGE 2000Norvolk Volcanodan, M.F., Hubbard, B., Bursik, M.I., Si ebe, C, Abrams, M., Macias, J.L., and H Delgado. 2001, Shortterm potential vol ca rman, B.W. 1986, Density Estimation for Statistics and Data Analysis 175 pp., Chap man and Hall, New York. aste Regul. Anal., 61 625, San Antonio, Texas. Tamu the th and Stamatakos,, J.A. and D.A. Ferrill. 1996, T ectonic processes in the central Basin and Range region, in NRC HighLevel Radioac tive Waste Research at CNWRA. JulyDecember 1995, CNWRA 9502S, edited by B. Sagar, Cent. for Nucl. W ra, Y.; Tatsumi, Y.; Zhao, D., Kido, Y.; and H. Shukuno. 2002, Hot fingers in mantle wedge: new insights into magma genesis in subduction zones. Ear Planetary Science Letters, v. 197, pp. 105116.
PAGE 40
33 Tsuboi, Chuji. Gravity George Allen and Unwin, 1979. London, UK. Von Mises, R. 1957, Probability, statistics, and truth Dover Publications, Inc. NY. al Winkler, K., and A. Nur. 1979, Pore fluids and seismic atte nuation in rocks. Geophysic Research Letters v. 6, pp. 14.
PAGE 41
34 Appendices
PAGE 42
35Appendix A: Armenia Vo lcano Location Data 43.56 41.09 43.59 41.08 43.68 41.13 43.67 41.09 43.67 41.04 43.73 41.04 43.75 41.11 43.82 41.05 43.86 41.03 43.95 41.17 43.96 41.15 43.94 41.13 43.97 41.13 43.96 41.11 43.92 41.10 43.93 41.08 43.95 41.09 43.95 41.08 43.95 41.07 44.00 41.07 43.99 41.10 44.02 41.08 44.00 41.13 43.82 40.57 43.92 40.60 44.01 40.62 44.05 40.59 44.04 40.59 44.05 40.60 44.11 40.60 44.11 40.60 44.12 40.58 44.13 40.57 44.09 40.57 44.09 40.57 44.10 40.55 44.13 40.55 44.09 40.54 44.14 40.49 44.15 40.49 44.13 40.49 44.10 40.47 44.19 40.51 44.16 40.41 44.16 40.40 44.07 40.38 44.08 40.39 44.04 40.37 44.00 40.38 43.96 40.37 43.94 40.47 43.91 40.43 43.90 40.44 43.90 40.44 43.79 40.38 43.79 40.36 43.83 40.35 43.62 40.45 43.59 40.46 43.58 40.48 43.56 40.48 43.64 40.51 43.66 40.51 43.66 40.52 43.64 40.53 44.10 40.69 44.29 40.44 44.30 40.44 44.41 40.37 44.44 40.37 44.45 40.41 44.45 40.39 44.32 40.26 44.32 40 44.34 44.41 40 44.41 40.28 44.41 40.29 44.43 40.29 44.46 40.29 44.45 40.29 44.55 40.21 44.55 40.21 44.55 40.22 44.56 40.22 44.56 40.20 44.57 40.23 44.60 40.23 44.59 40.26 44.12 40.24 44.12 40.23 44.12 40.23 44.12 40.22 44.12 40.22 44.13 40.22 44.13 40.22 44.15 40.23 44.16 40.23 44.15 40.23 44.15 40.23 44.15 40.22 44.16 40.23 44.16 40.21 44.17 40.21 44.17 40.22 44.17 40.22 44.17 40.22 44.18 40.22 44.18 40.22 44.18 40.22 44.19 40.22 44.19 40.22 44.19 40.22 44.19 40.22 44.19 40.21 44.15 40.19 44.15 40.19 44.16 40.19 44.16 40.20 44.15 40.19 44.12 40.18 44.13 40.18 44.13 40.18 44.13 40.19 44.12 40.18 44.13 40.18 44.13 40.19 44.79 40.43 44.87 40.48 44.88 40.49 44.91 40.50 44.68 40.38 44.69 40.37 44.70 40.36 44.69 40.37 43.84 40.26 43.85 40.26 43.85 40.27 43.85 40.26 43.85 40.26 43.84 40.26 43.84 40.27 43.85 40.27 43.79 40.21 43.79 40.22 43.79 40.22 43.79 40.22 43.79 40.22 43.80 40.22 43.76 40.20 43.77 40.20 43.77 40.21 43.77 40.21 .26 40.25 .29
PAGE 43
36 43.78 40.21 43.78 40.20 43.78 40.20 43.79 40.20 43.79 40.20 43.79 40.20 43.78 40.21 43.80 40.20 43.79 40.20 43.77 40.18 43.77 40.17 43.78 40.17 43.79 40.17 43.78 40.18 43.78 40.10 43.78 40.10 43.78 40.11 43.79 40.10 43.79 40.10 43.80 40.10 43.80 40.12 43.80 40.13 43.79 40.14 43.80 40.13 43.82 40.12 43.84 40.11 43.82 40.12 43.83 40.12 43.93 40.12 43.92 40.13 43.92 40.13 43.92 40.13 44.18 40.13 43.95 40.11 43.96 40.10 43.99 40.09 44.03 40.08 44.67 40.35 44.68 40.34 44.68 40.35 44.68 40.35 44.72 40.32 44.73 40.31 44.73 40.30 44.72 40.30 44.72 40.33 44.89 40.43 44.90 40.43 44.89 40.42 44.86 40.42 44.87 40.41 44.87 40.42 44.87 40.42 44.88 40.42 44.88 40.41 44.89 40.40 44.89 40.40 44.90 40.40 44.91 40.41 44.90 40.40 44.91 40.39 44.91 40.38 44.92 40.39 44.91 40.36 44.92 40.34 44.94 40.36 44.94 40.35 44.95 40.33 44.94 40.32 44.90 40.31 44.92 40.31 44.92 40.29 44.92 40.31 44.92 40.31 44.93 40.30 44.93 40.30 44.94 40.29 44.95 40.30 44.88 40.26 44.90 40.27 44.92 40.28 44.94 40.28 44.93 40.25 44.94 40.26 44.95 40.25 44.93 40.23 44.95 40.23 44.94 40.22 44.98 40.29 45.00 40.30 45.00 40.29 45.01 40.29 44.92 40.21 44.92 40.21 44.91 40.20 44.97 40.20 44.96 40.17 44.95 40.16 44.96 40.16 44.96 40.16 44.96 40.16 45.02 40.12 45.01 40.11 45.05 40.12 45.07 40.10 45.08 40.09 45.05 40.04 45.02 40.17 45.03 40.17 45.04 40.15 45.05 40.16 44.96 40.44 44.96 40.40 44.97 40.41 44.99 40.42 44.99 40.44 45.00 40.42 45.00 40.43 45.00 40.43 45.00 40.43 45.00 40.44 45.01 40.42 45.00 40.45 45.01 40.41 45.03 40.42 45.04 40.39 45.02 40.39 45.04 40.38 45.01 40.37 45.06 40.27 45.07 40.26 45.07 40.25 45.22 40.03 45.22 40.07 45.15 40.13 45.16 40.12 45.16 40.11 45.17 40.11 45.16 40.10 45.14 40.09 45.18 40.21 45.18 40.20 45.14 40.38 45.15 40.38 45.21 40.38 45.22 40.39 45.00 40.47 45.38 40.12 45.41 40.10 45.42 40.09 Appendix A (Continued)
PAGE 44
Appendix A (Continued) 45.42 40.08 45.42 40.06 45.46 40.06 45.46 40.06 45.46 40.06 45.51 40.08 45.51 40.09 45.48 40.11 45.57 40.06 45.64 40.05 45.58 40.10 45.47 40.16 45.50 40.13 45.56 40.12 45.59 40.11 45.64 40.13 45.63 40.14 45.63 40.06 45.61 40.03 45.61 40.01 45.72 40.05 45.73 40.06 45.75 40.03 45.74 40.02 45.77 40.03 45.77 40.03 45.79 40.08 45.81 40.07 45.78 40.05 45.81 40.05 45.28 40.02 45.28 40.00 45.30 40.01 45.30 40.01 45.31 40.00 45.39 40.00 45.47 40.00 45.13 40.00 45.32 39.97 45.36 39.98 45.39 40.00 45.40 39.99 45.40 39.99 45.41 39.99 45.47 39.96 45.39 39.93 45.50 39.93 45.52 39.97 45.59 40.00 45.59 39.99 45.61 39.99 45.61 39.97 45.58 39.96 45.56 39.93 45.58 39.91 45.60 39.94 45.61 39.92 45.61 39.93 45.68 39.96 45.70 39.96 45.71 39.96 45.73 39.87 45.49 39.80 45.59 39.89 45.59 39.88 45.61 39.88 45.59 39.82 45.59 39.82 45.63 39.84 45.63 39.89 45.69 39.85 45.69 39.83 45.70 39.81 45.70 39.78 45.79 39.86 45.78 39.75 45.69 39.76 45.70 39.76 45.71 39.70 45.71 39.69 45.71 39.69 45.82 39.86 45.82 39.82 45.83 39.79 45.84 39.79 45.87 39.80 45.86 39.79 45.83 39.78 45.86 39.77 45.86 39.77 45.86 39.76 45.85 39.76 45.84 39.76 45.83 39.76 45.84 39.74 45.88 39.73 45.89 39.74 45.90 39.74 45.87 39.70 45.88 39.69 45.91 39.72 45.89 39.71 45.90 39.71 45.91 39.70 45.91 39.69 45.92 39.67 45.94 39.67 45.94 39.67 45.95 39.65 45.95 39.63 45.96 39.62 45.97 39.59 45.93 39.59 46.14 39.53 46.09 39.54 46.09 39.56 46.06 39.59 46.07 39.59 45.95 39.78 45.94 39.77 45.96 39.76 45.95 39.75 45.97 39.75 45.97 39.73 45.96 39.72 45.99 39.77 45.99 39.77 45.99 39.76 46.00 39.76 46.01 39.77 46.00 39.76 45.99 39.74 45.99 39.76 45.99 39.75 45.99 39.75 46.01 39.76 46.01 39.75 46.01 39.75 46.02 39.75 46.03 39.74 46.01 39.73 46.01 39.72 46.01 39.72 46.01 39.71 46.02 39.71 46.02 39.72 46.02 39.71 46.02 39.71 46.03 39.71 46.03 39.70 37
PAGE 45
38A ppendix A (Continued) 46.03 39.70 46.03 39.69 46.03 39.70 46.04 39.69 46.06 39.71 46.05 39.69 46.05 39.70 46.05 39.70 46.05 39.70 46.06 39.70 46.06 39.69 46.06 39.69 46.06 39.69 46.05 39.67 46.07 39.68 46.08 39.67 46.10 39.68 46.07 39.66 46.02 39.66 45.99 39.66 46.05 39.66 46.05 39.65 46.04 39.64 46.03 39.62 46.06 39.64 46.07 39.64 46.08 39.64 46.07 39.63 46.07 39.63 46.09 39.64 46.09 39.64 46.09 39.65 46.10 39.66 46.11 39.65 46.12 39.66 46.12 39.65 46.10 39.64 46.10 39.60 46.13 39.66 46.12 39.63 46.12 39.62 46.13 39.62 46.11 39.63 46.12 39.63 46.13 39.61 46.14 39.61 46.14 39.60 46.14 39.62 46.14 39.62 46.14 39.62 46.22 39.58 46.24 39.56 46.24 39.58 46.25 39.57 46.24 39.57 46.17 39.54 46.21 39.53 46.20 39.53 46.22 39.54 46.23 39.53 46.23 39.53 46.24 39.53 46.24 39.52 46.21 39.52 46.22 39.52 46.22 39.51 46.25 39.51 46.19 39.50 46.20 39.49 46.22 39.49 46.25 39.49 46.26 39.49 46.23 39.48 46.27 39.49 46.21 39.47 46.26 39.45 46.26 39.44 46.25 39.46 46.25 39.47 46.26 39.47 46.27 39.47 46.28 39.48 46.31 39.46 46.34 39.47 46.31 39.44 46.28 39.43 46.29 39.43 46.28 39.42 46.33 39.40 46.48 39.32 46.46 39.28 46.40 39.26 46.49 39.26 46.50 39.28 46.50 39.28 44.20 40.38 46.16 39.57 46.20 39.58 45.84 40.06 44.95 40.22 44.68 40.39
PAGE 46
39 Codes is project were wL computing lat to run the Bay follows: spacing = to onlat_long_5000.xyple, the Gaussian kernel function code t grid of coordine e spacing like 5.pl converts the "normal" probaesults to pow.oupl redone_volc_gauss_3000.dat > pow.out" ne_volc_gauss_3000.dat > pow.out ethod is all abo he map has to integrate to 1 (unity) um_to_1.pl to fig_it_sum_to_1.pl _to_1.pl pow.oormalize the grility values to inte map (grid) regecho "perl normalize_data.pl pow.out > normalized_pow.out" perl normalize_data.pl pow.out > normalized_pow.out #just check one more time #remember this just prints the value of the summation echo "perl does_it_sum_to_1.pl normalized_pow.out" perl does_it_sum_to_1.pl normalized_pow.out #develop the weighting function. Use the input gravity data here #instead of small_grav001.dat echo "perl grav_wt_function.pl detrended_grav_Xyz.dat > normalized_grav.dat" Appendix B: Computer All codes for th ritten in the PER nguage. The shell scrip esian code is as #!/usr/bin/csh #make your grid e kilometer! #the file small_ z is output #from, for exam # It is the outpu ates and might b # in UTM at som 000 m # the script pow log of #probabilities to bility numbers #then feed the r t echo "perl pow. perl pow.pl redo #the Bayesian m ut "normalizing" #so this means t #run "does_it_s ure this out echo "perl does pow.out" perl does_it_sum ut #go ahead and n d file of probabi egrate to #unity across th ion
PAGE 47
40 perl does_it_sum_to_1.pl normalized_grav.dat echo "perl two_files.pl normalize d_grav.dat normalized_pow.out > v_x_pow.out" erl two_files.pl normalized_gr av.dat normalized_pow.out > .out" .out i zed_grav_x_pow.out > normalize_product.out" ow.out > normalize_product.out cho "perl does_it_sum_to_1.pl normalize_product.out" cho "take log for plotting..." .out > log_normalize_product.out bayes.gmt ipts that make up the shell script above: a, $b, $data) = split; !/usr/bin/perl ARGV[0] means read from the first file after .pl listed on command line, @data creates Appendix B (Continued) perl grav_wt_function.pl gra v_Xyz.dat > normalized_grav.dat normalized_gra p normalized_grav_x_pow.out echo "perl does_it_sum_to_1.pl normalized_grav_x_pow perl does_it_sum_to_1.pl normalized_grav_x_pow echo "perl normalize_product.pl normal perl normalize_product.pl normalized_gr av_x_p e perl does_it_sum_to_1.pl normalize_product.out e perl log_normalize_product.pl normalize_ product echo "./bayes.gmt" ./ The following are the individual scr Script: pow.pl #!/usr/bin/perl while (<>) { ($ $newdata = 10.0**$data; print "$a, $b, $newdata\n"; } Script: does_it_sum_to_one.pl # open (INPUT, $ARGV[0])die "cannont read file!/n"; # an array and assigns variables to each column
PAGE 48
41 hile () { $line = $_; data = split(" ", $line); $y=$data[1]; $sum_z=$sum_z+$z; cript: normalize.pl initializes variables ARGV[0])die "cannont read file!/n"; ARGV[0] means read from the first file after .pl listed on command line, @data creates ssigns variables to each column hile () { @data = split(" ", $line); a[0]; $y=$data[1]; Appendix B (Continued) w $N=$N+1; @ $x=$data[0]; $z=$data[2]; } print "$sum_z\n"; S #!/usr/bin/perl $sum_z=0; $N=0; # open (INPUT, $ # an array and a w $line = $_; $N=$N+1; $x=$dat $z=$data[2]; $sum_z=$sum_z+$z;
PAGE 49
42 #print "$sum_z\n"; )die "cannot read file!/n"; >) =~m/\d/) @data= split (" ", $line); $x=$data[0]; ta[1]; z=$data[2]; $z/$sum_z; print "$x $y $w\n"; rint "[$line]\n"; malize_product.pl Appendix B (Continued) } close (INPUT); open (INPUT, $ARGV[0] while () { ($a, $b, $data) = split;
PAGE 50
43 ewdata = log($data)/log(10.0); ta\n"; cript: grav_wt_function.pl sum_z=0; $N=0; initializes variables ARGV[0])die "cannont read file!/n"; first file after .pl listed on command line, @data creates gns variables to each column INPU ", $line); {$new_z = 0.01}; 2){$new_z = 0.1}; ew_z = 5.0}; 0){$new_z = 20.0}; if ($z < 22){$new_z = 30.5}; if ($z < 24){$new_z = 35.0}; ; 28){$new_z = 45.0}; if ($z < 30){$new_z = 50.0}; $new_z; Appendix B (Continued) $n print "$a $b $newda } S #!/usr/bin/perl $ # open (INPUT, $ #ARGV[0] means read from the an array and assi while () { $line = $_; $N=$N+1; @data = split(" $x=$data[0]; $y=$data[1]; $z=$data[2]; if ($z > 10) if ($z < 1 if ($z < 16){$new_z = 1.0}; if ($z < 18){$n if ($z < 2 if ($z < 26){$new_z = 40.0} if ($z < $sum_new_z +=
PAGE 51
44 Appendix B (Continued) rint "$sum_z\n"; lose (INPUT); pen (INPUT, $ARGV[0])die "cannot read file!/n"; if ($line=~m/\d/) x=$data[0]; data[1]; ]; 01}; 12){$new_z = 0.1}; = 1.0}; = 5.0}; = 20.0}; 22){$new_z = 30.5}; 24){$new_z = 35.0}; _z; } #p c o while () { $line=$_; { @data= split (" ", $line); $ $y=$ $z=$data[2 if ($z > 10){$new_z = 0. if ($z < if ($z < 16){$new_z if ($z < 18){$new_z if ($z < 20){$new_z if ($z < if ($z < if ($z < 26){$new_z = 40.0}; if ($z < 28){$new_z = 45.0}; if ($z < 30){$new_z = 50.0}; $w=$new_z/$sum_new print "$x $y $w\n"; } else { print "[$line]\n";
PAGE 52
45Appendix B (Continued) } !/usr/bin/perl pen(FILE1, $ARGV[0])  die "Cannot open <$ARGV[0]> fo r input: [$@]"; pen(FILE2, $ARGV[1])  die "Cannot open <$ARGV[1]> fo r input: [$@]"; ed) =0.50; ) { >; data1)+((1$a)*$data2); T code for mapping the output !/bin/ mtset R_FONT 5 5 _SIZE 16 e interpolation w ith a minimum curvature algorithm d spacing is 2 units in each direction ica tes the output file orth bounds of the map } Script: two_files.pl # o o Appendix B (Continu $a while ( ($a1, $b1, $data1) = split; $line =
PAGE 53
46 Appendix B (Continued) log_normalize_product.out Gnormalized_prod.grd 80000/662000/4299000/4573000 olor table (to color shade contours) e basic color table T specifies the range and interval #psmask clips or masks area of R gives range of data tickmark info 0000 Jx0.000022 B25000a50000/WSne K V P Y1.5 > 0.ps (color map) JX6.0i is the scale (must match the following) lat_long_5000.cpt is the colo r scale created with makecpt P portait mode (must match the following) (dots per inch) of the color shading t to be appe nded to cn.ps in the following dimage normalized_prod.grd Jx0.000022 Cnorma lized_grav.cpt P E100 K V O grdcontour draws the contours from the grid ap will be 6 inches wide C250 means there is a 250 nT contour interval are annotated every 500 nT frame, 25 m tic k with 50 m label, add 5 m ticks, label on h side only ge rdcontour normalized_prod.grd Jx0.000022 C.5 A2 L8 W0.25p P O K V >> surface I1000 V R3 # makecpt creates a c # Cseis specifies th # makecpt Cseis T8/3.5/. 25 V I > normalized_grav.cpt no data on a map # #B gives #I grid spacing psmask grav_utm.dat R I1 bayesian_product_5 # grdimage plots the image # # C # # E is the dpi # K more postscrip gr >> bayesian_product_50.ps # # JX6.0i means the m # # A500 means the contours # B25a50f5/WSne draw the # west and sout # W0.25p set line width # P draw in portrait mode # O overlay contours on the ima g bayesian_product_50.ps psmask C O K P V >> bayesian_product_50.ps
PAGE 54
47Appendix B (Continued) G0 makes the triangles solid K V >> sities plotted sscale D2.65/.5/5/0.2h Cnor malized_grav.cpt B1/:"log ( volc/km@+2@+)": O P I tm.dat Jx0.000022 R W2.0 O M V K P sxy lake_sevan_utm.dat Jx0.000022 R W1.25 P G255 V O K >bayesian_product_50.ps >> bayesian_product_50.ps oduct_50.ps OF sxy R Jx0.000022 O P G0/0/ 0 Sa0.15i W0.5 V <> ate the recurrence rate which is then mapped in Purpose: this script reads vol cano location data from a file and calculates the spatially onhomgeneous recurrence rate # psxy plots volcanic vents as solid black triangles # S specifies the size and shape # psxy volcano_type.dat R Jx0.000022 S t0.05i G0 O P bayesian_product_50.ps #add a color scale to show the range of lat_long_5000 in ten p V K >> bayesian_product_50.ps psxy armfaults_u >>bayesian_product_50.ps p > pstext R Jx0.000022 G0 O P V P K <> bayesian_pr 443000 4430000 12 0 24 BL Yerevan E p bayesian_product_50.ps 440000 4440000 EOF The following script is used to calcul GMT. #!/usr/bin/perl #this is a perl script by Jenn Weller #created on Feb. 24, 2003 # n #using a guassian kernal function.
PAGE 55
48 alue of the point density estimate is GMT. put file to contain two columns of numbers giving location of olcanoes. volcanoes=(); $line = $_; @data = split(" ", $line); a[0]; $volcanoes[1][$n]=$data[1]; #print "$volcanoes[0][$n] $volcanoes[1][$n]\n"; { 06; $y>4.29899e+06; $y=1000){ grid x<662001; $x+=1000){ points on grid Appendix B (Continued) #The location of the point density estimate a nd the v output (x y z); #this output can be contoured in #This code requires the in v #initialize variables used in calculations $n=0; $x=0; $y=0; $h=3000; @ open (INPUT, $ARGV[0])die "cannot access file!/n"; while () { if ($line=~m/\d/) { $volcanoes[0][$n]=$dat $n=$n+1; } else print "[$line]\n"; } } for ($y=4.573e+ #steps through a ll y points on for ($x=380000; $ #steps through a ll x
PAGE 56
49 $sum_k=0; or ($ct=0; $ct<$n; $ct++){ data points #calculates the distance from point x,y to point $ct $kernel=1/(2*3.14159)* exp(0.5*$dist**2/$h**2); $sum_k=$sum_k + $kernel; #calculates the kernel and the sum respectively } $lambda=1/($n*$h**2)*$sum_k; #lambda is calculated Appendix B (Continued) f #steps through all $dist=sqrt(($volcanoes[0][$ ct]$x)**2 + ($volcanoes[1][$ct]$y)**2); #print "$dist\n"; #print "$kernel\n";
PAGE 57
50 Appendix B (Continued) 0001){ lambd #also multiple by 1e6 to report answer in terms of volcanoes/km^2 rather an m if ($lambda>=0.0000000000 $a2=log($lambda*1e6)/log(10); #the log of lambda is calculated for ease of use in contouring #note that "log" in perl is ln so divid by log(10) th^2 } else{ $lambda2=13; } print "$x $y $lambda2\n"; } } close (INPUT);
PAGE 58
51 Map of spatial density h=1300 m. Appendix C: Additional Maps
PAGE 59
Map of spatial denstity h=6000 Appendix C (Continued) 52
PAGE 60
Appendix C (Continued) 53 Map of spatial denstity h=8000
