University of Southern California Dissertations and Theses

Numerical modeling of long period events at the Piton de la Fournaise by the use of object-oriented scientific computing

(USC Thesis Other)

Numerical modeling of long period events at the Piton de la Fournaise by the use of object-oriented scientific computing

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy subm itted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these w ill be noted. Also, if unauthorized copyright material had to be removed, a note w ill indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6" x 9" black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UM I directly to order. ProQuest Information and Learning 300 North Zeeb Road. Ann Arbor, M l 48106-1346 USA 800-521-0600 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. NUMERICAL MODELING OF LONG PERIOD EVENTS AT THE PITON DE LA FOURNAISE BY THE USE OF OBJECT-ORIENTED SCIENTIFIC COMPUTING by Jinbo Chen A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Geological Sciences) December 1999 Copyright 1999 Jinbo Chen Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UM I Number. 3054858 _ ___ ( 8 ) UMI UM I Microform 3054858 Copyright 2002 by ProQuest Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. ProQuest Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor. M l 48106-1346 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES. CALIFORNIA 90007 TTiis dissertation, written by .......................................... under the direction of ft.i£ Dissertation Committee, and approved by all its members, has been presented to and accepted by The Graduate School, in partial fulfillment of re- quirements for the degree of DOCTOR OF PHILOSOPHY Date August 23, 1999 DISSERTATION COMMITTEE Reproduced with permission of the copyright owner. Further reproduction prohibited without permission D edication This thesis is dedicated to my mom Wang Zong-ling, who passed away in December 6, 1998, and m y baby expected to arrive in August, 1999. ii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A cknowledgm ents First, I want to thank my wife for her tremendous support through these years, which are beyond any words I can use. I want to thank Kei Aki for a lot of things. Most of all, he is really the first one to show me what is science and what a real scientist is. I still vividly remember that he told me again and again “Alway compare with observation”. I feel my life has been lifted in many ways by a lot of inspiring persons. Kei is one of the inspiring persons. In the basement of Geological Sciences Building of USC. there were happy people down there, Dr. David Adams, Dr. Michael Forrest, Dr. Chen Xiaofei. Beltas Periklis, Dr. Ouyang Hongping, Dr. Wang Jin, Dr. Qu Jiang, Dr. Li Chengpin, Dr. Yonggong Li and me. V V e had been living there happily for many years until we faced inevitable fate - finding a job. Thank you guys for the time, support and companionship, especially Dr. Yong-Gong Li as a true friend. I truly respect people who educate me and bring good out of me , for this purpose, I want to thank Ms. Jin (my elementary teacher), Mr. Gao (my undergraduate math instructor at University of Science and Technology of China, showed me true passion on m ath), Dr. Huang Peihua (deceased in 1999), Dr. Xu Guoming, Dr. Charles iii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Sammis, Dr. Steve Lund and Dr. Leon Teng and Dr. Wellford. W ithout these educators, my life would be totally different. God bless teachers. Finally I want to thank my friends at Caltech - Dr. Ding Xiaoming, Dr. Wang Yichun and friends at Herman Chinese Church for caring and supporting me through hardship. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C ontents D edication ii Acknowledgm ents iii List O f Tables vii List O f Figures viii A bstract xi 1 Introduction 1 1.1 Geological and Geophysical Observation ................................................. 1 1.1.1 Geological Observation of Exposed Volcanoes ........................... 1 1.1.2 Surface Displacement M easurem ent.............................................. 3 1.1.3 Tom ography....................................................................................... 5 1.1.4 Volcanic Tremor and Long Period E v en ts..................................... 7 1.2 Volcano Source M odeling.............................................................................. 10 1.3 O bjectives....................................................................................................... 11 2 Long Period Events at Piton de la Fournaise 13 2.1 Introduction.................................................................................................... 13 2.2 Identification of long-period e v e n t s ........................................................... 19 2.3 Summary of observations on LP events at Piton de la Fournaise . . . 25 2.4 Finite Element M odel................................... 31 3 Finite Elem ent M odeling Im plem entation 34 3.1 Introduction to Finite Element M ethod..................................................... 34 3.1.1 Basic Concepts in Finite Element M e th o d ................................. 34 3.1.2 Finite Element vs Finite Difference Method .................................. 44 3.2 Object-oriented scientific com puting........................................................... 46 3.2.1 Object-oriented technology.............................................................. 46 3.2.2 Object-oriented language C + + vs Procedure-oriented languages 49 3.2.3 Objective and m otivations.............................................................. 50 v Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.3 Introduction to D iffp ack ............................................................................... 52 3.4 Long Period Event Simulation System ...................................................... 55 3.4.1 ModelPack.............................................................................................56 3.4.2 Generic Elastic Wave S im ulator.......................................................60 3.4.3 Result R e p o rt...................................................................................... 61 3.4.4 Documentation and future development..........................................61 4 Num erical Sim ulation R esults 63 4.1 Homogeneous half-space Line Pressure S o u rc e ........................................... 63 4.2 Vertically elongated magma c h a m b e r..........................................................69 4.3 Horizontally elongated magma cham ber.................................................. 73 4.3.1 Effect of the w i d t h .............................................................................73 4.3.2 Effect of the d e p t h .............................................................................80 4.3.3 Effect of the height(thickness)..........................................................86 4.4 Refined m o d el...................................................................................................93 4.5 Discussion and Sum m ary................................................................................97 A ppendix A Finite element procedures, Object-oriented concepts and C + + header files 101 A.l Basic Procedures in Finite Element M e th o d ............................................101 A.2 Basic concepts in the Object M odel............................................................108 A.3 C + + System Header files.............................................................................. 115 A.3.1 Elastic W av e........................................................................................ 115 A.3.2 Numerical Integration .....................................................................121 A.3.3 Elastic Material H a n d lin g ...............................................................123 R eference List 138 vi Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. List O f Tables 2.1 Seismic crises at the Piton de la F o u rn aise.............................................. 29 2.2 Coda amplification factor at stations near the su m m it........................... 29 2.3 Table of normalized amplitude for 31 LP events from 1990 - 1995 . . 30 4.1 The normalized vertical surface ground velocity amplitude distribu tion of the homogeneous half-space model with line pressure source at different depths, the bracketed values in bold are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km. 68 4.2 The normalized vertical surface ground velocity amplitudes for the following cases: (I) a vertically elongated magma chamber (2) ho mogeneous half-space, (3) a horizontally elongated magma chamber (this case is studied in the next section). The bracketed values in bold are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km........................................................................ 70 4.3 The normalized vertical surface velocity amplitude distribution of the horizontally elongated magma chamber with fixed height 400m and different width, (*)the normalized amplitude is based on P waves. The bracketed values in bold are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km......................................... 75 4.4 The normalized vertical surface ground velocity amplitude distribu tion for the case of a horizontally elongated magma chamber at dif ferent depth. The bracketed values in bold are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km. 80 4.5 Parameters for the Crosson-Bame‘s model, where Q = Im(n)/Re(fi) 94 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. < M O l Cl Cl Cl Cl Cl Cl Cl Cl List O f Figures 1.1 A LP event at the Piton de la Fournaise of La R eu n io n ....................... 8 .1 Location of the Piton de la Fournaise of La Reunion.............................. 16 .2 Eruption of the Piton de la Fournaise of La Reunion ........................... 16 .3 A relation between the precursor time of the seismic crisis and ele vation of the eruption site. Only those crises preceding the rift-zone eruptions are accompanied by the occurrence of LP e v e n ts ................ 17 .4 Piton de la Fournaise of La Reunion ....................................................... 18 .5 Stations at the Piton de la Fournaise, La R eunion.....................................21 .6 Amplitude distribution of seismic waves in the frequency band from 1 to 3 Hz from a Mahavel rock fall .......................................................... 23 .7 Distribution of median amplitude of LP events normalized to the mean of the three summit stations; obtained from (a) 14 events oc curred in the period from 1990 to April, 1991, (b) 11 events during the intrusion crisis of December, 1991, (c) 6 events after the crisis till April, 1995. (d) shows the distribution of the amplitude of a large volcano-tectonic event during the crisis of December, 1991, band-pass filtered between 1 and 3 Hz. The number attached to each station is normalized to the mean of the three summit stations..................................26 .8 The histogram of annual occurrence of LP events with a given domi nant frequency. The vertical scale is common to all years, exception for 1991, for which it is 4 times greater than the rest even though the year is separated into two p a r t s ................................................................ 27 .9 Finite element model of Long Period e v e n t s ........................................... 32 .10 Uniform mesh used in finite element modeling, which is better than other kinds of m eshes................................................................................... 33 3.1 Long Period Event Simulation S y s te m .........................................................56 3.2 Schematic diagram of constructing a model with complex velocity stru ctu re..............................................................................................................60 4.1 The comparison between the observed arrival times and those for a line source at different d e p th s ........................................................................67 viii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.2 Vertical surface ground velocity for the case of a line source at depth 2.0 k m .................................................................................................................68 4.3 Vertical surface ground velocity for the case of a vertically elongated magma chamber with height 2.0 km and width 0.5 k m ........................... 71 4.4 Amplitude comparison of vertical surface ground velocity for the case of a vertically elongated magma chamber with height 2.0 km and width 0.5 k m ................................................................................................ 72 4.5 Vertical surface ground velocity for the case of a horizontally elon gated magma chamber at depth 2.0 km with height 500m, width 2.0 km and rise time 0.3 sec. This case is used to compare with the case of vertically elongated magma chamber..................................................... 75 4.6 Vertical surface ground velocity for the case of a magma chamber height 400m and width 0.577km with rise time 0.3 se c.......................... 76 4.7 Vertical surface ground velocity for the case of a magma chamber height 400m and width 1.04km with rise time 0.3 s e c ...............................77 4.8 Vertical surface ground velocity for the case of a magma chamber height 400m and width 1.62km with rise time 0.3 s e c ........................... 78 4.9 Vertical surface ground velocity for the case of a magma chamber height 400m and width 2.1km with rise time 0.3 sec .............................. 79 4.10 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 1.0km with use time 0.8 s e c .................................81 4.11 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 2.0km with rise time 0.8 s e c .................................82 4.12 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 2.5km with rise time 0.8 s e c .................................83 4.13 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 3.0km with rise time 0.8 s e c .................................84 4.14 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 4.0km with rise time 0.8 s e c .................................85 4.15 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.1 km) with rise time 0.3 s e c ....................................................................86 4.16 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.2 km) with rise time 0.3 s e c ....................................................................87 4.17 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.6 km) with rise time 0.3 s e c ....................................................................88 4.18 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) with rise time 0.3 s e c ................................................................ 89 4.19 Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.5 km) with rise time 0.3 s e c ....................................................................90 ix Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.20 Spectra of the vertical surface ground velocity for the case of a magma chamber with 500m height with different dominant frequencies of source pressure time-function. From left to right, the dominant fre quencies of the line pressure are 0.5, 1.0, 1.5 2.0 and 3.0Hz....................... 92 4.21 Schematic diagram of the complementary m o d e l...................................... 93 4.22 Amplitude spectrum for the case of inner radius of 1.7m and outer radius of 200m, the dominant frequency is 1.0H z........................................95 4.23 Amplitude spectrum for the case of inner radius of 0.8m and outer radius of 200m, the dominant frequency is 2.0 H z .....................................96 4.24 The relation between the inner radius (m) and the dominant frequency and inferred relation between the depth (km) and the dominant fre quency ..............................................................................................................98 A.l Sparse Pattern of the Stiffness matrix of a 3D m o d e l ........................... 106 x Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A bstract Numerous long period(LP) events with the dominant frequency in the range of 1- 3Hz have been observed during active periods of the Piton de la Fournaise in the past several years with the following characteristics: • their arrivals at the stations located along the rim of the summit caldera are simultaneous and the waveforms are similar. • their amplitudes(vertical component) are the largest at the summit area, decay to 0.4 to 0.5 of the maximum at a distance about 5 km from the summit, using station correction by the coda method. We started with the simplest model of a pressure line-source buried in a homoge neous half-space, and found that the observed arrival times of the first cycle of the event can be explained very well with the source depth of 2 km using seismic ve locities appropriate for the vicinity of the summit area. This depth is close to the focal depths of volcano-tectonic earthquakes preceding an eruption which occurred as a clustered swarm beneath the summit. This simple model, however, overesti mates the amplitudes at the summit stations relative to the stations at 5 km away from the summit. In order to eliminate the above discrepancy, we introduced a low xi Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. velocity region( magma reservoir) to house the pressure source. We found that a ver tically elongated reservoir would increase relative amplitude at the summit, thereby increases the discrepancy with observation. On the other hand, the horizontally elongated body tends to reduce the discrepancy. We thus infer that the low velocity magma chamber responsible for the generation of LP events under the Piton de la Fournaise is horizontally elongated. For horizontally elongated low velocity magma chamber, we investigated how the dimension and depth of the chamber affect the surface amplitude distribution. We found that for a fixed depth and height( thickness), the wider magma chamber generally results in greater normalized amplitude at 3.0 km or farther. As our simulation results show that, because of the introduction of the low velocity magma chamber, the depth effect for a horizontally elongated magma chamber is quite different from the depth effect in the case of line pressure source in a homogeneous half-space. For a fixed height (thickness) and width, when the depth is below a critical value, the normalized amplitude at epicenter distance 5.0 km increases at a linear rate, but beyond the critical depth, the normalized amplitude at 5.0 km stays constant with increasing depth. While in the case of homogeneous half-space, the normalized amplitude at distance 5.0 km increases linearly with the depth. The combined effect of depth and width is a very important factor in determining the amplitude distribution. By examining the surface synthetic seismogram, we found that amplitude of S wave is much greater than the amplitude of P wave at station at distance about 3 km or farther. From comparison of simulation with observation, we conclude that the LP events source is located at a depth 2.5 km and has width of about 2.5 km. xii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Our simulation results show that the dominant frequency is not related to the thickness of the low velocity magma chamber, but is determined by the dominant fre quency of the line pressure source. At the Piton de la Fournaise, we have observation supporting that the 1Hz LP generators are located above the 2Hz LP generators. To explain the different dominant frequency of the LP events at different depths at the Piton de la Fournaise, we propose a refined model with a small gas chamber unbound within the fluid chamber. Crosson and Bame [18] investigate the characteristics of such a model in a homogeneous space. V V e analyze Crosson-Bame’s model using parameters appropriate for the LP event source at the Piton de la Fournaise. We suggest that the different inner radii of the gas chamber at different depths could explain different dominant frequencies at different depths. In this thesis, we used object-oriented scientific computing frame work to facili tate collaboration in seismological community. We developed our program based on a well-recognized free scientific framework - Diffpack which has been developed at Oslo University in C + + . We benefited greatly by utilizing many modularized objects from Diffpack, which also make our program easy to be used by other researchers . xiu Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C hapter 1 Introduction u All this time, the guard was looking at her, first through a telescope, then through a microscope, and then through an opera glass” -Lewis Carroll, Through the Looking Glass 1.1 G eological and G eophysical O bservation In this chapter, we shall first summarize geological and geophysical observations on active volcanoes relevant to the present thesis work. 1.1.1 Geological Observation of Exposed Volcanoes The geology of eroded volcanoes gives ample illustrations of the complex structures of the volcanoes. A short encyclopedia of such geologic forms and structure is presented by Frankel [34]. The four general forms are sills, dikes, laccolith and plugs. The thickness of a sill can be remarkablely uniform across an exposure of 10 km or more. The range of the thickness is from about a meter to a few hundred meters. A sill, like other intrusive bodies described here, may form by successive injections 1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. of magma. Because the horizontal geometry of a sill does not intersect the surface, a sill is usually considered only as a means of storing magma, perhaps temporarily. The geometry of joints in a parallel set of sills at Katmai, Alaska, shows that sills are also a means of transporting magma horizontally [63]. If horizontal extension of a sill is retarded, so that continued injection of magma causes vertical growth of intrusive body and arching of overlying layers, then the sill evolves into a laccolith. Laccoliths are less common than sills, though they are not rare. More than 600 have been described in the United States [17]. Field observations suggest that sills are the forerunners of laccoliths, which form only after magma has spread far enough laterally to gain leverage to bend the overlying layers upward [48]. On the basis of cooling rate, the time required for a sill or laccolith to form within cold rock is less than a few years [17]. A volcanic plug, also called a “volcanic neck” is a vertical, pipelike body of magma that may be conduit of a former volcanic vent. Such a feature has no exposed or inferred floor. Instead, a plug connects to a lower magma body that was the source of eruptions. A clear example of this field relation is exposed near the Cloudy Pass batholith, Washington [13]. The southeastern end of the batholith plunges beneath older crust pierced by a line of plugs. A dike is a thin intrusive body with a vertical or nearly vertical orientation . Often, dikes are crowded in swarms, either as radial pattern around a central point or along a long, linear system. The central point of a radial pattern may be a complex of plugs or laccoliths, such as at Sunlight volcano, Wyoming [74] and the island of Rhum [10]. Radial dikes on the dormant volcano Mount Rainer, Washington, radiate from the summit that overlie a central plug [33]. A radial dike pattern may result o Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. from the tensions produced when pressure increases within a central conduit, such as a plug. The expansion causes the surrounding rock to split in vertical cracks that extend outward from the conduit [41]. In Iceland, erosion has exposed thousands of dikes aligned parallel to the axis of the mid-ocean ridge system. In a survey of 700 basaltic dikes in eastern Iceland, Gudmundsson [36] determined an average width of 4 m and followed the longest exposed dike for more than 22 km. In Hawaii, basaltic dike widths measured at the eroded Koolau range from 0.05 m to 6.7 m; the average width is 0.7m [90]. 1.1.2 Surface Displacement Measurement The surface displacements near volcanoes have been studied extensively in order to understand the physics of active volcanoes before and after the eruptions. Historically the patterns of the surface displacement were determined mainly by using high-precision land surveying techniques, originally developed to make topo graphic surveys. These techniques measure changes in elevation, distance, or angle among a set of reference points called bench marks. Recently new techniques have been introduced into this field, such as Global Positioning System (GPS). Regardless of the techniques used, the strategy in making field measurements is simple. The relative positions are measured of a network of bench marks scattered across the surface. Some time later, the same points axe measured again to determine changes since the previous measurement. A balance is sought between infrequent measure ments of a large number of bench marks and frequent measurements of a few bench marks. 3 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Surface displacement has been measured at several dozen volcanoes worldwide, though only a half dozen or so have been studied intensely for a few decades. Usually the measurements do not begin at a volcano until a crisis has developed, which is mostly signaled by a shallow earthquake swarm, and so, usually only unknown fraction of the total surface displacement associated with a volcanic event is recorded. The amount and rate of surface displacement can be astounding at an active volcano. A meter or more of displacement records in a day or less is not unusual. Smaller amounts and rates are also common, over an area of a hundred square kilo meters or more. Such a wide ranges of possible surface displacement illustrates the difficulty in making measurements. To make useful measurements, techniques must have (1) a broad dynamic range to ensure that an entire event is recorded. (2) suffi cient sensitivity to record slight and slow changes, and (3)adequate spatial coverage to capture an event. A single technique cannot satisfy all three requirements, and so multiple techniques are often used. Detailed discussions of survey errors and ac curacies and of the installation and stability of benchmarks are given by Ewert and Swanson [27] and Murray [67]. A description of the installation and operation of tiltmeters is given by Toutain et al. [68] The surface displacement measurement on some volcanoes including Sakura- jima,Campi flegrei, Long Valley, and Kilauea showed that the pattern and rate of the surface displacement can often be fitted quite well by a model in which a few simple pressure sources embedded in the elastic material, which might seem to be at odds with internal complexity evident at the eroded volcanoes. The displace ment pattern reveals such characteristics as depth and rate of magma accumulation. 4 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Near the eruptive site, the elastic limit is exceeded, and the pattern of the surface displacement becomes chaotic as the faults and fissures develop. 1.1.3 Tomography The volcanoes usually have low velocity bodies beneath the caldera. The region beneath the Yellowstone caldera has a lower seismic velocity and shear wave shad ows are considered as evidence for a shallow body of partial melt. A 3-D model of seismic velocity beneath Yellowstone shows a large low velocity zone about 100km in diameter that extends from surface to a depth of about 250 km [43]. The velocity reduction is 10% throughout most of the low velocity zone, and the reduction is as much as 30% in two small areas within 10km of the surface beneath northwest and southwest sections of the caldera [59]. Major changes in temperature and composi tion across the caldera boundary, without the presence of magma, could account for the 10% reduction of velocity [59]. But the 30% velocity reduction, does require a partial melt at depths greater than about 2 km, with melt fraction that ranges from 10 to 50% [70]. By analogy with Yellowstone, the 30% velocity reduction beneath the Valles caldera, New Mexico, is also evidence of a body of partial melt [76]. The body is about 17 km across, is centered at depth of 10 to 13 km, and probably have a vertical height of at least 8 km. Elsewhere, the geophysical evidence is not as compelling for the presence of shallow magma, even though continued eruption indicates that shallow magma must be present and shallower magma storage location is more critical to the prediction of volcano eruptions. 5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. An aseismic, low velocity zone beneath the summit area of the basaltic volcano Kilauea, Hawaii, represents a zone of temporary magma storage [87]. Because the velocity reduction is only a few percent, the zone is probably mixture of magma and solid. Clear identification of magma storage is complicated by high-velocity zones within Kilauea, which represent dense intrusions [78]. On the basis of the pattern of movement of Kilauea’s south flank after a magnitude 7.2 earthquake in 1975, Delaney et al., 1990 [20] proposed that a magma storage zone extends from 3 to 9 km beneath the summit region and both of Kilauea’s rift zones. At Kraft, Iceland, another basaltic system, a magma body is inferred from atten uation of shear waves [25]. The zone of strong attenuation is about 2 km wide, lies within 3 km of the surface, and extends no deeper than 7 km. A similar zone may represent magma storage shallower than 3 km beneath the basaltic volcano Plosky Tolbachik, Kamchatka [7]. Beneath Long Valley caldera, California, the velocity reduction is only about 10% [86], and so is not unambiguous indicator of a partial melt. Instead, attenuation of shear waves in two small zones [80] and refraction of seismic waves are better indicators of possible shallow magma storage. Results of seismic reflection surveys are used to infer structures that could be the tops of small magma bodies. Such a body lies about 4 km beneath the floor of Campi Flegrei caldera, Italy, at the base of an intense zone of shallow earthquakes[32]. A 20-km-deep body lies beneath Etna [83]. Results of other reflection surveys reveal a sill-like body, about 8 km wide and less than 0.5 km thick, about 20 km beneath Socorro, New Mexico, along the Rio Grande Rift Zone [1]. Continuous seismic profiles across the several segments of the mid-ocean ridge system reveal the top 6 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. of crustal magma chambers 1-2 km below the seafloor[22]; Collier and Sinha, 1992 [16]]- No evidence of either low-velocity or attenuation zones has been found beneath isolated andesitic or dacitic volcanoes in the Cascade Range, in particular Mount Hood, Oregon, and Mount Shasta and Mount Lassen, California [43]. The lack of evidence means that magma rises from a deep source through a narrow conduits that are too small to be detected by available seismic methods. A small magma body is inferred from seismic studies to lie about 5 km beneath the summit caldera of the mixed basaltic and rhyolitic Medicine Lake volcano in California [26], and Newberry caldera, Oregon [92]. 1.1.4 Volcanic Tremor and Long Period Events During an active period, earthquakes similar to usual tectonic earthquakes occur in and near a volcano. These volcano-tectonic earthquakes are triggered by the increased stress caused by magma injection to magma reservoir. We can determine the locations of these earthquakes easily, because they show clear P and S wave. But these signals are not directly caused by the magma activities. Omori [72] [73] published the first observations of the volcanic trem or recorded on Usu-San and Asamayama volcanoes in Japan. Sassa [81] [82] gave the first very detailed quantitative analysis of the volcanic tremor on Mt. Aso. He discovered tremors with period ranging from 3.5 to 7s, and he suggested that this kind of tremor is related to oscillations of the magma reservoir because of the very long period observed. 7 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. BOR A A \ I \yVvv^vVv y ' U / ITA ' V ^ v / 'yV'VVxA'^' C H R f t NTR NCR V \/V \W ^ // ^ / v A , 1 s Figure 1.1: A LP event at the Piton de la Fournaise of La Reunion 8 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Minakami [64] found small seismic events characterized by indistinct S phase and dominated by surface waves that he called B-type earthquakes. He distinguished B- type events from explosion earthquakes, but noted their waveform similarity and considered the possibility that a B-type event could be an explosion earthquake of small amplitude. Explosion earthquakes are often accompanied by an air phase. On Langila volcano, in New Guinea, Mori et al., [66] observed all kinds of events from air phase to B-type events showing gradual changes in wraveform, starting with an air phase and a B-type alone, passing through a B-type event accompanied with an air phase and a B-type event alone. They suggested that the observed waveform depends on whether the explosion occur in an opened or closed system, a similar interpretation as Minakami [64]. On Sakurajima, Kamo [49] reported that the foci of explosion earthquakes are distributed at depths ranging from 1 to 3 km beneath the crater, and Ishihara [42] mentions that the first motion is always the compression independent of the azimuth with respect to the tremor source. B-type earthquakes are also called long period events. While Minakami [64] defines B-type earthquakes as having an extremely shallow focus, long period] LP) event is a more generic term. LP event has the characteristics of having no clear P and S wave and dominant frequencies much lower than the dominant frequency of volcano-tectonic event of a similar amplitude, and an example of LP events with such characteristics is shown in Fig 1.1. Fehler and Chouet [29] present convincing evidence that the unusual character of these events is primarily due to a source rather than path effect, and they further indicate that the frequency stability suggests that in contrast to volcano- tectonic earthquakes, the source size is relatively constant. Most LP events appear to originate at shallow depths] 1 to 3 km beneath the volcanic edifice) [64]. 9 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Furthermore, Latter [57] proposed that the tremor of volcanoes is just a sequence of randomly occurring long period events at Ruapehu and Ngauruhoe, New Zealand. This is supported by Fehler [28] who found out that the volcanic tremor and LP events on Mount St. Helens share the same spectral shape. Aki and Ferrazzini [5] have been conducting extensive observation and research on LP event at Piton de la Fournaise of La Reunion. These LP events have dominant frequencies in a range of 1 to 3 Hz. This thesis will be concerned with these LP events observed at the Piton de la Fournaise. Chapter 2 gives a detail description of long period events at Piton de la Fournaise of La Reunion. 1.2 Volcano Source M odeling Mogi [65] used a small spherical pressure source buried in a homogeneous half space to explain the measured surface displacements at several volcanoes. The Mogi model successfully reproduces the pattern of vertical displacement at many volcanoes during either uplift or subsidence. The depth of pressure source has been determined at more than a dozen volcanoes. Some shallow sources is located less than 2 km below the surface. Reliable depths deeper than 6 km have been determined only for a large caldera. More than one point source is needed to explained the displacement patterns for some volcanoes. Recently, Okada's [71] formulas for dike intrusions are commonly used for interpreting displacement associated with fissure eruptions. A number of models have been suggested to explain the source of LP events. Shima [85] [84] considered the free oscillation of a spherical and ellipsoidal magma chamber buried in an unbounded elastic space as a possible source for the different 10 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. kinds of tremor observed by Sassa on Mt, Aso. He provided for this model with very complete analysis of seismic radiation. Based on Shima’s results, Kubotera [51] used the fundamental mode of a spherical liquid-filled magma chamber to explain very long period volcanic ground motion. Riuscetti et al. [79] suggested that free oscillations in the cylindrical magma chamber explained the two-peaked spectra that they observed at Mt. Etna. Aki et al. [3] proposed that volcanic tremor may originate by the repeated opening of tensile cracks by fluid injection. Chouet [14] has applied this model to predict the details of the ground motion near a tremor source. The model of resonant oscillations of a fluid-filled cylinder has been explored by St. Lawrence and Qamar [58], Ferrick et al. [31] and Ferrick and St. Lawrence [30]. Chouet [15] explored a complex model consisting of a gaseous chamber over a fluid-filled cylinder which is encased in an elastic half space. Crosson and Bame [18] used a spherical liquid-filled source to study the relation between the dominant frequency and size of the source. 1.3 O bjectives Our main goal is to interpret the observed characteristics of LP events at the Piton de la Fournaise described in 1.1 and chapter 2 by using a simple source model. V V e adopt the finite element method which allows a flexible choice of geometry of the source and velocity distribution of the surrounding area. It has a potential for incorporating the topography in ground motion calculation. Our model, however, is restricted to 2 dimensions because of the computational cost. Therefore we cannot make detailed comparison of the observed and simulated 11 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. seismograms. Rather, we shall extract essential features of the observations on LP events at the Piton de la Fournaise, and interpret them using a simplified model. We shall use a line pressure source to simulate the LP events, and try to determine the geometry and location of the source from the records of LP events in a similar manner as the location of the magma chamber has been estimated traditionally from static displacement using Mogi’s method. Our proposed modeling may appear straight forward and easy, but because of the simulation flexibility and functionality required. LP event simulation is actually very involved. Our modeling by the finite element method involves mesh generation, source simulator, matrix assemblage, matrix bandwidth reduction, linear equation solving and result visualization. In order to obtain reliable, numerical results, it is crucial to have a high-quality, flexible framework that can serve as a basis for extensions, and the design and implementation of generic software tools, which can be used in a wide range of applications. We base on the excellent framework of Diffpack [56] to develop a generic program to study seismic wave field. By specifying a generic interfaces, it allows users to extend the program to meet their own specific needs. This is another major goal I want to achieve in this thesis in addition to the seismological research of LP events. In seismological researches, wave propagation is the central subject with a va riety of applications such as fault rupture, strong motion, volcanic seismic signal. Although these problems seem quite different, its core is to study the wave field. The main difference among these problems is the source. If we implement a program that allow people to plug-in a different source simulator, it can allow reseaxchers to study all these problems with quite a ease without writing code from scratch. 12 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C hapter 2 Long P eriod E vents at P ito n de la Fournaise 2.1 Introduction Piton de la Fournaise (Fig 2.4) on the island of La Reunion in the Indian ocean is a basaltic volcano like Kilauea, Hawaii. Both are intra- plate volcanoes on the oceanic lithosphere with a supply of magma from hot spots in the mantle (Fig 2.1). The rate of erupted magma is not much different. An estimate for the Kilauea by Dzurisin et al. [24] for the period from 1956 to 1983 is 30 x 106m3 per year, while Lenat and Bachelery [62] [61] gives 10 x 106m3 per year for the Piton de la Fournaise during the period from 1930 to 1985. Since Dzurrisin et al. [24] estim ate that 65% of the magma supplied from the upper mantle remains in the rift zone under Kilauea, while there is no evidence for such a steady accumulation of magma in the rift zone of Piton de la Fournaise, the difference in magma supply from the upper mantle may be about an order of magnitude between the two volcanoes. The rift zone clearly exists in Piton de la Fournaise, as demonstrated by eruptions along the northeast rift zone outside the Enclos caidera in April,1977, and a series of fissure eruptions extending from the summit to the south coast in March, 1986. It is also clear that 13 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the rift zone is not developed here as well as in the Hawaiian volcanoes. From a sea beam bathym etry survey, Lenat at al., [62] found that the active rift zones of Piton de la Fournaise widen downslope unlike typical Hawaiian rift zones, which form narrow ridges extending tens of kilometers offshore. Nakamura [69] was the first to ask why long rift zones develop in Hawaiian vol canoes. Comparing with the Galapagos islands, another intra-plate basaltic volcano without a rift zone and located on a young lithosphere, he suggested that the rift zone intrusion in Hawaii is sustained by repeated Kalapana-type earthquake caused by a slip along the top of the oceanic crust which may have anomalously low friction due to a thick sediment and a high pore pressure ( Hubbert and Rubey, 1959). This idea has been supported by more recent works( [23], [88]). Piton de la Fournaise is located on an old oceanic lithosphere with the flexural rigidity (Walcott [89], Bonneville et al. [9]) comparable to the Hawaiian lithosphere, and the oceanic crust upon which the volcano was built must have had a thick sediment cover. There is, however, a major difference in the speed of lithosphere relative to the mantle under these two volcanoes. In Hawaii, the high plate speed always gives purely oceanic crust in front of the volcano, allowing the Kalapana-type earthquake along the weak plane due to the deep sea sediments, which sustains the intrusion into the rift zone. On the other hand, the African plate, on which Piton de la Fournaise is located, moves so slowly that a fluctuation in the path of magma ascent in the crust can move the volcano back and forth, allowing the remains of proto Fournaise volcano, as evidenced from the old magma chamber encountered by deep drilling (Rancon et al., [75]) as well as the gravity highs deviating from the current su m m it 14 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. area (Rousset et al., [77]; Malengray, 1996), to prevent the full development of the rift zone. Thus it is more difficult for magma to intrude into the rift zone here than in Hawaii, and that explains why Piton de la Fournaise has a remarkable central cone as high as 500m from the floor of the Enclos caldera due to repeated summit eruptions, which is absent in the Kilauea volcano. From the above geological and geodynamical setting, we expect that the erup tion at the summit of the central cone and that at the rift zone may have distinctly different seismic signatures. In fact, Aki and Ferrazzini (1998) [4], found a remark ably systematic relation between the elevation of the eruption site and the duration of the swarm of volcano-tectonic events under the summit area which occur as a precursor to an eruption. As shown in Table 2.1 and Fig 2.3, the precursor time ranges from 17 minutes to more than 9 hours, systematically longer for eruption sites at lower elevations. An important observation relevant to this thesis is that the LP events were observed during the precursory swarm only for the four events which clearly belong to the rift zone eruption from the location of their eruption sites. This suggests that the LP events are associated with motions of magma in horizontal direction but not with those in vertical direction. Piton de la Fournaise offers an unusual opportunity for a study of LP events because the volcano went through a variety of eruptive activities since early 1980’s during which a modem seismic network of currently 19 stations has been operated by the Institute de Physique du Globe de Paris (IPGP) for monitoring the seismicity associated with eruptions. The volcano was quite active from 1981 to 1992, during which 28 eruptions occurred (Fig 2.2). It erupted again on March 9, 1998. 15 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 2.1: Location of the Piton de la Fournaise of La Reunion Figure 2.2: Eruption of the Piton de la Fournaise of La Reunion Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. z/uu , 2500 - 2300 - V3 e o S 2100 - u 'o o c 3 > JO Hi 1900 A 1700 - 1500 ? s | s i ° - 5 ? ~ 8 S Summit £ 3 s s © £ : . i Intermediate • Highest Point of Eruption Site • With Long-Period Events £ M 'T I y 3 Rift Zone i n i I i i‘i i’i i i n 11 i'i t ri i r ri i i i i i i i i i r 11 i i m 0 1 2 3 4 5 6 7 8 Precursor Duration in Hour 3 a e 10 Figure 2.3: A relation between the precursor time of the seismic crisis and eleva tion of the eruption site. Only those crises preceding the rift-zone eruptions are accompanied by the occurrence of LP events 17 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 2.4: Piton de la Fournaise of La Reunion 18 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2.2 Identification o f long-period events Seismological investigations on eruptions at Piton de la Fournaise since the establish ment of the observatory in early 1980’s were described in journal papers by individual scientists ( Bachelery et al.,[6]; Lepine, 1987; Lenat et al., [61],[60]; Delorme et al. [21]; Hirn et al., [38] [39] [40]; Grasso and Bachelery [35] ) as well as in the annual report published in the Bulletin of Volcanic Eruptions by the observatory personnel. It was found that almost all the eruptions were preceded by the so-called “seismic crisis” which is a swarm of volcano-tectonic earthquakes occurring beneath the sum mit of the central cone at about the sea level. The precursor swarm lasts 17 minutes to more than 9 hours depending on the elevation of the eruption site as mentioned earlier. The seismic signature of the eruption is also quite regular characterized by steady continuous tremors originating from the eruption site. The eruption tremor is not purely harmonic, but contains broad spectrum of frequencies from about 1 to several Hz. The amplitude is the largest at stations closest to the eruption site and decays quickly with the distance from it. In addition to the above swarm of shallow volcano-tectonic events and the eruption tremors, there have been occasional mention of “long-period events”, but they have not received much attention partly because their occurrences axe not simply related to eruptions, and partly because of difficulty in discriminating them from other seismic events such as numerous rock falls along steep valleys in the island and P, S and T phases from earthquakes along active ridges in the Indian ocean. Since we cannot apply the conventional method for locating tectonic earthquakes to LP events because of the lack of clear onset of P and S waves, Aki and Ferrazzini 19 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (1998) [4] developed a method for discriminating a LP event from other events with similar seismic signature by using the spatial distribution of their amplitude across the network. For this purpose it is essential to eliminate the local amplification effect of recording site from the observed amplitude. Recently, Koyanagi et al. [50] have demonstrated that the coda amplification factor applicable to S waves (e.g., Kato et al., 1995) can be effectively used to smooth the amplitude distribution of T phase observed at the seismic network of the Island of Hawaii. The apparent dominance of S waves in seismic waves converted from the acoustic waves in the oceanic waveguide may be understood from the order of magnitude greater conversion coefficient from P to S as compared to that from S to P by localized heterogeneities ([2]; [91]). The coda method works well also at La Reunion, not only for T phases, but also for seismic waves from rock falls and LP events. In order to obtain the station amplification factor by the coda method, we need a large enough tectonic event which gives a good signal to noise ratio at a lapse time late enough after all the primary waves and forward scattered waves having passed the station. The choice of lapse time twice the S wave arrival time originally adopted for the central Asia by Rautian and Khalturin(1975) has been found adequate for most part of the world. Aki and Ferrazzini (1998) [4] found, however, that a longer lapse time may be required for La Reunion because of a zone of strongly localized heterogeneities associated with the presumed path of magma ascent and intrusion. This severely limits the availability of data for measuring the coda amplification factor in this low seismicity region, where the maximum magnitude is around 3.5. A particular mode of delayed trigger recording adopted at the observatory also limits the availability of high quality coda data. For these reasons, the amplification factor 20 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. • Sf-R m • DSR Figure 2.5: Stations at the Piton de la Fournaise, La Reunion measured for a single M = 3 event occurred at 13:28 on February 1, 1995 off the coast near the northern end of the island is used. The peak-to-peak amplitude averaged within the lapse time interval from 60 to 70 seconds is normalized to a summit station BOR (Fig 2.5) as listed in Table 2.2. The dominant frequency of coda in the above interval was 2 to 3 Hz at all stations. The above amplification factors were tested by applying them to a seismic event with known source location, namely, a rock fall at Cascade de Mahavel. Rock falls are very common sources of seismic signals in this island, partly because of heavy rainfall (annual average of several meters) and partly because of steep cliffs and 21 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. deep valleys originated from caldera formations. The depth of the valley of Cascade de Mahavel is about 800 meters. The dominant frequency was around 3Hz at all stations. The variation of the amplitude corrected for the station site effect is smooth enough for drawing iso-amplitude contours at intervals of a factor of 2 as shown in Fig 2.6. The center of contours agrees well with the site of the rock fall. The more rapid decrease in amplitude toward west than east is expected because the rock falls occurred on the eastern slope of the valley. The decrease of amplitude with distance r toward east can be fitted with the formula for surface waves, namely, where / is frequency, v is the wave velocity and Q is the attenuation quality factor. A good agreement with the observation is obtained for f = 3Hz, v = 1 km /s, and Q = 40 to 60. This range of Q agrees well with the value obtained by Koyanagi et al. (1995) for the summit area and rift zones of Kilauea, Hawaii from the amplitude decay of T phase. Aki and Ferrazzini (1998) first applied the above station amplification factors to LP events observed during a seismic crisis of December, 1991. This crisis is called an intrusion event, because there was no eruption of magma, although both seismic and geodetic observations showed activities similar to those preceding eruptions. The seismic crisis started with a swarm of shallow volcano-tectonic earthquakes with fre quencies higher than 10Hz beneath the summit at about 5:05 GMT on December 7, 1991. After a vigorous occurrence of these events for about 20 minutes, low fre quency motions appeared forming the background over which high frequency events 22 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. CaMeaWall or CBff Relative Aagriitude 10$ « T la Figure 2.6: Amplitude distribution of seismic waves in the frequency band from I to 3 Hz from a Mahavel rock fall 23 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. axe superposed. This low-frequency background reached a peak at 5:41 GMT. As the activity declines discrete LP events became distinguishable at about 5:56 GMT. An example of records of a LP event observed during the crisis is shown in Fig 1.1 The low frequency background and LP events must come from a common source be cause their amplitude distributions across the network were indistinguishable. Their waveforms are very similar among the three summit stations as shown in Fig 1.1, BOR, SFR, and DSR (Fig 2.5) located near the rim of the summit caldera. Their arrival times and amplitude after correction for the site effect are also very similar among them. This excludes the possibility that they may be due to rock falls near the summit, because we know that the latter generates waves of much higher fre quency which varies among the 3 summit stations in the arrival time, amplitude and frequency content. Thus, the similarity of arrival time, amplitude and waveform among the three summit stations becomes an important criterion for discriminating LP events from rock falls near the summit area. For their discrimination from rock falls outside the summit area and events outside the island, we need to calculate amplitude distribution. In general, LP events show the largest amplitude at the summit and are rarely observed at stations more than 10 km away from the summit. There are no obvious difference in amplitude distribution among events with different dominant frequencies. In order to estimate the depths of LP events, Aki and Ferrazzini (1998) [4] compared their amplitude distribution with that of volcano-tectonic events after band-pass filtering. According to A. Hirn of the IPGP, the latter were clustered beneath the summit area at a depth of about 2 km. We filtered the records of one 24 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. of the largest volcano-tectonic event occurred at 05:14, December 7, 1991 through a pass-band between 1 and 3 Hz. The amplitude contours for the band-pass filtered volcano-tectonic event (see Figure 2.6) is similar to those for LP events, indicating the location of the source of LP events is close to the hypocenter of the cluster of the volcano-tectonic event. Aki and Ferrazzini (1998) [4] extended the study of amplitude distribution to 31 LP events occurred from 1990 to 1995, and found that they are all quite similar to those of the December 91 crisis as shown in Fig 2.7. The mean and standard error of amplitude of LP events normalized to the mean amplitude of the three summit stations are given in Table 2.3. 2.3 Sum m ary o f observations on LP events at P ito n de la Fournaise An important observation relevant to the question of magma involvement in the source process of LP event may be found in their occurrence after the last eruption in August, 1992. The temporal pattern of the disappearance of LP events since 1992 shows a systematic dependence on the dominant frequency of LP events as shown in Fig 2.8. Those around 1 Hz disappeared in 1993, while those around 2 Hz persisted until 1995. Not a single LP event was found large enough for confirmation by the amplitude distribution method since the summer of 1995 until June 4, 1997. The above observation suggests a 2 tier magma reservoir, one generating 1 Hz LP events on top of the other generating 2 Hz LP events. Aki and Ferrazzini (1998) [4] suggested that with the withdrawal of magma back to a deeper reservoir, the former 25 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C a iio a W a U o r C S f f JL5 a C itd m WiMorCBff C • S— rtc Sarto* = : CiUam ViMarCaff / — v No— fagrt 5 U > l • W ir tr Sabo* = > C itim Watt or CHIT M r— tiwrt AayUmrtt / ft/2 , U Figure 2.7: Distribution of median amplitude of LP events normalized to the mean of the three summit stations; obtained from (a) 14 events occurred in the period from 1990 to April, 1991, (b) 11 events during the intrusion crisis of December, 1991, (c) 6 events after the crisis till April, 1995. (d) shows the distribution of the amplitude of a large volcano-tectonic event during the crisis of December, 1991, band-pass filtered between 1 and 3 Hz. The number attached to each station is normalized to the mean of the three summit stations. 26 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. itTOWffBii li ii ii li li l l99IJan-Af> 3 8 1 8 8 0 1 } n m n c a t BIEH I Q Q l A f w O c c BSflSd k ll li nasBBfaBa a n n n f i B i I M h O O ! « H 2 2 a 9 & 0 8 10 12 14 I ft 1.8 2 0 2 2 2 4 2ft 2 8 <0 Predom inant Frequency (Hz) Figure 2.8: The histogram of annual occurrence of LP events with a given dominant frequency. The vertical scale is common to all years, exception for 1991, for which it is 4 times greater than the rest even though the year is separated into two parts 27 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. is emptied first followed by the latter. The seismicity of volcano-tectonic earthquake during this period was very low and decreasing each year, yet maintained one half of the 1992 level in 1995. In fact, their activity is revived in 1996 with a swarm of deep earthquakes (focal depth 16 km below sea level) occurred about 5 km NW of the summit in September, 1996, followed by an increase in seismicity under the summit and an intrusion event accompanied by an inflation detected by tiltmeters and a swarm of shallow volcano-tectonic events on 26 November, 1996. The first LP events since 1995 confirmed to originate from beneath the summit area occurred on June 4, 1997. The dominant frequency was 2Hz. Thus, the last disappeared is of the same frequency as the first reappeared. The first LP events with dominant frequency 1 Hz reappeared only 5 hours before the eruption on March 9, 1998. Although LP events are absent during the seismic crisis preceding a summit eruption as mentioned earlier, they occur, during a summit eruption ( for example, eruption of September, 1985) and after a summit eruption (for example, eruptions of December 2, 1985: January to February, 1987; and August to September, 1992). Aki and Ferrazzini (1998) [4] interpret these LP events as due to diversion of magma flow from the summit-eruption path to the rift-zone reservoir. There was a peculiar swarm of LP events in January to April, 1991, which appears to have no direct relation to any eruption. It was distinct from other LP events in that the distribution of dominant frequency was symmetric centered at 1.4 Hz. This swarm was terminated at the time of the occurrence of one of the largest earthquakes in Madagascar, about 800 km away from La Reunion. This reminds us the coincidence between the commencement of an LP swarm under Mammoth 28 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 2.1: Seismic crises at the Piton de la Fournaise Number site Date Precursory time Elevation of the eruption 1 14-06-85 less than 1 hour 2500 m 2 10-07-85 - no eruption 3 05-08-85 2 hour 37 min 2100-2000 4 01-12-85 17 min. 2500-2100 5 18-03-86 9 hour 24 min. 1750-1720 6 02-06-87 - no eruption 7 19-07-87 2 hour 13 min. 2050-2000 8 07-02-88 2 hour 5 min. 2000 9 18-05-88 31 min. 2300-2000 10 31-08-88 2 hour 25 min. 2250-2100 11 14-12-88 4 hour 31 min. 2100-2000 12 18-01-90 47 min. 2500 13 18-04-90 6 hour 45 min. 1800 14 18-07-91 52 min. 2510-2400 15 07-12-91 no eruption 16 27-08-92 57 min. 2500-2100 17 26-11-96 - no eruption Table 2.2: Coda amplification factor at stations near the summit Site station Coda amplification factor BOR 1.00 DSR 0.36 SFR 0.75 CHR 0.72 NCR 0.94 NTR 0.70 Mountain, California and killing of trees by carbon dioxide reported by Hill [37]. It is interesting to note that the time also agrees with the occurrence of the Loma Prieta earthquake of November, 1989. Thus apparently some LP events generators are in a delicate balance and minor disturbances can start or stop their occurrence. Let us summarize now key observations that should be considered in our modeling study. 29 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 2.3: Table of normalized amplitude for 31 LP events from 1990 - 1995 Station Distance Normalized amplitude CHR NCR NTR 2.5km 5.0km 5.0km 0.66 ±0.12 0.44 ±0.10 0.53 ±0.15 1. The fact that LP events always occur during the precursory crisis for eruptions in the rift zone, but none during those for the summit eruptions suggest that the LP generator is related to the horizontal movement of magma, but not to its vertical movement. 2. The observed disappearance of LP events after 1993 (first 1Hz LP followed 2Hz LP) suggests that the 2 Hz LP generator is located deeper than the I Hz generator. This is confirmed by the first reappearance of 2 Hz LP followed by I Hz LP before the March, 1998 eruption. 3. Comparison of amplitude distribution of LP events with volcano-tectonic events band-pass filtered in the same frequency range as those of LP events suggests that both types of events are originated at a similar depth, namely around 2 km. 4. We found that amplitude distribution of LP events shows nearly the same pattern during the period from 1990-1995. We could not detect any difference in the amplitude distribution for events with different dominant frequency within the experimental error. 5. The amplitude, arrival time, and waveform axe nearly identical at three summit stations. This strongly suggests that the LP source is symmetrically located beneath the summit area. 30 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6. The mean and standard error of amplitude normalized to the mean amplitude at the three summit stations given in Table 2.3 will be directly compared with model simulation. 7. For some cases in which waveforms are relatively simple such as shown in Fig 1.1, we are able to estimate the difference in arrival time between different stations. For example, Fig 1.1 shows the delay time relative to the summit stations is 0.65 s at CHR, 2.00 s at NCR and 2.10 at NTR. 2.4 F in ite Elem ent M odel The structure of a volcano system is usually very complex. To simplify the LP source modeling, we shall focus on a single dike or sill in which magma injection or degassing happens. When magma injection or degassing happens inside a magma chamber, it will generate localized high pressure inside the chamber. This localized high pressure therefore propagate the pressure outwards, acting like the explosion source. In 2D case, to further simplify the modeling of the process of magma injection or degassing inside a dike or sill, an embedded line pressure source inside the involved dike or sill is utilized. We assume that the source pressure will increase linearly from 0 to a certain value. Since the magma chamber is usually considered as a mixture of fluid and solid, we use elastic material with high Possion’s ratio and low wave velocities to approximate the characteristics of this kind of magma body. 31 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ground Magma cham btr 2km EmbaddadLina Praaaura Sourca Absorbing boundary Absorbing boundary Absorbing boundary Figure 2.9: Finite element model of Long Period events We shall use the model of the magma chamber elongated vertically or horizontally (Fig 2.9) to study the effect of magma chamber orientation on the seismic amplitude distribution at the surface. In our simulation, the ground velocities are simulated for comparison with the ob servations. The surface synthetic seismograms at 0.5 km, 2.5 km and 5 km away from the center are compaxed with the observations at (SFR, BOR, DSR) 0.5 km,(CHR) 2.5 km and (NTR and NCR) 5.0 km away from the summit. The depth of the magma chamber is inferred around 2.0 - 2.5 km from comparison with earthquake swarms proceeding the eruptions. But we also try various depths to examine the effect of the depth on amplitude distribution. It is believed that 2.0 km/s for P wave velocity and 1.1 km /s for S wave velocity are reasonable values at Piton de la Fournaise within 2km beneath of the surface. The P and S wave velocities of the magma chamber are assumed to be in the rang e 32 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 2.10: Uniform mesh used in finite element modeling, which is better than other kinds of meshes of 0.8-1.2 km /s and 0.5-0.6 km /s respectively. V V e use the finite element model (Fig 2.9) to approximate an infinite model by introducing absorbing boundary at both sides and bottom to reduce the reflection from the artificial boundaries. The surface displacement, velocity and acceleration are recorded during the sim ulation for analysis. V V e used uniform mesh (figure 2.10), proper time step and small element size to minimize wave dispersion. 33 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C hapter 3 F in ite E lem ent M odeling Im plem entation ‘ ‘ To get works done. Get tools done first. -Confucius 3.1 Introduction to F inite E lem ent M ethod 3.1.1 Basic Concepts in Finite Element Method Finite element method is a numerical technique to solve engineering analysis prob lems. The first step is to break a complicated structure into smaller elements. Sec ondly, analyze each element based on the physical law using numerical computing techniques and assemble ail elements into a big matrix of algebraic equations. As the third step, this matrix is solved. In this thesis work with regard to the im plementation of finite element modeling, the boundary conditions and variational formulation, integration scheme and numerical error of finite element method are discussed in details here. 34 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Boundary Conditions and Variational Formulation In finite element model ing, we can identify two classes of boundary conditions, called essential and natural boundary conditions. The essential boundary conditions axe also called geometric boundary conditions; namely, in structural mechanics, the essential boundary condi tions correspond to prescribed displacement and rotations. The order of the deriva tive in the essential boundary conditions is, in a Cm~l problem, at most m — 1 and only these boundary conditions must be satisfied when invoking the stationarity of II, the total potential as defined in ( 3.1). The second class of boundary conditions , natural boundary conditions, are also called force boundary conditions; namely in structural mechanics, the natural boundary conditions correspond to prescribed boundary forces and moments. It is very desirable that the governing finite element equations include the es sential boundary conditions for practical applications. Variational method is ideal for this purpose. The essence of the variational approach is to calculate the total potential II of the system and to invoke the stationarity of II, ie. £11 with respect to the state variables. This approach provide a powerful mechanism for the analysis of continuous systems. The main reason for the effectiveness lies in the way by which essential boundary conditions can be generated and taken into account when using variational approach. The total potential II is also called the functional of the prob lem. When the functional contains the derivative of a state variable (with respect to a space coordinate) of order up to m, the problem is called a C m~l variational -problem. 35 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For a linear elastic continuous system, the potential can be presented as: n = \ J eTCedV - J UTf BdV - j U j f s dS F{ (3.1) " V V s * where C is the stress-strain matrix of the material, U is the displacement of the system, f B and f s are the body force and surface force of the natural boundary condition, U { and Us are the displacement at the essential boundary and c is the strain. V V e can get a generic representation of the above equation 3.1: KU = R (3.2) where [\ is the system stiffness matrix and R is the system nodal force vector. K = f BTkBdv (3.3) t; R — J N Tf Bdv (3.4) V where N is the shape function, k is the material matrix and B is the strain-displacement matrix. The detailed derivation of the equation 3.2 can be found in most finite el ement textbooks. For the linear elastodynamic system, the generic representation can be formu lated as: MU + C U + KU = R (3.5) 36 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where M and C are the mass and damping matrices; and U, if, if are acceleration, velocity and displacement vectors of the finite element assemblage. M and C of an element are M = j N TpNdv V C = J N TciVdv (3.6) (3.7) where p and c are material density and damping coefficient. The relation between the shape function and strain-stress matrix can be found in most finite element textbooks. Material matrix k for 2D plane strain is k = E ( l - u ) (! + *)(! -2 u ) 1 1 / 1 - 1 / — I 1 - 1 / * ■ 0 0 0 0 l-2i/ 2 ( 1- 1 /) (3.8) Im position o f Essential Boundary C onditions Now we consider the imposi tion of the essential boundary conditions in equation 3.1. Two widely used procedures axe available, namely the Lagrange multiplier method and the penalty method. Here only the description of the Lagrange multiplier method is shown below. For example, assume that we have a system represented as: Maa Ma b Mb a Mb b Ua U b + Kaa Kab Kba Kbb U a Ra U b Rb (3.9) 37 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the Ua are the unknown displacements and the U b are the known, or prescribed displacements. Using a integration schema (for example, Newmark) to represent U based on U, the equation 3.9 can be put into the form (3.10) By using Lagrange multiplier method to impose the essential boundary conditions Ub, the equation 3.10 now become fiaa f^ab Ua Ra 1 --- i r p 2r- O' f t . Ub Rb l v aa 0 1 i [\ab&b 0 / i 1 — i i (3.11) where I is a unit matrix. The equation 3.11 takes the essential boundary conditions into account. The actual implementation of essential boundary conditions is carried out at element level before assemblage. N ew m ark In teg ra tio n M e th o d To solve Equation 3.5, we use integration method ( Newmark method), which is represented as: (7*+A‘ = Ul + [(1 - 8)0* + 8Ut+M]&t Ut+At = C /'f + AtU 1 + [ ( i - a)U* + aU t+*‘]At2 (3.12) (3.13) where 8 > 0.5 and a > 0.25(0.5 + 8 )2. The complete algorithm using Newmark integration method is: 1. Form stiffness matrix K, mass matrix M and damping matrix C. 38 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2. Initialize U°, U°, U° 3. Select the time step size Af, parameter a and S 4. Form effective stiffness matrix K = K + a0M + aiC . 5. For each time step, calculate effective loads: # + At = + < 1 2 0 *+ azUt) + C{alUt +aA Ut +<*(/*) (3.14) 6. Impose the essential boundary conditions and obtain the displacement at t + A£ from the linear equations: KUt + = Rt+At (3.15) 7. Calculate accelerations and velocities at time t + At: Ut+^ = a0 (Ut+At - W) - a2 Ul - azU 1 (3.16) (jt+±t = (jt + a6[ft + a7(ft+At (3>17) 8. based on acceleration, velocity and displacement of the time step t + At, recalculate the effective load vector at t + 2At. Then repeat the step 5,6,7 and 8. where a o ^i, < 1 2 , 0 . 3 , 0 . 4 , 0 . 5 , 0 . 6 , 0 . 7 are: 1 °0 = . ,2 a A t c (3.18) 0 ° l ~ aA t (3.19) 39 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 a2 = . . aAt 1 a .3 = ------1 2a c (3.20) (3.21) 6 a4 = -----1 a (3.22) A ir 5 n as = T ( - - 2 ) (3.23) a6 = A f(l — a) (3.24) a7 = 8 At (3.25) N um erical Error In order to obtain an solution of a dynamic response with adequate accuracy, it is important to choose an appropriate time integration scheme. This choice depends on the finite element idealization upon which the integration scheme is operating. However, the finite element idealization to be chosen depends in turn upon the actual physical problem to be analyzed. It follows therefore that the selection of an appropriate finite element idealization of a problem and the choice of an effective integration scheme for the response solution are closely related and must be considered together. The finite element model and time integration scheme are chosen differently depending on whether a structural dynamics or a wave propagation problem is solved. The basic consideration in the selection of an appropriate finite element model of a structural dynamics problem is that only the lowest modes( or only a few intermediate modes) of a physical system are being excited by the load vector. If a Fourier analysis of the dynamic load input shows that only frequencies below lju are contained in the loading, then the finite element mesh should at most represent accurately the frequencies to about u jc o = 4wu of the actual system. There is no 40 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. need to represent the higher frequencies of the actual system accurately in the finite element system, because the dynamic response contribution in those frequencies is negligible. For the analysis of wave propagation problems, a large number of frequencies are excited in the system so that we need to use a sufficiently high cut-off frequency u > co to obtain enough solution accuracy. Identifying the cut- off frequency to be used is not easy. Instead of using the above approach to obtain an appropriate finite element mesh for the analysis of a wave propagation problem, it is generally more effective to employ the concepts used in finite difference analysis and the method of characteristics for establishing an appropriate finite element mesh and time step A t for the analysis. If we assume that the critical wave length to be represented is Lw, the total time for this wave to travel past a point is: tw = — (3.26) c where c is the wave speed. Assuming that n time steps are necessary to represent the traveling wave, we have A t = — (3.27) n and the “effective length" of a finite element given by Le = cAt (3.28) 41 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. This effective length and corresponding time step must be able to represent the complete wave travel accurately and is chosen differently depending on the kind of element idealization and time integration scheme used. Although a special case, the effectiveness of using ( 3.28) becomes apparent when a one-dimensional analysis is performed with a lumped mass idealization and the central difference method. For example, a bar is idealized as an assemblage of two-node truss elements each of length cAt. the exact wave propagation response is obtained in the solution of the model, because the spatial and temporal discretization are completely equivalent to integrating along the characteristics. It is also interesting to note that assuming a uniform bar free at both ends the time step A t given in ( 3.28) corresponds to the stability limit Tn/;r, and the nonzero (highest) frequency of a single unconstrained element is u)n. Hence, the most accurate solution is obtained by integration with a time step equal to the stability limit and the solution is less accurate when a smaller time step is employed! In more complex two-and three-dimensional analyses, the exact solution is gener ally not obtained, and Le must be chosen depending on whether the central difference method or an implicit method is employed for solution. If the explicit central dif ference method is used, a lumped mass matrix should be employed and in this case low-order finite elements in uniform meshes axe probably most effective. Using these elements we can choose Le to be equal to the smallest distance between any two of the nodes of the mesh employed, and this length determines A t as given by ( 3.28). However, for higher-order (parabolic or cubic)continuum elements the time step has to be further reduced, because the interior nodes “are stiffer” than the corner nodes. Also, if structural (beam,plate or shell)elements are included in the mesh, the time 42 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. step size At may be governed by the flexural modes in these elements so that the distances between nodes do not alone determine At. Since the condition is always that A£ < 2*, where T„ is the smallest period of the mesh, we need to establish a lower bound on Tn. This bound is given by the smallest period of any element, considered individually.measured over all elements in the mesh. For some elements, this period can be established exactly in closed form, whereas for more complex (distorted and curved) elements, a lower bound on may have to be employed. The choice of the effective length and hence time step A t is considerably simpler if an implicit unconditionally stable time integration method (including Newmark Method) is used. In this case, Le can be equal to the smallest distance between any two of the nodes that lie in the direction of the wave travel. Wave dispersion is another issue in the numerical simulation. For non-uniform grids, the exact solution of wave dispersion can not be solved. But for a uniform grid, the good estimation comes from the analysis of finite difference by Dalain (1986). On the uniform grid [19], the relative phase speed for the time approximation is c / 92 * = V 1 + l2 (3-29> where 9 = wAf. c is the phase speed and C o is the wave velocity. Similarly, the relative phase speed for the space approximation is 5 = v1 - n (3 - 3 0 ) where (3 = kAt. 43 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Hence error in the time approximation makes that high frequencies advanced with respect to the phase speed. Error in the space approximation tends to de lay the higher frequencies. In this thesis, we use the parameters which make the wave dispersion in the time approximation defined by Equation 3.29 and space approximation defined by Equation 3.30 under 0.003. 3.1.2 Finite Element vs Finite Difference Method Finite element and finite difference method have very close relationship. As an simple example, consider the analysis of the uniform bar with the governing differential equation A 2u + / = 0 (3.31) Using an equal mesh point spacing h and the central difference formula we obtain as the difference equation at point i = / (3.32) This same equation is also obtained using a finite element idealization of the bar, in which the stiffness matrix of a single bar element is 1 -1 - L 1 (3.33) Hence, finite difference equations can be viewed as stiffness relations in finite element procedures. 44 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. However, the major difficulty of using conventional finite difference procedure is the incorporation of the boundary conditions. Since in the analysis the differential equations of equilibrium of the system are approximated directly by the difference scheme, it is necessary to satisfy in the differencing both the essential and the natural boundary conditions. This could be difficult to achieve at arbitrary boundaries. Because of this difficulty in conventional finite difference, finite difference energy method is developed based on the principle of minimum potential energy to include the essential boundary conditions in the differencing. As might well expected, the finite difference energy method is very closely related to the displacement-based finite element method, and in some cases same algebraic equations are generated. The specific differences between finite difference energy method and finite element method lie in essentially in the choice of generalized displacement components. An advantage of finite difference energy method lies in the effectiveness with which the coefficient matrix of the algebraic equations can be generated. This is due to the simple scheme of energy integration employed. However, when comparing finite difference energy method with finite element method, it must be noted that finite difference energy method has been applied to specific problems and finite element programs can also be unusually effective for specific classes of problems. The main advantage of the finite element method is that the procedure can be used effectively in general-purpose analysis programs. This feature is extremely valuable to engineering and scientific research. Thus we prefer using finite element method for scientific computing. 45 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.2 O bject-oriented scientific com puting 3.2.1 Object-oriented technology To realize the simulation by finite element method, we must choose a way to build the system. The traditional way is the procedure-oriented technology. This development scheme is straight forward and simple. Also due to the limitation of the computing power in the past, this is the predominant way to do scientific computing. Another way is called Object-oriented technology. With the dramatic development of com puting technology, the computing power has increased so much that now a personal computer can carry a rather large simulation job. Object-oriented languages also are becoming the predominant languages for software development in the indus try. Therefore object-oriented technology is becoming one of the most important technologies in software development. It has enormous potential for significantly increasing programmer productivity and code maintainability. Why is object-oriented technology gaining such a widespread acceptance among commercial software developers and large corporations? Procedural code of indus trial strength (thousands of lines of code) is often unstable under small perturbations, e.g., trying to fix a single problem might create five new ones. When the code size is larger than some threshold (about 105 lines of code), procedural techniques appear to break down: design and implementation complexity makes projects quite diffi cult and expensive, and eventually maintenance cost dominates development cost. Object-Oriented techniques promise to allow more stable and larger applications (several millions and tens of millions of lines of code) to be built. How can it be done? Why should a natural scientist care? 46 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Natural scientists model and simulate natural phenomena. Computers simulate “real” (e.g., balancing a checkbook, mimicking the behavior of a physical system,...) and “imaginary” (e.g., games) worlds. Rapid advances in hardware and software are making possible the treatment of interesting and complex systems (e.g., tur bulent flows, critical dynamics, lattice quantum chromo-dynamics, weather predic tion, cellular metabolism, social phenomena,...). Simula, a language that pioneered object-oriented technology was specially conceived for simulations. Hence, it is not surprising that object-oriented programming may make simulation and visualiza tion of complex phenomena easier than procedural techniques (e.g., Fortran). For fully exploring models of complex systems, displaying results and data in a powerful manner (e.g., scientific visualization, multimedia, virtual reality, ...), for enhanc ing collaboration fully exploiting networking facilities, object-oriented technology appears better suited than procedural languages. The principle “divide and conquer” has been known for a long time as an efficient technique to conquer large empires (e.g.. large codes). But, how to divide? Both procedural and object-oriented programming employ the principle. The fundamental difference between these two models is the choice of building blocks. In procedural programming, procedures are the fundamental units (usually called subroutines or subprograms). In the object model, basic units are “cells” or “atoms” called objects which contain both data (i.e., state variables associated with the “atom” state at a given time) and methods (i.e., dynamical rules, rules that explain how “atoms” interact in the outside world). In object-oriented terminology, objects encapsulate data and behavior. These objects are individual weakly interacting blocks. Objects interact (exchange messages) to produce collective behavior. 47 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Interactions axe mediated by messages that objects exchange with each other. An illustrative example for this kind of interactions can be shown by stiffness matrix formation in our system. Object ElasticWave wants to assemble the system stiffness matrix. To do that, object ElasticWave needs to get the stiffness matrix from each element in the system. There are many different types of elementfe.g. ElrnB4n2D, ElmB9n2D, EimB3n2D). Each kind of element is quite different from the others. But to the object ElasticWave, only behavior it is interested in each kind of element is to expose the stiffness matrix to the ElasticWave object. So by defining an interaction protocol ( interface), a collaboration can be achieved. Because of this separation, “divide and conquer” becomes much easy, and the modularity of the system is much better. At first glance, object-oriented design may seem straightforward: just mimic the objects of the outside world and their interactions, relevant to the piece of reality to be captured in softwaxe. However, this is not the case in real object-oriented design. In general, reality is too complex to be fully captured in a design. On the other hand, many of the outside objects/interactions may be totally or partially irrelevant to our problem. In the design of a complex system one has to identify key abstractions, propose some candidate classes and objects to mimic the system, and conjecture some class taxonomy. Most time spent on the development of our system is on the identification of key abstractions, which will be discussed later. 48 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.2.2 Object-oriented language C + + vs Procedure-oriented languages Object-oriented systems need to be implemented in object-oriented languages. Ba sically there are two kind of languages in the computer sciences. One is called procedure oriented language. Well known FORTRAN and C are both procedure- oriented languages. FORTRAN 77 is simple, easy to learn, and indeed suitable for many small size scientific programs counting only a few thousand code lines. To most researchers, they really don’t care about the languages as long as they get jobs done. In the past, researchers in the academic fields mostly deal with simple or simplified problems. Either these problems have analytical solutions, or numerical simulations are not very complex. Therefore the programs were usually not very big. So in scientific computing, FORTRAN 77 is dominating because of these features. The problem with FORTRAN 77 for large-scale problem (for example, 2D or 3D wave propagation simulation) are two-fold. First, programs using FORTRAN 77 axe prone to error because its compiler cannot detect many of frequent typing error and cannot do the run-time type-checking. Researchers probably will spend much more time verifying the programs than using C + + [12]. While C + + compiler will help you detect a lot of errors and save the duplication of code. Second, large FORTRAN 77 programs make very strong demands on the devel opers in order to avoid a spaghetti-like design as time evolves and function grows. C + + make it much easy to design a system with good extensibility and modularity because of its object-oriented nature [8]. 49 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FORTRAN 77 is just an example of such procedure languages. FORTRAN 90 has added a lot of useful features to FORTRAN language. Matlab also has many interesting computing features and graphic functions. But still the expressiveness, brevity and modularity of those languages does not even come close to C + + . Long Period event modeling is a very complicated one. To simulate the seismic wave from volcano activity using finite element method, is going to involve mesh generation, source simulator, matrix assemblage, matrix bandwidth reduction, lin ear equation solving and result visualization, which will be discussed in details in later section. Each component is fairly complex and big. Using FORTRAN to im plement the system will have two major problems. One is difficulty to design a modularized system. Second is difficult to utilize other people’s codes and be uti lized by other people. Especially the first problem of using FORTRAN will lead to hard time changing code for new physical models, which already hunt me during the development of my FORTRAN code. By using C ++ or other object-oriented languages, we can overcome these prob lems. However,it requires researchers ( including seismologist) to spend fair amount of time to understand object-oriented concepts as well as the object-oriented lan guages and operating system. 3.2.3 Objective and motivations Today, researchers and engineers have access to ever growing computing power and software that can solve numerical problems which are not fully understood in terms of existing theory. Thus, computational sciences can be viewed as an experimental discipline. A firm foundation for this experimental view is therefore crucial in order 50 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. to obtain reliable, numerical results. As a consequence, there is a demand for high- quality, flexible software that can serve as a basis for future extensions, and the design and implementation of generic software tools, which can be used in a wide range of application fields. In the coming decade, as the world comes to rely more and more on programming to solve real-world engineering, scientific and social problems, the importance of new languages, tools, environments, and compiler technology to support scientific programmers will increase rapidly. By focusing attention on practical aspects of this emerging technology , scientific programming will become mandatory reading for all researchers in their fields. In the past, people in the academia have spent considerable time and energy on the scientific computing field and made significant achievement using language C or FORTRAN. But one pitfall behind those success stories is that huge amount of work is redundant and a lot of nice efforts are wasted because of their poor extensibility. Most researchers have to spend (waste) a lot of time and energy to write the similar kind of codes over and over again. Also to most researcher, due to limitations of their computer knowledge, they probably cannot write good code which make it easy to allow their works extended because of lack of code generality and code modularity, even though the ideas behind the codes are excellent. It will be very beneficial to academic researchers if we can build a free well-tested common library of generic numeric codes, so other researchers can have access to the library and build their application on it very easily with high quality. Also it will make it easy to accumulate efforts from a lot of people working in the same area. 51 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. In seismological researches, the wave propagation is the central with a variety of applications such as fault rupture, strong motion, volcanic seismic signal. Although these problems seem quite different, its core is to study the wave field. The main difference among these problems is the source. If we implement a program that allow people to plug-in a different source simulator, it can allow researchers to study all these problems with quite a ease without writing code from scratch. Diffpack is an excellent scientific computation framework developed by Oslo Uni versity. It consists of a collection of object-oriented libraries for solving partial differential equations, and several Unix utilities for general software management and numerical programming. In particular, this piece of software is aimed at rapid prototyping of simulators based on PDEs, while offering a high level of efficiency. Based on Diffpack, I developed a generic application program to study wave field of Long Period events. By specifying a generic interfaces, it allows users to extend the program to meet their own specific needs. This is what I want to achieve in addition to the seismological research of LP events. 3.3 Introduction to Diffpack Diffpack consists of a collection of object-oriented libraries for solving partial dif ferential equations, and several Unix utilities for general software management and numerical programming. In particular, this piece of software is aimed at rapid pro totyping of simulators based on PDEs, while offering a high level of efficiency. The libraries, which are implemented in C + + , are organized into several layers: 52 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. • BasicTools: Non-numerical utilities such as memory handling, basic array structures, menu system and generic I/O. • LaTools: Class hierarchies for representation of different vector and matrix formats and linear algebra tools, such as direct and iterative solvers. • DpKernel: Basic building blocks for finite element programming such as rep resentation of grids, fields, element types and numerical integration. • DpUtil: High level utilities for finite element programming including prepro cessors, bandwidth reduction, filters to different visualization programs and easy-to-use interfaces to low-level classes. Diffpack offers a set of C-|— classes to support finite element programming in Object-oriented fashion. These major classes for finite element programming are described briefly here: F E M Class FEM offers the standard finite element algorithms for the element- loop, numerical integration over elements and smoothing of discontinuous deriva tives. Usually, FEM is used in conjunction with a simulator class and serves as basis class for the simulator. All the standard finite element algorithms are hence inherited in the simulator class. Some virtual functions must be implemented in the simulator class (typically setting the essential boundary conditions, and defining the integrands in the weak form of the boundary value problem). The reports Getting Started with Finite Element Programming in Diffpack and in particular Details of Finite Element Programming in Diffpack give information on the application and purpose of class FEM. 53 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The makeSystem functions that computes the linear system arising from appli cation of the finite element method, measures the CPU time spent in the function. This CPU time is written to the log-file (CaseName.dp) if more than half a second is spent in the makeSystem function. Similar CPU time information is provided by LinEqAdm::solve and written on the same line in the log-file. The CPU time is available in an internal (protected) variable that the programmer can use in a derived simulation class if desired. D egreeFE The class relates degrees of freedom in a linear system to the nodal values of discrete finite element fields. By degrees of freedom we mean the unknowns in a linear system associated with some finite element discretization of a system of partial differential equations. For example, when solving a scalar partial differential equation (e.g. Poisson’s equation) by the finite element method, there is usually one degree of freedom associated with each node in the finite element mesh. When solving a vector equation in d dimensions there are d component fields. The relation between, for example, component no j at node k and the corresponding global degree of freedom in the linear system is computed by a DegFreeFE object. The class is simple to use for both isoparametric and non-isoparametric element meshes. We refer to the document A Brief Documentation of Finite Element Fields over Non- isoparametric Elements in Diffpack for the definition of basic concepts used in and the general ideas behind class DegFreeFE. This document defines the important concepts general and special numbering used below. The DegFreeFE class may also mark the degrees of freedom that are subjected to essential boundary conditions. Class ElmMatVec uses DegFreeFE for this purpose. 54 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Moreover, the bandwidth of the coefficient matrix in the associated system of linear equations can be calculated. Note that it is assumed that the size of this system equals the total number of degrees of freedom in the associated DegFreeFE object. That is, if a vector equation is solved sequentially, component by component, a DegFreeFE object with one degree of freedom per node should be used for each component equation. GridFE The class represents a finite element grid with information on nodal co ordinates, the element topology and the nodes that are marked with boundary in dicators. Class GridFE objects are usually declared by an empty constructor and filled by preprocessor classes. For preprocessor programmers and for educational purposes direct use of the assignment functions in the GridFE class is necessary. Finite element computations require access to the read-only functions in class GridFE. 3.4 Long P eriod E vent Sim ulation S ystem By using the modules provided by Diffpack for finite element programming, we build a system with great modularity and ease of use. The modularity allow people to extend the system in a systematic object oriented way to meet the increasing needs for the system. For example, if we want to use a different grid generator, we can easily set up a new object type which specify some defined interface so that the extension is effortless and proper. The ease of use will benefit users from the extensive functionality and comprehensive documentation of the system. 55 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Grid Generator Boundary Conditioner Material Property Handler ModelPack G ) ® £ 5- I © CO 3 c £ o a © C O 3 © * o o a Long Period Event Simulation System Figure 3.1: Long Period Event Simulation System The system consists of three main components (Fig 3.1): ModelPack, Generic Elastic Wave Simulator and Result Report. 3.4.1 ModelPack ModelPack is the core of the system. It consists of three sub-components: Grid Generator, Boundary Conditioner and Material Property Handler. Grid generator not only generate the grid, but also mark the nodes in the grid with boundary indicators which help the system to set the boundary conditions in Boundary Conditioner. There are four class of grid generators that you can plug into the system. • PreproBox in Diffpack 1.4 [56], simple and easy to generate rectangular and triangular element in ID, 2D or 3D without holes. PreproBox is good for the problems with simple geometries. 56 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. • PreproSupElSet in Diffpack 1.4 [56], more complex and powerful than PreproBox. PreproSupElSet allow user to create a grid which consists of a set of sub- domains. The grid of each sub-domain is generated by PreproBox with its own specifications. By combining a set of grids generated by PreproBox, we can create more complex grid. • GeomPack in Diffpack 1.4 [56], developed by Barry Joe [44] [45] [46] [47], very complex and powerful than the above two grid preprocessors. By using GeomPack, user can create a grid with holes, but the input is rather complex [56]. • CrackPack, for handling special grids. It is a product by personal communi cations between myself and H. P. Langtangen who is one of chief designers of Diffpack(I promised him a report on this CrackPack development) and com prehensive study of the source code in Diffpack. CrackPack allow users to create a special grid which some nodes have same coordinates but belong to different elements. This feature is crucial to simulate the explosion source and crack. The grid generated by these grid preprocessors usually do not lead to formulation of a coefficient matrix with the minimum bandwidth. By re-numbering of the nodes in the grid using band reduction algorithms provided in Diffpack, we can reduce the bandwidth of the coefficient matrix to save the storage size and computation time later. This is a good example how a good and comprehensive framework can benefit the system development. 57 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. B oundary C onditioner associate the boundary indicators with the actual bound ary condition value, either displacement or nodal forces. This component handle very flexible boundary conditions. For example, one kind of boundary condition is time- dependent. User can specify the time history for one boundary condition of either displacements or forces by simple inputs. In summary, this component can handle the following boundary conditions: • Fixed displacement boundary condition • Fixed force boundary condition • Time dependent displacement boundary condition • Time dependent Force boundary condition M aterial Property Handler use a very sophisticate method to assign material properties to each element in the simulation system. In this component, several material regions are defined independent of the element meshing. Each region has its coverage, z-index and interpolation of material property with respect to spatial location. Z-index is a term used in Computer Graphics to determine the visibility of objects. We extend the concept as described in the following. By checking whether the center of an element is inside a material region, we can determine the material property of this element. If the center of an element is inside multiple regions ( regions overlap), the concept is that the material region with greater z-index con tains this element. Inside a material region, the material property can be same or interpolated by the spatial location. By allowing the property to be interpolated in 58 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. a region, it becomes very easy to define a region whose velocities increase with the depth. These two features are very useful for people to construct a model with complex material property distribution with combination of several material regions. In Summary, the material region can be the following: • Rectangular shape with same material property; • Rectangular shape with material properties defined by a function of location; • Elliptic shape with same material property; • Elliptic shape with material property define by a function of spatial location; • Polygon with same material; • Polygon with material property define by a function of spatial location; The material property consists of parameters : mass density p. damping coeffi cient, elastic coefficient vectors. The isotropic elastic coefficient vector has Young’s Modules E and Poisson’s ration u. The anisotropic coefficient vector can be repre sented as Cijim,i.j,l,m = 1,2,3 The anisotropicial is useful to study the wave-split in seismology. In addition to provide the great flexibility in handling material property, we also make the specification of such material properties easy and modular in a consistent way so that it is easy to extend our system to add their own additional property and interpolation of material property in a region. An example of constructing a model with complex velocity structure can be shown in the following. Suppose we want to define a model which contains a basin. 59 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Material Region 1 zrindexsj _ .... / Material Region 1 z-index=1 Figure 3.2: Schematic diagram of constructing a model with complex velocity struc ture In the basin, the velocity structure is different from the one outside the basin. To achieve this, we define two material regions. The first is an elliptic region and the second is a rectangular region. Each of these two material region can have different velocity structures. The first region has z-index 1 while the second region has z-index 0. The resulting velocity structure is shown in Fig 3.2. Inside this component, there are four major parts. • Initialization - initialize the system to a proper state before carrying any cal- • Calculation of coefficient matrices M, C and K of each element; For linear systems, matrices M, C and K are time-independent. We take advantage of this fact to speed up the simulation. 3.4.2 Generic Elastic Wave Simulator culations; 60 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. • Applying Essential boundary condition - By utilizing Degree of Freedom Man ager DegFreeFE, the essential boundary conditions are included in the linear equation m atrix during matrix assemblage. • Solving the equations -This is the core of this component which provide effi cient solution of wave simulation. By using generic Linear Equation Admin istration - LinEqAdm, users can plug in different linear solvers. Currently users can use direct solvers ( Gaussian Elimination Method, Algebraic Multi grid), iterative solvers (Jacobi Iterations, Successive over-relaxation SOR and Symmetric successive over-relaxation SSOR), or Krylov subspace methods. For elastic wave propagation problems, we can improve the efficiency greatly by taking advantage of the fact that these coefficient matrices are time independent ([52],[54]). Therefore the coefficient matrices are calculated once, then are used in the later computation. 3.4.3 Result Report This component export the simulation results into different formats that can be visu alized in a variety of tools ( for example, Matlab, Gnuplot etc). Detailed description can be found in technical reports of [55], [11] and [53]. 3.4.4 Docum entation and future development All these system design and implementations described above are useless without extensive documentation. The documentation not only provide users the only way 61 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. to learn to use the system, and give the developers to track the system development. Therefore, the documentation is a crucial part of the system. Our documentation consists of two parts. The first part is to provide the doc umentation for explanation of the system design. This part is for tracking system development and further development. The second part is to provide the examples how the system can be used. Although the work needed for the documentation is huge and a long time ef fort , we realized that this is very important. This thesis will provide part of the documentation, and we will put more usage examples on World Wide Web. 62 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. C hapter 4 N um erical Sim ulation R esults '“ Just do it. -Nike” We now describe results of applying the computer program described in Chapter 3 to the model described in Chapter 2 in order to model LP events at Piton de la Fournaise. 4.1 H om ogeneous half-space Line Pressure Source The first case investigated is a homogeneous half-space with a line pressure source buried at depth of 2.0 km. The time-dependence of pressure of line source is a ramp function, increasing from 0 to a certain value linearly for some time(called rise time) and stays constant afterward. This is to approximate the degassing or injection process inside the magma chamber. The resultant waveforms of vertical component of ground velocity at surface at various epicentral distance are shown in Fig 4.2. 63 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. In this thesism, we measured the maximum peak to peak amplitude of synthetic LP event at each stations to calculate the normalized amplitudes to compare qual itatively with the observation. To compare quantitatively with the observation on normalized amplitude distribution shown in Table 2.3 and Fig 2.7, we convert 2D normalized amplitudes into approximate equivalent 3D normalized amplitudes using ray theory. For our two dimensional model with the line source at depth d km, the vertical surface ground velocity at epicentral distance x can be represented in polar coordinates (r,9) as V 4 dW = ^ (4.1) v r where r = \/x 2 + d? is the distance between the source and the station, and 0 = sin~l(x/r). f(0) = s(0)t(0) is a function of 6 representing the amplitude adjustment factor, which consists of two contributing factors. The first one s(0) is the the source radiation pattern of LP event. The second one t{0) is the amplitude of vertical free surface ground velocity caused by an incident cylindric wave of LP event with a unit of amplitude and angle of incidence 0. If we normalize amplitudes to the one of the station at epicentral distance xa km, the normalized amplitude at epicentral distance x is ^ (I)= ( 4 -2 ) where ra = y d 2 + x2 and 0O = sin l(x0/r0). 64 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For a 3 dimensional model with the pressure source at depth d km, the verti cal surface ground velocity at epicentral distance x can be represented in spherical coordinates (r,0,<p) as V3D(x) = (4.3) r where F{0,o) = S{0,o)T{0) represents the amplitude adjustment factor, which also consists of the two contributing factors. The first one S(0,4>) is the source radiation pattern. The second one T{0) is the amplitude of vertical free surface ground velocity caused by an incident spherical wave of synthetic LP event with a unit of amplitude and angle of incidence 6. which is independent of d > . Similarly, if we normalize amplitudes of the 3D model to the amplitude of the station at epicentral distance xp km, then the normalized amplitude at epicentral distance x is F ( d ’ <t>)r P ( A JX V 3£>(*) = (4-4) r (0p,0p)r where rp = yjd1 + x and 0P = sin~l(xp/rp). To convert the 2D normalized amplitudes to the equivalent 3D normalized am plitudes, we need to find the relation between F{0,<p) and f{0). Because a cylindric wave can be viewed as integral of spherical waves along one axis, for a cylindric wave and a spherical wave incident upon the free surface at a same incident angle, the ratio of t(0) and T(0) is independent of 0, if they are of a same type wave and a same amplitude. Specifically for the cylindrically symmetrical source in 3 dimen sions, which is corresponding to the bilaterally symmetrical source in 2 dimensions, the source radiation pattern S(0, < p ) is independent of c f > and the same as the 2D 65 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. source radiation pattern s(0), Therefore F(0,<j>) = constant * f(9). Thus, the nor malized amplitude based on epicentral distance xp at epicentral distance x can be represented by the normalized amplitude 3o( 1 " V 7 (4-5) In this way, we can convert the normalized amplitudes obtained by our 2D fi nite element model to the approximate equivalent 3D values. Since the normalized amplitudes observed (Table 2.3) are based on the stations (SFR,BOR and DSR) at epicentral distance 0.5 km, for each case studied in this thesis, we converted 2D normalized amplitudes to approximate equivalent 3D normalized amplitudes using xp = 0.5. By comparing arrival time of the main P wave peak, we obtained the arrival time difference between different stations. V V e find that the time delay between receiver at 0.5 km and 2.5 km is 0.6 second, and the time delay between receiver at 0.5 km and 5.0 km is 1.8 second, and is in a good agreement with the observation. In order to see the effect of the depth of line source on arrival time, we changed the depth from 1.0 km to 4.0 km at interval of 1.0 km. The results is shown in Fig 4.1, which indicates a better agreement between the observation and calculated times for the shallower depth. On the other hand, the calculated amplitude normalized to the summit station is much smaller than the observed listed in Table 2.3, the agreement is poorer for shallower source. Thus, we conclude that the line pressure source embedded in a homogeneous half-space cannot explain the observation on amplitude and arrival time simultaneously. The line pressure source embedded in 66 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. depth 1.0km D i s t a n c e d e p t h 3 . 0 k m D i s t a n c e 6 5 4 2 1 0 0 2 4 6 8 10 depth 2.0km 6 5 4 ® E 3 2 1 0 2 4 6 8 10 D i s t a n c e d e p t h 4 . 0 k m D i s t a n c e 6 5 4 ® E (- 3 2 1 2 0 4 6 8 10 6 5 4 ® E 3 2 1 0 2 6 8 10 4 Figure 4.1: The comparison between the observed arrival times and those for a line source at different depths a homogeneous half-space is too simplified to explain observations on LP events at Piton de la Foumaise. The lines shown in Fig 4.2 corresponds to the arrival times of P and S waves calculated by the ray theory. It is clear that the records are dominated by P waves in the case of a pressure line source. 67 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -6 -4 -2 0 2 4 6 D i s t a n c e ( k m ) Figure 4.2: Vertical surface ground velocity for the case of a line source at depth 2.0 km Table 4.1: The normalized vertical surface ground velocity amplitude distribution of the homogeneous half-space model with line pressure source at different depths, the bracketed values in b o ld are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km.__________ _________ ___ Depth (km) Distance (km) 0.0 0.5 1.5 2.5 3.5 4.5 5.0 1.0 1.0 0.90 0.48 0.19(0.14) 0.14 0.09 0.07 (0.04) 2.0 1.0 0.91 0.64 0.45(0.40) 0.32 0.23 0.18 (0.12) 3.0 1.0 0.98 0.82 0.64(0.57) 0.41 0.27 0.23 (0.17) 4.0 1.0 0.99 0.92 0.67(0.63) 0.54 0.42 0.33 (0.26) 68 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6134 4.2 V ertically elongated m agm a chamber We shall now consider a low velocity magma chamber enclosing the pressure source with a vertically elongated shape. The line pressure source is embedded at the center of the magma chamber at depth of 2.0 km. The material property of active magma is believed to be mixed solid and liquid with Poisson’s ratio higher than 0.25. As mentioned earlier, we assume that P wave velocity 1.1 km /s and S wave velocity 0.518 km/s for the magma chamber. We assume also that the magma chamber is 2.0 km high and 0.5 km wide. The resultant waveforms of vertical component ground velocity at the surface at various epicentral distances are shown in Fig 4.3 and 4.4. The most distinctive difference from the case of the line source, in a homogeneous half-space is that the result shows the strongly trapped waves around the epicen ter ( epicentral distance 0 - 1.5 km). While at 5.0 km away from the epicenter, the normalized amplitude becomes very small. With regard to the amplitude dis tribution, compared to the homogeneous half-space case, the discrepancy between the observation and the simulation becomes even larger. The vertically elongated magma chamber magnifies the amplitude of the center station partly because of the trapped wave, therefore the normalized amplitude ratios at station 2.0 km and 5.0 km become smaller. Table 4.2 shows that comparison of the amplitude between the case of line source and the case of vertically elongated magma chamber at various epicentral distances. The depth of line source is 2.0km for both cases. Table 4.2 also shows that the amplitude for the case of horizontally elongated magma cham ber with the same dimension and the same depth from the result obtained in the next section (Fig 4.5). It is clear that the horizontally elongated magma chamber 69 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 4.2: The normalized vertical surface ground velocity amplitudes for the follow ing cases: (1) a vertically elongated magma chamber (2) homogeneous half-space, (3) a horizontally elongated magma chamber (this case is studied in the next sec tion). The bracketed values in bold are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km.__________________________________ Distance (km) Vertically elongated Homogeneous Half-space Horizontally elongated 0.0 1.00 1.00 1.00 0.5 0.64 0.91 0.94 1.5 0.57 0.64 0.98 2.5 0.29 (0.36) 0.45 (0.40) 0.97 (0.82) 3.5 0.22 0.32 0.61 4.5 0.14 0.23 0.36 5.0 0.11 (0.11) 0.18 (0.12) 0.33 (0.21) is the preferred model. Therefore in the next section, we will focus our numerical simulation on horizontally elongated magma chamber. The lines shown in Fig 4.3 are again the arrival times of P and S waves calculated by the ray theory. The inclusion of a low-velocity magma chamber generated not only trapped waves in and near the chamber, but also converted S waves propagating outside the chamber. 70 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. T i m e ( s ) Surface vertical velocity of vertically elongated magma chamber (0.5 x2.0 km) -2 0 D i s t a n c e ( k m ) Figure 4.3: Vertical surface ground velocity for the case of a vertically elongated magma chamber with height 2.0 km and width 0.5 km 71 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 001 535353535323232348485348485353 x 1 Q "5 Surface vertical velocity of vertically elongated magma chamber (0.5 x2.0 km) i k m & ■ 10 2 . 5 3 . 5 5 . 0 k m 0 1 2 3 5 4 6 10 7 8 9 Figure 4.4: Amplitude comparison of vertical surface ground velocity for the case of a vertically elongated magma chamber with height 2.0 km and width 0.5 km 72 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.3 H orizontally elongated m agm a cham ber In this section, we study the case in which the low velocity magma chamber is horizontally elongated, and the line pressure source is located at its center. The material property of the magma is the same as in the case of vertically elongated magma chamber. We shall investigate how the depth and the dimension of the low velocity magma chamber affect the ground motion at the surface. 4.3.1 Effect of the width First we look at how the width (horizontal extent) of the magma chamber might affect the surface ground motion of the model while the height (thickness) of the magma chamber is fixed at 400 m. We consider the cases in which the width is 0.577km, 1.04 km, 1.62km and 2.1km. The resultant vertical component ground velocity is shown in Fig 4.6, 4.7, 4.8 and 4.9, respectively. In section 4.1, we found that the observed arrival times of the LP events can be fitted by a line source in a homogeneous half-space. In that case, the LP events are primarily P waves as shown in Fig 4.2. But in case of horizontally elongated magma chamber, amplitude of synthetic P waves decay quite quickly with distance. In the synthetic seismograms (Fig 4.6, 4.7, 4.8, 4.9), S waves becomes dominant at a distance greater than about 3 km. Table 4.3 summarizes the dependence of amplitude distribution on the width of the magma chamber. The width of the magma chamber have significant effect on the amplitude distribution. When the width is 0.55 km, the magma chamber is approximately of a square shape. In this case, there are no strong S waves in the 73 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. synthetic seismograms. and the normalized amplitude is based on P waves. When the width is 1.04 km, 1.62km and 2.10km, strong S waves appear in the synthetic seismograms at distant stations, the normalized amplitude is based on S waves. In case of width 0.577 km, we found that the distribution of the normalized amplitude is close to the distribution of the normalized amplitude in the homogeneous half space case, which does not explain the observation. When the width increases, the normalized amplitude does not simply increase at all locations. For example, when the width is changed from 1.04 km to 2.10 km, the normalized amplitude at 1.5 km decreases from 0.95 to 0.80, while the normalized amplitude at 5.0 km increases from 0.20 to 0.50. The wider low velocity magma chamber makes more flattened amplitude distribution. Comparing the converted 3D normalized amplitudes listed in Table 4.3 with the observed normalized amplitude listed in Table 2.3, we suggest that the width of the low velocity magma chamber is around 2.0 km. Since P waves still dominate at shorted distance, the arrival time measured on the synthetic records no longer agree with the observed. To reduce the discrepancy on arrival time, we vary the depth of the magma chamber. A better agreement is obtained when depth is around 2.5 km (Fig 4.12). The calculated arrival delay time with respect to the summit stations is 0.61 second at 2.5 km and 2.2 second at 5.0 km compared to the arrival delay time observation 0.65 second at 2.5 km and 2.1 second at 5.0 km. This result will put the magma chamber slightly below the earthquake swarm occurring before the eruption. 74 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 4.3: The normalized vertical surface velocity amplitude distribution of the horizontally elongated magma chamber with fixed height 400m and different width, (*)the normalized amplitude is based on P waves. The bracketed values in b o ld are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km.______ ___________________________ ___ W idth (km) Distance [km) 0.0 0.5 1.5 2.5 3.5 4.5 5.0 0.577H 1.04 1.62 2.10 1.0 1.0 1.0 1.0 0.90 0.85 0.99 0.99 0.85 0.95 0.82 0.80 0.42(0.40) 0.90(0.89) 0.70(0.60) 0.80(0.69) 0.30 0.45 0.65 0.75 0.20 0.35 0.45 0.55 0.18(0.14) 0.20(0.15) 0.35(0.23) 0.50(0.34) x 1 0 ,-5 T h e m a g m a c h a m b e r w i t h h e i g h t 5 0 0 m a t d e p t h 2 . 0 k m 0 . 5 k m S ’ . O k m Figure 4.5: Vertical surface ground velocity for the case of a horizontally elongated magma chamber at depth 2.0 km with height 500m, width 2.0 km and rise time 0.3 sec. This case is used to compare with the case of vertically elongated magma chamber. 75 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 o"5 The magma chamber — height 400m width 0.577 km 0 . 0 0 . 5 k m 1 . 5 0 1 2 3 5 4 10 6 7 8 9 Figure 4.6: Vertical surface ground velocity for the case of a magma chamber height 400m and width 0.577km with rise time 0.3 sec 76 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s x t o '5 The magma chamber — height 400m width 1.04 km -2 -2 -2 -2 -2 4 . 5 -2 5 . 0 k m -2 0 1 2 3 5 6 7 4 8 9 10 Figure 4.7: Vertical surface ground velocity for the case of a magma chamber height 400m and width 1.04km with rise time 0.3 sec 77 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 Q -5 The magma chamber — height 400m width 1.62 km t - 2 -2 -2 -2 -2 0 1 2 3 5 4 6 7 8 10 9 Figure 4.8: Vertical surface ground velocity for the case of a magma chamber height 400m and width 1.62km with rise time 0.3 sec 78 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 o'5 The magma chamber — height 400m width 2.1 km (*10 O . J L k m (*10 3 . 5 4 . 5 5 . 0 k m 2 1 3 5 4 6 7 8 9 10 Figure 4.9: Vertical surface ground velocity for the case of a magma chamber height 400m and width 2.1km with rise time 0.3 sec 79 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 4.4: The normalized vertical surface ground velocity amplitude distribution for the case of a horizontally elongated magma chamber at different depth. The brack eted values in bold are the converted 3D amplitudes normalized to the amplitude at epicentral distance 0.5 km.____________________________________ Depth (km) Distance (km) 0.0 0.5 1.5 2.5 3.5 4.5 5.0 1.0 1.0 0.65 0.59 0.40(0.39) 0.30 0.20 0.10(0.07) 2.0 1.0 0.94 0.79 0.84(0.72) 0.63 0.42 0.31(0.20) 2.5 1.0 0.94 0.78 0.71(0.64) 0.67 0.56 0.55(0.40) 3.0 1.0 0.99 0.98 0.89(0.79) 0.79 0.77 0.71(0.52) 4.0 1.0 1.06 0.94 0.99(0.85) 0.89 0.78 0.72(0.54) 4.3.2 Effect of the depth We now investigate the effect of the depth for the case of the horizontally elongated magma chamber. We assume that the dimension of the magma chamber is 400m high and 2.1 km wide. We vary the depth of the magma chamber from 1.0 km to 4.0 km. The simulation results (Table 4.4) show that the depth have significant effect on amplitude distribution. The effect for the horizontally elongated low velocity magma chamber, is different from that for the homogeneous half-space case. The converted 3D normalized amplitude at 5.0km increases with depth until the depth is about 3.0km. When the depth exceeds 3.0 km, the converted 3D normalized amplitude at 5.0 km stays constant at 0.52 - 0.54. Thus, as the depth increases, the distribution of normalized amplitude becomes more Battened. For depth 2.5 km, the converted 3D normalized amplitude is 0.64 at 2.5 km, 0.40 at 5.0 km compared to observation 0.66 ± 0.12 (CHR) at 2.5 km and 0.45 ±0.10 (NCR) and 0.53 ±0.15 (NTR) at 5.0 km. As shown in the preceding section, at depth 2.5 km, the arrival delay time is also very close to the observed. As mentioned earlier, this will put the magma chamber slightly below the earthquake swarm occurring before the eruption. 80 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s x 10' Horizontally elongated magma chamber (400m) at depth 1.0km 0 : 5 k m 5 . 0 k m Figure 4.10: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 1.0km with rise time 0.8 sec 81 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s A p p c m / s x 10 Horizontally elongated magma chamber (400m) at depth 2.0km 0 . 5 k m 5 . 0 k m Figure 4.11: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 2.0km with rise time 0.8 sec 82 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x io~5 The magma chamber at depth of 2.5 km < h 0 -2 1 . 5 -2 -2 (fJO 3 . 5 -2 #10 4 . 5 -2 5 . 0 k m -2 0 2 3 1 4 5 8 6 7 9 10 Figure 4.12: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 2.5km with rise time 0.8 sec 83 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Horizontally elongated magma chamber (400m) at depth 3.0km 0 . 5 k m # 1 0 3 , 5 0 1 2 3 7 10 4 5 6 8 9 Figure 4.13: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 3.0km with rise time 0.8 sec 84 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 10 Horizontally elongated magma chamber (400m) at depth 4.0km _..... —r -----—i ---- ----r---------t — — — r--------- r t----------1 —----— : /0.0 a y — - ........ - -----C ----1 ----------1 ----------1 _ J - 1 1 1 1 > 2 - I I I 5 . 0 k m Figure 4.14: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) at depth 4.0km with rise time 0.8 sec 85 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x i o '5 The magma chamber with height 100m at depth 2.0 km 0.0 10 *10 : 1 . 5 3 . 5 4 . 5 0 2 3 5 1 4 6 7 8 10 9 Figure 4.15: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.1 km) with rise time 0.3 sec 4.3.3 Effect of the height (thickness) From the results in previous sections, we fix the depth of the low velocity magma chamber at 2.5 km and the width of the magma chamber 2.1 km. V V e now change the height(thickness) of the magma chamber from 100m to 500m at an interval of 100m to study how the amplitude distribution and the dominant frequency of the long period events depend on this parameter. 86 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 o'5 The magma chamber with height 200 m at depth 2.0 km -2 < 1 .5 -2 -2 3 . 5 -2 -2 5 . 0 k m -2 2 0 1 3 5 4 6 7 8 9 10 Figure 4.16: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.2 km) with rise time 0.3 sec 87 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 o '5 The magma chamber with height 300 m at depth 2.0 km 0.0 | - 2 <*!£ 3 . 5 5 . 0 k m 0 1 2 5 3 4 6 7 8 10 9 Figure 4.17: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.6 km) with rise time 0.3 sec 88 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The magma chamber with height 200 m at depth 2.0 km t f ’0 1 O f - a k 0 . 5 k m 1 U I , 1 . 5 1 o f - a •S5 1 2 l - a 3 . 5 f* 10 1 O f - a 4 . 5 #10 w 1 2 f - a 5 . 0 k m 0 2 3 1 4 5 6 7 8 10 9 Figure 4.18: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.4 km) with rise time 0.3 sec 89 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 o * 5 T h e m a g m a c h a m b e r w i t h h e i g h t 5 0 0 m a t d e p t h 2 . 0 k m 0 . 5 k m *10 3 . 5 *10 4 . 5 *10 Figure 4.19: Vertical surface ground velocity for the case of a magma chamber (2.1 x 0.5 km) with rise time 0.3 sec 90 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The simulation results (Fig 4.15, 4.16, 4.17, 4.18, 4.19) indicate that the height(thickness) does not have significant effect either on the amplitude distribution or the dominant frequency. Only the case of 100 m appears to be an exception, showing different waveforms from the other cases. In order to examine what determines the dominant frequency of LP events, we first considered a source time function which has a peak in the spectrum, namely p = sin(2, Kft)e~2^t (4.6) where e~ i s to make duration of the wavelet short so that the analysis of the results is easier. The dominant frequency of source pressure will be approximately f Hz. We varied f to see how it affect the dominant frequency and velocity observed at the free surface. We use a line pressure source whose dominant frequency is 0.5, 1.0, 1.5, 2.0 or 3.0 Hz, and put it in the magma chamber. The height of the magma chamber is 500m. Inside the chamber, P and S wave velocity are l.lk m /s and 0.55 km/s respectively. P wave length ranges from 0.37 km to 2.1 km, and S wave length ranges from 0.18 km to 1.1 km. We calculated the spectra of the vertical surface velocity at each location for different dominant frequency of line pressure source (Figure 4.20). There is very little variation among locations, and the dominant frequency at each location is generally close to the dominant frequency of line pressure source. From the above simulations, we conclude that the height (thickness) of magma chamber has little effect on the dominant frequency of LP events, which is primarily determined by the dominant frequency of the line pressure source. We can adjust 91 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x 1 0 - ® 5 H z 0.0km x10-^OHz #10 a > * 2 o a . 0 § 2 Q . 5 , (t!0 O k m x 1 0 i 1 . 5 H z ! L u " ~ “ A r - ; j k L ™ ______________________ 5 : * T O 5 ° J f ^ Z 5 ;« j£ i 0 . 0 k m _____ i --------------j 1 Y ? 2t — A 1 . 5 k m 0 . 5 i l l 1 . 5 k m 1 | 1 . 5 k m J W . o J V C o - J u l . I t 10 c ^ (jc 10 * u*10 © I- 0 . 5 k m L A f c . # 1° I' 1£ J l 3 . 5 k m 1 U — ! 10 2' 1 °rtT T O ' 5 Y --------35km- 5 k m 5 *£!£ 2 . 5 k m s § 0 . 5 C L 0 1 § § 0 . 5 CL 0 L M 5 . 0 k m 3 . 5 k m 5 4l£ 10 I t 4 . 5 k m 2 A . 4 5 * " 2 f _ - d " ° 4 ? 1 > ' 0 ______________________ 5 ? l f 1 ° 4 . 5 k m .Air- I iLii 5 0 5 . 0 k m x io-SOH z X 1 Q - 9 . 0 H Z 0 41 2 0 1 1 0 . 5 i r ™ #10' 5 M 5km 10 f #10 5 H i t #10 • 5 I 3 . 5 k m - J l - ■ 0 1 1 0 . 5 10 L 5 k m °4t 2* 1 0 2( 1 0 5 W 10 J L k m 10 J j \ 5kftv J i | | U . 5 k m # 1 0 5 r 5 0 5 0 1 ' 0 . 5 0 1 « 0 . 5 0 10 5 : 4 1 5 k m 5 . 0 k m H z H z H z H z HZ Figure 4.20: Spectra of the vertical surface ground velocity for the case of a magma chamber with 500m height with different dominant frequencies of source pressure time-function. From left to right, the dominant frequencies of the line pressure are 0.5, 1.0, 1.5 2.0 and 3.0Hz. 92 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Rock Fluid Chamber Low Velocity magma chamber Gas Chamber Figure 4.21: Schematic diagram of the complementary model the pressure source to generate the LP events with the dominant frequency observed at Piton de la Fournaise. But such adjustment does not provide any insight into the mechanism of the long period events. In the next section, we will discuss a possible mechanism that may control the dominant frequency of LP events. 4.4 R efined m odel In order to study what control the dominant frequency of the LP events at Piton de la Fournaise, we shall combine our model with the model of Crosson and Bame [18]. Crosson and Bame [IS] proposed a spherical source model for the long period event ( called “low frequency volcanic earthquakes” in their paper). In their model, there is a spherical fluid magma chamber in which a small gas chamber is embedded 93 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 4.5: Parameters for the Crosson-Bame’s model, where Q = [m(n)/Re(ft) Property value Radius of the fluid chamber (m) Mass density of the fluid (gfcm3) Q of the fluid 200 2.7 10 at the center. Our refined model contains the Crosson-Bame7 s model inside the horizontally elongated low velocity magma chamber as shown in Fig 4.21. The low velocity magma chamber still has P wave velocity 1.1 km /s and S wave velocity 0.55 km /s ( same values as in all studied cases in this thesis), while the fluid magma has higher compression wave velocity and extremely low S wave velocity. We choose P wave velocity 2.5 km /s and S wave velocity 0.005km/s. At the inner surface of the gas chamber, pressure suddenly increases in a step function. This model is too expensive to apply our finite element method, because of the need for many small elements to accommodate very short wave length in the fluid magma. Hence, we approximate our model with the spherical model of Crosson and Bame, for which a compact analytic solution exists. We compute the frequency domain solution using the velocities defined above and additional paxameters shown in table 4.5. Our numerical simulation show that when the inner radius is 1.7m and the outer radius is 200m, the dominant frequency of the fundamental model is 1.0 Hz (Fig 4.22); when the inner radius is 0.8 m and the outer radius is 200m, the dominant frequency is 2.0 Hz (Fig 4.23). The dominant frequency increases with the decreas ing inner radius as shown in Fig 4.24. 94 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Am plitude spectrum 10 .6 10' .5 1 0 - A 10 .3 10" 0 5 10 1 5 20 2 5 3 0 3 5 4 0 4 5 5 0 Frequency (Hz) Figure 4.22: Amplitude spectrum for the case of inner radius of 1.7m and outer radius of 200m, the dominant frequency is 1.0Hz 95 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 7 5 © Q . to © T9 3 ; 10 .6 10' .5 1 0 - ,4 10 0 5 10 1 5 2 5 20 3 0 3 5 4 0 4 5 5 0 Frequency (Hz) Figure 4.23: Amplitude spectrum for the case of inner radius of 0.8m and outer radius of 200m, the dominant frequency is 2.0 Hz 96 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 52 As mentioned in Chapter 2, observations at Piton de la Fournaise indicates that the 2 Hz LP generator may be deeper than 1 Hz LP generator. This may be at tributed to smaller gas chamber at deeper magma reservoir. Assuming the gas inside the gas chamber is ideal gas, the pressure P and volume V are related by PVjnT = constant (4.7) where n is molar quantity and T is temperature. Assuming the temperature of the gas is rather constant because it is inside the melt magma, the gas volume will decrease with depth because of increasing pressure for the same molar quantity. Fig 4.24 shows relation between the depth and the dominant frequency of fundamental mode under this assumption. Although the result is not supporting observation at La Reunion quantitatively ( 2Hz LP event cannot be as deep as 10 km), it gives the range of the parameters if the model is applicable to LP events at Piton de la Fournaise. 4.5 D iscussion and Sum m ary In section 2.3, we summarized key observations on LP events at Piton de la Fournaise in 7 items. In our finite element simulation, the item [5] about the symmetry of the LP event source with respect to the summit stations was the starting point. We modeled the source first as a line pressure source in a homogeneous half-space, and found that it cannot explain the observed amplitude distribution (item [6]) and the observed arrival time (item [7]) simultaneously. 97 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Relation between the inner radius and predominate frequency 2.5 £ 1 . 5 0 . 5 I n n e r r a d i u s R e l a t i o n b e t w e e n t h e d e p t h a n d p r e d o m i n a t e f r e q u e n c y 2 . 5 £ 1 - 5 0 . 5 D e p t h Figure 4.24: The relation between the inner radius (m) and the dominant frequency and inferred relation between the depth (km) and the dominant frequency 98 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. We then improved our model by placing the line pressure source in a low velocity magma chamber. We tried first vertically elongated magma chamber, simulating a dike, a common form of magma conduit as described in Chapter 1. We found that this model makes the discrepancy from the observation even greater than the case of a line pressure source in a homogeneous half-space. We found that horizontally elongated magma chamber ( geological term is a sill as described in Chapter 1) explains well observed amplitude distribution as well as arrival times. The best model is characterized by a chamber at the depth of 2.5 km and width of 2 km. This depth is close to the one inferred from comparison with the band-pass filtered tectonic events ( item[3]). Our conclusion that the magma chamber is elongated horizontally rather than vertically is in harmony with our observation that LP events are related with horizontal movement of magma (item[l]). Finally we suggest that the dominant frequency of LP events may be explained by a model of gas chamber embedded in a magma chamber as proposed by Crosson and Bame [18]. This mechanism explains, at least quantitatively, the observation that the 2Hz LP source is located below the 1Hz LP source (item[2]). Thus, our numerical simulation by the finite element method explained key ob servations described in Chapter 2, and gave new information about the geometry and depth of the source of LP events at Piton de la Fournaise. In addition to the information about the source, our simulation clarified the type of seismic waves dominating in different cases. P waves dominate in a homoge neous half-space with a line pressure source, but the presence of low velocity magma chamber generates trapped waves in and near the chamber as well as converted S 99 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. waves outside the chamber. We found that S waves become dominant for the case of horizontally elongated magma chamber at distant stations. Although the simulation presented in this thesis is only 2D, our program is capable of 3D modeling. We are cautious about 3D modeling, however, because without a thorough study of 2D problems and guidance from 2D investigation, 3D modeling does not necessarily produce more meaningful results. This is exactly our goal here to investigate the problems in 2D to have better understanding of the seemingly simple but actually very complex problem before we study the long period events in 3D. We also emphasize the methodology on modeling. In this thesis work, we promote object-oriented scientific computing to facilitate the collaboration in the seismology community, which already widely practiced in industry and other academic commu nities and really lacks in our community. We have collaborated by publishing papers in the past, but we have not done that with regard to tools. This kind of collabora tion will greatly benefit seismologist in many ways. We recommend Diffpack for the free collaborative framework, which is developed at Oslo University using object- oriented language C + + . We believe that object-oriented scientific computing and collaboration are the future trends for the seismologica! numerical simulation. 100 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A p p en d ix A F in ite elem ent procedures, O bject-oriented concepts and CH— |- header files A .l Basic P rocedures in F inite E lem ent M ethod P re-processing Finite element analysis is comprised of pre-processing, solution and post-processing phases. The goals of pre-processing are to develop an appro priate finite element mesh, assign suitable material properties, and apply boundary conditions in the form of restraints and loads. The finite element mesh subdivides the geometry into elements, upon which are found nodes. The nodes, which are really just point locations in space, are generally located at the element corners and perhaps near each midside. For a two-dimensional (2D) analysis, or a three-dimensional (3D) thin shell analysis, the elements are es sentially 2D, but may be “warped” slightly to conform to a 3D surface. An example is the thin shell linear quadrilateral; thin shell implies essentially classical shell the ory, linear defines the interpolation of mathematical quantities across the element, and quadrilateral describes the geometry. For a 3D solid analysis, the elements have 101 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. physical thickness in ail three dimensions. Common examples include solid linear brick and solid parabolic tetrahedral elements. In addition, there are many spe cial elements, such as axisymmetric elements for situations in which the geometry, material and boundary conditions are all symmetric about an axis. The model’s degrees of freedom (DOF) are assigned at the nodes. Solid elements generally have three translational dof per node. Rotations are accomplished through translations of groups of nodes relative to other nodes. Thin shell elements, on the other hand, have six dof per node: three translations and three rotations. The addition of rotational dof allows for evaluation of quantities through the shell, such as bending stresses due to rotation of one node relative to another. Thus, for structures in which classical thin shell theory is a valid approximation, carrying extra dof at each node bypasses the necessity of modeling the physical thickness. The assignment of nodal dof also depends on the class of analysis. For a thermal analysis, for example, only one temperature dof exists at each node. The geometry is meshed with a mapping algorithm or an automatic free-meshing algorithm. The first maps a rectangular grid onto a geometric region, which must therefore have the correct number of sides. Mapped meshes can use the accurate and cheap solid linear brick 3D element, but can be very time-consuming, if not impossible, to apply to complex geometries. Free-meshing automatically subdivides meshing regions into elements, with the advantages of fast meshing, easy mesh-size transitioning (for a denser mesh in regions of large gradient), and adaptive capa bilities. Disadvantages include generation of huge models, generation of distorted elements, and, in 3D, the use of the rather expensive solid parabolic tetrahedral element. It is always important to check elemental distortion prior to solution. A 102 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. badly distorted element will cause a matrix singularity, killing the solution. A less distorted element may solve, but can deliver very poor answers. Acceptable levels of distortion are dependent upon the solver being used. In simple cases (homogenous and half-space homogenous), material properties can be easily assigned without providing complicated interfaces (API). To handle hetergenous cases (e.g. material properties vary with depths) , it is better to separate the material properties handling from mesh generation. Solution While the pre-processing and post-processing phases of the finite element method are interactive and time-consuming for the analyst, the solution is usually a batch process, and is demanding of computer resource. The governing equations are assembled into matrix form and are solved numerically. The assembly process depends not only on the type of analysis (e.g. static or dynamic), but also on the model’s element types and properties, material properties and boundary conditions. In the case of a linear static structural analysis, the assembled equation is of the form Kx = 6, where K is the system stiffness matrix, x is the nodal degree of freedom (dof) displacement vector, and b is the applied nodal load vector. To appreciate this equation, one must begin with the underlying elasticity theory. The strain-displacement relation may be introduced into the stress-strain relation to ex press stress in terms of displacement. Under the assumption of compatibility, the differential equations of equilibrium in concert with the boundary conditions then determine a unique displacement field solution, which in turn determines the strain and stress fields. The chances of directly solving these equations are slim to none for 103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. anything but the most trivial geometries, hence the need for approximate numerical techniques presents itself. A finite element mesh is actually a displacement-nodal displacement relation, which, through the element interpolation scheme, determines the displacement any where in an element given the values of its nodal dof. Introducing this relation into the strain-displacement relation, we may express strain in terms of the nodal dis placement, element interpolation scheme and differential operator matrix. Recalling that the expression for the potential energy of an elastic body includes an integral for strain energy stored (dependent upon the strain field) and integrals for work done by external forces (dependent upon the displacement field), we can therefore express system potential energy in terms of nodal displacement. Applying the principle of minimum potential energy, we may set the partial derivative of potential energy with respect to the nodal dof vector to zero, resulting in: a summation of element stiffness integrals, multiplied by the nodal displacement vector, equals a summation of load integrals. Each stiffness integral results in an element stiffness matrix, which sum to produce the system stiffness matrix, and the summation of load integrals yields the applied load vector, resulting in Kx = 6. In practice, integration rules are applied to elements, loads appear in the b vector, and nodal dof boundary conditions may appear in the x vector or may be partitioned out of the equation. Solution methods for finite element matrix equations are plentiful. In the case of the linear static K x = 6, inverting K is computationally expensive and numerically unstable. A better technique is Cholesky factorization, a form of Gauss elimination, and a minor variation on the the “LDU” factorization theme. The K m atrix may 104 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. be efficiently factorized into LDU, where L is lower triangular, D is diagonal, and U is upper triangular, resulting in LDUx = b. Since L and D are easily inverted, and U is upper triangular, x may be determined by back-substitution. Another popular approach is the wavefront method, which assembles and reduces the equations at the same time. Some of the best modern solution methods employ sparse matrix techniques. Because node-to-node stiffnesses are non-zero only for nearby node pairs, the stiffness matrix has a large number of zero entries ( Figure A .l). This can be exploited to reduce solution time and storage by a factor of 10 or more. By taking advantage of the sparsity of stiffness matrix, more efficient iterative and multigrid solvers are developed, which are discussed in Chapter 4. Post-processing After a finite element model has been prepared and checked, boundary conditions have been applied, and the model has been solved, it is time to investigate the results of the analysis. This activity is known as the post-processing phase of the finite element method. Post-processing begins with a thorough check for problems that may have oc curred during solution. Most solvers provide a log file, which should be searched for warnings or errors, and which will also provide a quantitative measure of how well-behaved the numerical procedures were during solution. Next, reaction loads at restrained nodes should be summed and examined as a “sanity check”. Reaction loads that do not closely balance the applied load resultant for a linear static anal ysis should cast doubt on the validity of other results. Error norms such as strain energy density and stress deviation among adjacent elements might be looked at 105 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Y-Axis Tbn O ct 29 1J J « :! 0 19M Sparsity pattern 8 4 7 8 0 0 - 7 0 0 - 6 0 0 5 0 0 - 4 0 0 - 3 0 0 - 200- 3 100-3 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 X-Axis Figure A.l: Sparse Pattern of the Stiffness m atrix of a 3D model 106 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. next, but for h-code analyses these quantities are best used to target subsequent adaptive remeshing. Once the solution is verified to be free of numerical problems, the quantities of interest may be examined. Many display options axe available, the choice of which depends on the mathematical form of the quantity as well as its physical meaning. For example, the displacement of a solid linear brick element’s node is a 3- component spatial vector, and the model’s overall displacement is often displayed by superposing the deformed shape over the undeformed shape. Dynamic viewing and animation capabilities aid greatly in obtaining an understanding of the deformation pattern. Stresses, being tensor quantities, currently lack a good single visualization technique, and thus derived stress quantities are extracted and displayed. Principle stress vectors may be displayed as color-coded arrows, indicating both direction and magnitude. The magnitude of principle stresses or of a scalar failure stress such as the Von Mises stress may be displayed on the model as colored bands. When this type of display is treated as a 3D object subjected to light sources, the resulting image is known as a shaded image stress plot. Displacement magnitude may also be displayed by colored bands, but this can lead to misinterpretation as a stress plot. An area of post-processing that is rapidly gaining popularity is that of adaptive remeshing. Error norms such as strain energy density are used to remesh the model, placing a denser mesh in regions needing improvement and a coarser mesh in areas of overkill. Adaptive remeshing is a recent demonstration of the iterative nature of h-code analysis. Optimization is another area enjoying recent advancement. Based on the val ues of various results, the model is modified automatically in an attem pt to satisfy 107 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. certain performance criteria and is solved again. The process iterates until some convergence criterion is met. In its scalar form, optimization modifies beam cross- sectional properties, thin shell thicknesses and/or material properties in an attem pt to meet maximum stress constraints, maximum deflection constraints, and/or vibra tional frequency constraints. Shape optimization is more complex, with the actual 3D model boundaries being modified. This is best accomplished by using the driving dimensions as optimization parameters, but mesh quality at each iteration can be a concern. A .2 Basic concepts in the O bject M odel O bject is a software bundle of variables and related methods. As the name object- oriented implies, objects are key to understanding object-oriented technology. You can look around you now and see many examples of real-world objects: your dog, your desk, your television set, your bicycle. These real-world objects share two char acteristics: they all have state and they all have behavior. For example, dogs have state (name, color, breed, hungry) and dogs have behavior (barking, fetching, and slobbering on your newly cleaned slacks). Bicycles have state (current gear, current pedal cadence, two wheels, number of gears) and behavior (braking, accelerating, slowing down, changing gears). Software objects are modeled after real-world objects in that they, too, have state and behavior. A software object maintains its state in variables and implements its behavior with methods. 108 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. You can represent real-world objects using software objects. You might want to represent real-world dogs as software objects in an animation program or a real- world bicycle as a software object within an electronic exercise bike. However, you can also use software objects to model abstract concepts. For example, an event is a common object used in GUI window systems to represent the action of a user pressing a mouse button or a key on the keyboard. Everything that the software object knows (state) and can do (behavior) is ex pressed by the variables and methods within that object. A software object that modelled your real-world bicycle would have variables that indicated the bicycle’s current state: its speed is 10 mph, its pedal cadence is 90 rpm, and its current gear is the 5th gear. These variables and methods are formally known as instance variables and instance methods to distinguish them from class variables and class methods (described later in What Are Classes?). The software bicycle would also have methods to brake, change the pedal cadence, and change gears. (The bike would not have a method for changing the speed of the bicycle, as the bike’s speed is really just a side-effect of what gear it’s in, how fast the rider is pedaling, whether the brakes are on, and how steep the hill is.) Anything that an object does not know or cannot do is excluded from the object. For example, your bicycle (probably) doesn’t have a name, and it can’t run, bark, or fetch. Thus there are no variables or methods for those states and behaviors in the bicycle class. The object’s variables make up the center or nucleus of the object. Methods sur round and hide the object’s nucleus from other objects in the program. Packaging 109 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. an object’s variables within the protective custody of its methods is called encapsu lation. Typically, encapsulation is used to hide unimportant implementation details from other objects. When you want to change gears on your bicycle, you don’t need to know how the gear mechanism works, you just need to know which lever to move. Similarly in software programs, you don’t need to know how a class is implemented, you just need to know which methods to invoke. Thus, the implementation details can change at any time without affecting other parts of the program. This conceptual picture of an object-a nucleus of variables packaged within a protective membrane of methods-is an ideal representation of an object and is the ideal that designers of object-oriented systems strive for. However, it’s not the whole story. Often, for implementation or efficiency reasons, an object may wish to expose some of its variables or hide some of its methods. In many languages, an object can choose to expose its variables to other objects allowing those other objects to inspect and even modify the variables. Also, an object can choose to hide methods from other objects forbidding those objects from invoking the methods. An object has complete control over whether other objects can access its variables and methods and in fact, can specify which other objects have access. Variable and method access in Java is covered in Controlling Access to Members of a Class. Encapsulating related variables and methods into a neat software bundle is a simple yet powerful idea that provides two primary benefits to software developers: • Modularity-The source code for an object can be written and maintained independently of the source code for other objects. Also, an object can be 110 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. easily passed around in the system. You can give your bicycle to someone else and it will still work. • Information hiding-An object has a public interface that other objects can use to communicate with it. But the object can maintain private information and methods that can be changed at any time without affecting the other objects that depend on it. You don’t need to understand the gear mechanism on your bike in order to use it. C lass is a blueprint or prototype that defines the variables and methods common to all objects of a certain kind. In the real world, you often have many objects of the same kind. For example, your bicycle is just one of many bicycles in the world. Using object-oriented terminology, we say that your bicycle object is an instance of the class of objects known as bicycles. Bicycles have some state (current gear, current cadence, two wheels) and behavior (change gears, brake) in common. However, each bicycle’s state is independent of and can be different from other bicycles. When building bicycles, manufacturers take advantage of the fact that bicycles share characteristics by building many bicycles from the same blueprint-it would be very inefficient to produce a new blueprint for every individual bicycle they manufactured. In object-oriented software, it’s also possible to have many objects of the same kind that share characteristics: rectangles, employee records, video clips and so on. Like the bicycle manufacturers, you can take advantage of the fact that objects of the same kind are similar and you can create a blueprint for those objects. Software “blueprints” for objects are called classes. I l l Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For example, you could create a bicycle class that declares several instance vari ables to contain the current geax, the current cadence, and so on, for each bicycle object. The class would also declare and provide implementations for the instance methods that allow the rider to change gears, brake, and change the pedaling ca dence. The values for instance variables are provided by each instance of the class. So, after you’ve created the bicycle class, you must instantiate it (create an instance of it) before you can use it. When you create an instance of a class, you create an object of that type and the system allocates memory for the instance variables declared by the class. Then you can invoke the object’s instance methods to make it do some thing. Instances of the same class share the same instance method implementations (method implementations are not duplicated on a per object basis), which reside in the class itself. In addition to instance variables and methods, classes can also define class variables and class methods. You can access class variables and methods from an instance of the class or directly from a class-you don’t have to instantiate a class to use its class variables and methods. Class methods can only operate on class variables-they do not have access to instance variables or instance methods. The system creates a single copy of all class variables for a class the first time it encounters the class in a program-all instances of that class share its class variables. For example, suppose that all bicycles had the same number of gears. In this case defining an instance variable for number of gears is inefficient-each instance would have its own copy of the variable, but the value would be the same for every instance. In situations such as this, you could define a class variable that contains the number 112 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. of gears. All instances share this variable. If one object changes the variable, it changes for all other objects of that type. Instance and Class Members discusses instance variables and methods and class variables and methods in detail. You probably noticed that the illustrations of objects and classes look very similar to one another. And indeed, the difference between classes and objects is often the source of some confusion. In the real world it’s obvious that classes are not themselves the objects that they describe-a blueprint of a bicycle is not a bicycle. However, it’s a little more difficult to differentiate classes and objects in software. This is partially because software objects are merely electronic models of real-world objects or abstract concepts in the first place. But it’s also because many people use the term ^object” inconsistently and use it to refer to both classes and instances. Objects provide the benefit of modularity and information hiding. Classes pro vide the benefit of reusability. Bicycle manufacturers reuse the same blueprint over and over again to build lots of bicycles. Software programmers use the same class, and thus the same code, over and over again to create many objects. In h e rita n c e Generally speaking, objects are defined in terms of classes. You know a lot about an object by knowing its class. Even if you don’t know what a penny- farthing is, if I told you it was a bicycle, you would know that it had two wheels, handle bars, and pedals. Object-oriented systems take this a step further and allow classes to be defined in terms of other classes. For example, mountain bikes, racing bikes, and tandems are all different kinds of bicycles. In object-oriented terminology, mountain bikes, racing bikes, and tandems are all subclasses of the bicycle class. 113 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Similarly, the bicycle class is the superclass of mountain bikes, racing bikes, and tandems. Each subclass inherits state (in the form of variable declarations) from the super class. Mountain bikes, racing bikes, and tandems share some states: cadence, speed, and the like. Also, each subclass inherits methods from the superclass. Moun tain bikes, racing bikes, and tandems share some behaviors: braking and changing pedaling speed, for example. However, subclasses are not limited to the state and behaviors provided to them by their superclass. What would be the point in that? Subclasses can add variables and methods to the ones they inherit from the super class. Tandem bicycles have two seats and two sets of handle bars; some mountain bikes have an extra set of gears with a lower gear ratio. Subclasses can also override inherited methods and provide specialized imple mentations for those methods. For example, if you had a mountain bike with an extra set of gears, you would override the “change gears” method so that the rider could actually use those new gears. You axe not limited to just one layer of inheritance. The inheritance tree, or class hierarchy, can be as deep as needed. Methods and variables are inherited down through the levels. In general, the further down in the hierarchy a class appears, the more specialized its behavior. Subclasses provide specialized behaviors from the basis of common elements pro vided by the superclass. Through the use of inheritance, programmers can reuse the code in the superclass many times. Programmers can implement superclasses called abstract classes that define “generic” behaviors. The abstract superclass defines 114 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. and may partially implement the behavior but much of the class is undefined and unimplemented. Other programmers fill in the details with specialized subclasses. A .3 C +-b System H eader files In this section, we list some of the major header files to illustrate the object-oriented approach used in the system developed. A .3.1 ElasticWave / / C++ /* Jinbo Chen’s ElasticWave */ / * * / /* (c) Copyright 1998 */ /• Department of Earth Sciences */ /* Univeristy of Southern California * / /* All rights reserved. */ / * * / / • * / //#indude<Copyright. h> tifndef ElasticWave.h.IS.INCLUDED tdefine ElasticWave_h_IS_INCLUDED ♦include <MenuUDC.h> ♦include <FEM.h> ♦include <LinEqAdm.h> ♦include <TimePrm.h> ♦include <GridFE.h> ♦include <DegFreeF£.h> ♦include <Arrays.real.h> ♦include <Store4Plotting.h> ♦include <FieldFormat.h> //= = ===========«=,-■!------------ ^rr=== === = s== = s= = ==s= // Design Goal : Flexibility and Modularity, Efficiency // - Implmentation should be not specific to any grid // ElasticMaterialProperty. / / = = = = = = = = = = = = = = = * = = = = = = = = = = = = « = = = = = = = = = = „ — = = = = — //ToDo: ElasticMaterialProperty, ModelPack, BoundauryConditions / / // ElasticMaterialProperty: // -represent the material matrix for 2D in Plain Strain or 3D // -also make it easy to construct the stiffness matrix 115 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. // ModelPack: // - generate grid based on PreproBox, PreproSupElSet, GeomPack // - generate duplicate node, crack node if needed. // - associate the boundary Indicators with the actual boundary // condition Impl, either essential or natural // Boundary Condition: // - act as a collection of boundary conditions. // Report Pack; // - dump the results class ElasticVave : public FEM, public MenuUDC, public Store4Plotting / / ----------------------------------------------------------------------------------------------------------- < //to create Nesh vith explosion source and crack source // ModelPack have possible plugin to generate grid / / friend class ElasticMaterialProperty; friend class ModelPack; friend class BoundaryConditions; friend class ReportPack; protected: / / ==== = = = ======r=r -------------- = ==s==^ =s=^ // Grid / /=================~ ......... Handle(GridFE) grid; // Mesh generated by special preprocessor / / ================ ===================== ==== // Displacement Field and Solution / / ========================T====J============ Handle(FieldsFE) u; // u at time step n+1 Handle(FieldsFE) u.prev; // displacement at time step n Handle(FieldsFE) v.prev; // velocity at time step n Handle(FieldsFE) a.prev; // acceleration at time step n //===a==== = =s===s= =s===== ==== =s= s==s====s= = ==ss=^ //Degree of Freedom hander and Linear Equation Administrator / / ===g= = s ==== s===s =r ======== ==, — Handle(DegFreeFE) dof; // mapping: field <-> equation system Handle(LinEqAdm) lineq // Linear Equation Manager //======== = ==== = = = : ========== = = •= === ~ ^ = ~ = // Time Parameter //===S========S==SJ==============a=s===s=s=s= TimePrm tip; // time interval, time step length etc /! 116 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. // Material Properties: Isotropic or Non-isotropic / / = = = ^ - . . s==B= = a ^ rT a r. - ^ 1 .T=__TJ Handle(ElasticMaterialProperty) property; / / = = = = «= = = « « ,= , ~ " - Z ~ T // Crack Nodes Condition Handler / / « = = = « ------------------------8 - T .T.™ — M - i M . ----------------------■— I L - T // Placeholder here Handle(BoundaryConditions) conditions; / / = = = S= = = == = = = t : r ------------ // IntegrationSchema / / = = = = = = = = = = = " S S S S S S S S S S B T i S B ----------------------------------------- Handle(IntegrationSchema) schema; public: ElasticWave (); virtual 'ElasticWave (); //================= menu system =============== ===== virtual void adm(MenuSystemk menu); static void defineStatic (MenuSystem* menu, int level = MAIN); virtual void define (MenuSystem* menu, int level = MAIN) - C defineStatic (menu, level); > virtual void scan (MenuSystem* menu) ; //we inherit the standard version of virtual void adm from class MenuUDC // ============== overide FEM: :makeSystem() ====== // Purpose: since the Left Matrix is time independent, Right Side Vector // can aslo exploit this characteristics to make it more // efficient, virtual void makeMassMatrix (GridFE* .grid, Matrix(NUMT)* mm, BooLean lumped); virtual void makeSystem(DegFreeFE* _dof, LinEqAdm* _lineq, BooLean compute.*, BooLean compute_Rhs); // similar to makeMassMatrix, since there are case, damping // coefficients are different along x, y, z direction. // this is useful to study wave splitting virtual void makeOampingMatrix(GridFE* .grid, Matrix(NUMT)* dm); // =======-====-= Time Step handling ============== virtual void solveProblem (); virtual void resultReport (); protected: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. virtual BooLean ok () const; virtual void init(); // set Initial condition for this problem void setIC (); //set the Essential Boundary conditions, void fillEssBC (); void timeLoop (); void solveAtThisTimeStep (); virtaul void CalcElmMatVec(int e, ElmMatVeck elmat, FiniteElementk fe); //volume integral virtual void integrands (ElmMatVeck elmat, FiniteElementk fe); //side integral virtual void integrands4side(int side, int boind, ElmMatVeck elmat, FiniteElementk fe); void saveResults (); void updateDataStructures (); >; /*>ElasticWave:*/ /*Class:ElasticWave NAME: ElasticWave - tem) SYNTAX: «ElmOef KEYWORDS: finite elements, ModelPack, ElasticMaterialProperty, BoundaryConditions DESCRIPTION: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. NOTE: Some comments in the documentation of this class refer to mixed finite element methods. Support for such methods is only to a very CONSTRUCTORS AND INITIALIZATION: This base class is initialized completely by its only constructor (which is protected and only called from the subclass objects). MEMBER FUNCTIONS: "refill" - on basis of the element type (and the number of space dimensions EXAMPLES: SEE ALSO: DEVELOPED BY: Department of Earth Sciences University of Southern California AUTHOR: Jinbo Chen End: ♦ / // classes for mass matrix /*<MassMatIntg:♦/ class VaveMassMatlntg : public IntegrandCalc BooLean lumped; ElasticMaterialProperty property; public: VaveMassMatlntg (BooLean lumped. = dpTRUE, ElasticMaterialProperty* _p) ■ C this->lumped = lumped.; this->property=_p; >; "VaveMassMatlntg () {} virtual void integrands (ElmMatVeck elmat, FiniteElementk fe); >; /♦>MassMatIntg:♦/ /♦Class:MassMatlntg NAME: MassMatlntg - an IntegrandCalc functor for mass matrices Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. SYNTAX: CMas sMatIntg KEYWORDS: finite element method, mass matrix DESCRIPTION: The class offers the integrand of the Beak form for the mass matrix in finite element computations. It is used by the function "FEM::makeMassMatrix" (which makes it easy for a programmer to compute a mass matrix). Besides that, it has probably few applications. SEE ALSO: class "FEM" DEVELOPED BY: SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo, Dept, of Mathematics, Norway AUTHOR: Hans Petter Langtangen, SINTEF/UiO End: ♦ / // classes for damp matrix /*<MassMatIntg:*/ class VaveDampMatlntg : public IntegrandCalc ElasticMaterialProperty property; public: VaveDampMatlntg (ElasticMaterialPropertyk _p) { this->property = _p;> 'VaveDampMatlntg () {} virtual void integrands (ElmMatVeck elmat, FiniteElementk fe); >; /*>MassMatIntg:*/ /♦Class:MassMatIntg Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. NAME: MassMatlntg - an IntegrandCalc functor for mass matrices SYNTAX: QMassMatlntg KEYWORDS: finite element method, mass matrix DESCRIPTION: The class offers the integrand of the weak form for the mass matrix in finite element computations. It is used by the function "FEM:: makeMassMatrix" (which makes it easy for a programmer to compute a mass matrix). Besides that, it has probably feu applications. SEE ALSO: class "FEM” DEVELOPED BY: SINTEF Applied Mathematics, Oslo, Norway, and University of Oslo, Dept, of Mathematics, Norway AUTHOR: Hams Petter Langtangen, SINTEF/UiO End: * / tendif A.3.2 Numerical Integration // C++ /* Jinbo Chen’s ElasticWave * / * * /* (c) Copyright 1998 * /* Department of Earth Sciences * / * Univeristy of Southern California * /* All rights reserved. * /* * /* * 121 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. /« //#include<Copyright.h> ♦ifndef Nevmark_h_IS_INCLUDED ♦define Newmark_h_IS_INCLUDED ♦include <MenuUDC.h> ♦include <FEM.h> ♦include <LinEqAdm.h> ♦include <TimePrm.h> ♦include <GridFE.h> ♦include <DegFreeFE.h> ♦include <Arraya.real.h> ♦include <Store4Plotting.b> ♦include <FieldFormat. h> // classes for NevMark integration scheme class IntegrationSchema : public Handleld protected: Vec(real) Vec(real) Vec(real) Vec(real) Vec(real) tmpl; tmp2; scratch_u; scratch_v; scracth.a; // obtained from lineq // vector for velocity // vector for acceleration DegFreeFE TimePrm dof; tip; Handle(Matrix(real)) mm; Handle(MatrixCreal)) dc; Handle(VectorCreal)) vb; public: IntegrationSchemaCDegFreeFe* _dof, TimePrm k tip); *IntegrationSchema(); BooLean ok(); void initHatricsCLineEqAdm k lineq); void sync(FieldsFE* u, FieldsFE* v. FieldsFE* a); void setHassHatrizCMatrix(real) * _mm) { mm.rebind(_mm);> void setDampHatrix(Matrix(real) * _dc) { dc.rebind(.dc);> void copyFields2Vec(FieldsFE* u, FieldsFE* v, FieldsFE* a); void copyVec2Fields(FieldsFE* u, FieldsFE* v, FieldsFE* a); Matrix(NUMT)* M(); Matrix(HOMT)* C(); Vector(NUMT)* B(); virtual void calcEffectiveLoad (LineEqAdm * lineq); 122 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. virtual void update(Vector f t .solution, FieldsFEft u, FieldsFEft v, FieldsFEft a); virtual void update(Vector f t .solution); virtual void eompute£ffectiveStiffnessNatrix(Nat ftA, Nat f t eK, Nat fteH,Nat f t eC); protected: virtual void setlntegrationCoefficientsO; void init(DegFreeFEft dof, TimePrm f t tip); > class Nevmark : public IntegrationSchema protected: VecSimple(real) ca(8); //coefficients based on delta T real alpha; real delta; public: Nevmark CrealDegFreeFE f t _dof, TimePrmft .tip, real alpha = 0.25, reed delta=0.5 ); 'Nevmark (); virtual void calcEffectiveLoad (LineEqAdm f t lineq); virtual void update(Vector f t .solution, FieldsFEft u, FieldsFEft v, FieldsFEft a); virtual void update(Vector f t .solution); virtual void computeEffectiveStiffnessHatriz(Hat ftA, Nat f t eK, Nat fteH,Hat f t eC); >; ftdefine ClassType IntegrationSchema ^include <Handle.h> *undef ClassType tendif A .3.3 Elastic Material Handling / / C++ Jinbo Chen’s Blastictfave (c) Copyright 1998 Department of Earth Sciences Univeristy of Southern California All rights reserved. > * / */ */ */ */ * / */ * / 123 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. /* /* */ ' * / //#include<Copyright.h> #ifndef IsotropicMaterialProperty_h_IS_INCLUDED #def ine IsotropicMaterialProperty_h_IS_INCLUDED / / = = == = = = --==== = = = == ---------------- // classes to get Material property / / ------------------------------------------------------------------------------------------ class IsotropicMaterialProperty : public ElasticMaterialProperty protected: / / = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = // one way to specify the elastic material property is to // use wave velocity and density // it is intuitive to the seismologists. / / = = = = = = = = = = = = = = = = = = = = = = = == = = = = = = = = = = = = // m U* * + cU’ + k U =f / / / / / / = = = = = = = = = = = = = = = = = = = = = = = = = = =- = = = = = = = = == = = = = = = real pVelocity; real sVelocity; public: IsotropicMaterialProperty(); "IsotropicMaterialProperty(); / / = = = = = = = = = = = = = = = = = = = == = = = = = = = == = = = = = = // virtual function: allow to overwrite on sub class / / = = = = = = = = = = = == === ==== = = = ■ ------------------------------- virtual Hat* getMaterialMatrizO; virtual void getMaterialMatrix(Mat t mat); / / = = = = = = == = = = = = = ======== = = ===== //operation + and *, / / = = ===== ===== = = = = = = == = = = = // result = this*a + adder*b virtual void add(real a=l.0, ElasticMaterialProperty k adder, real b=l.0, ElasticMaterialProperty f t result); virtual void scanCls is); 124 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. virtual void print(Os os); / / = = r -r - — -------------------, „ . ;1J= SS=S=S===S 5 =„ SS=;=S //only function: set elements within this domain of grid //to have the sub domain be the region.no .......... - ■ = s= a s s ==s=5, s= r ^ g = > #endif // C++ /* Jinbo Chen’s ElasticVave /* /* (c) Copyright 1998 /* Department of Earth Sciences /* Univeristy of Southern California /* *11 rights reserved. /* / * //#indude<Copyright .h> ftifndef AnisotropicMaterialProperty_h_IS_INCLUDED #def ine AnisotropicMaterialProperty_h_IS_INCLUDED / / = = = S==gS==SS := = = == = S :S=== ===ss===s==aB = =^ = s=s=s===s= // classes to get Material property / / ------------------------------------------------------------------------------------------ class AnisotropicMaterialProperty : public ElasticMaterialProperty { protected: public: AnisotropicMaterialProperty(); "AnisotropicMaterialProperty(); / / = = SS========.:=J. , == = = = = = « = = // virtual function: allow to overwrite on sub class //==== ====: — .-==============r===r==-.:====== Virtual Mat* getMaterialMatrixO; virtual void getMaterialMatrixCMat k mat); / / ==r«=„ r r ---------- - :sa= = = sss= ss=s=s« s5s===s= //operation + and *, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. / / ======= : ■ . ■ // result = this*a + adder*b virtual void addCreal a=l.0, ElasticMaterialProperty k adder, real b=1.0, ElasticMaterialProperty k result); / / = = = = = = r.- - - - - - - - - - - - - - - - - - - - -T-. . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . // MenuSystem / / = = = = = r ^ = = = = = s = s = = = r s = = = = = = = = = = s = r = ..................... virtual void scanCIs is); virtual void print(0s os); > •endif // C++ -*- /* Jinbo Chen's ElasticWave /* /* (c) Copyright 1998 /* Department of Earth Sciences /* Univeristy of Southern California /* All rights reserved. /* /* //#indude<Copyright. h> •ifndef ElasticMaterialHandler_h_IS_INCLUDED #define ElasticMaterialHandler_h_IS_INCLUDED *include <MenuUDC.h> •include <FEM.h> •include <LinEqAdm.h> •include <TimePrm.h> •include <GridFE.h> •include <DegFreeFE.h> •include <Arrays_real.h> •include <Store4Plott ing.h> •include <FieldFormat.h> class ElasticMaterialHandler : public Handleld, public MenuUDP protected: //GridFE grid; Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. //========= ====== „ ==as======ss=====sa: ==s=aB=„ s==- //An Array of Material Region / / = = = " = -- ............... :--------------- - ■ ■ -------------- r — r - int nregions; ArrayGenSimplest(MaterialDomain) regions; //== = = sa===»===r»=az r „ - .= , , ==s« = «==r.=====«s== sa== » a=r. //material property for regions not specified in above regions / / = = s= =s==== = = s===s= a=== =^======s===== ^=^^==^=^====== = ElasticMaterialProperty defaultValue; public: / / = = = J = = = == = ^ = = « = r = = SS = == = = S S a j = ^ z = = == ==^ = = = // constructor //= =====^s======s=g= a=s= === =J= J= ^==^^=====s==== = ElasticMaterialHandlerO; 'ElasticMaterialHandlerO; BooLean ok(); //= === ==== ========= === ==== = =-=~~~==~~==~ ^===== //MenuSystem / / = = = = ============ = ===«=== ===== == = = = = = static void def ineStatic (MenuSystem ftmenu, int level = MAIN); void define(MenuSystem fcmenu, int level = MAIN) { defineStatic(menu,level);} void scan(MenuSystemft menu); void adm(MenuSystem f t menu); //======;== = ===« =^=== = = = = = ^ = = := «=== = ====== // getter / / = === S=== ====== = = =S= S===== = SJ========S ===== real getElmOamping(int elm_no) const; real getElmMassOensity(int elm.no) const; void getMaterialMatriz(Ptv f t p, Matrix f t mat); real getDamping(Ptv f t p) const; real getMassOensity(Ptv ftp) const; int getMaterialOomain(Ptv ftp) const; protected: //usually get called by scan() void addMaterialRegion(MaterialRegionft region); void removeMaterialRegion(MaterialRegionft region); void removeMaterialRegion.(String regionld); > •define ClassType ElasticMaterialHandler •include <Handle.h> •undef ClassType Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ten d if // C++ / * Jinbo Chen's ElasticWave / * / * (c) Copyright 1998 /* Department of Earth Sciences /* Univeristy of Southern California /* All rights reserved. /* /* /t******************,***,**,,**,*,*****,,,**,,**,,,*,,*,**, //#include<Copyright.h> #ifndef Elast icHaterialProperty_h_IS_INCLUDED #def ine Elast icMaterialProperty_h_IS_IMCLODED #include <NenuUDC.h> tinclude <Store4Plotting.h> ^include <Complex.h> // classes to get Material property class ElasticMaterialProperty : public Handleld { protected: int nsd; //Dimension / / ====s===«»==r a = = :r = = = ======3= ====s= ==r= = s=s= // one way to specify the elastic material property is to // use wave velocity and density // it is intuitive to the seismologists. //========================== ■ // m 0" + cD' + k U =f / / / / / / = ==== ==== == == = = ==== = ======^ ------------------ real density; real damping; //===========.==: ^ - - ===== = = ==== - //Question ?: how to represent Q in terms of c // q = c/density; / / S===g===S=S=S===a===a= S= - //derivative value from c real q; Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. public: Elast icMat or ialPropertyO; "ElasticHaterialPropertyQ; // getters //= = = ====r =s====s==nn g------------------------------------------------,= = real getMassDensityO const - C return density;} real getDampingO const { return damping; } real getQvalue const { return q; } / / ------^ - ^ ========^======================^ ====== // virtual function: allow to overwrite on sub class / / ...... =3===^ =.^================= : =^==:======= virtual Matt getKaterialNatrizO const; virtual void getMaterialMatrixCMat teat); / / ====SJ= S= = ======== =S===SS=======S==== = ===S===== == //operation + and ♦ , //=== = := = = = = = = -: ■ ■ : J= r . ,= = = = = = = = = = = // result = this+a + adder*b virtual void add(real a=l.0, ElasticMaterialProperty k adder, real b=l.0, ElasticMaterialProperty k result); void times(real a,ElasticMaterialProperty k result) { add(a,♦this,0.0,result);} virtue void scan(Is is) {}; virtual void print (Os os) O; •define ClassType ElasticMaterialProperty •include <Handle.h> •undef ClassType •endif / / - ♦ - C++ - ♦ - /« Jinbo Chen's Elastietfave /♦ /♦ (c) Copyright 1998 / * Department of Earth Sciences /* Univeristy of Southern California Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. /* All rights reserved. /* / * //#include<Copyright.h> •ifndef RegionGeometry_h_IS_INCLUDED •define RegionGeometry_h_IS_INCLUDEO •include <MenuUDC.h> •include <macros.h> •include <gendass.h> class Is; class Os; class RegionGeometry : public Handleld protected: int nsd; int zlndex; //z-index : concept borrowed from Computer Graphics Handle(RegionMaterial) property; RegionGeometry(); public: "RegionGeometry(); virtual BooLean contains(Ptv t p) const; int getZindexC) const { return zlndex;} //= ====== ==== === ==== = = ====^ = = ===^ ^ ^ ==== = ===== // MenuSystem / / ==============================================r ====^ virtual void scanCIs is); virtual void printCQs os); > •define ClassType RegionGeometry •include <Handle.h> •undef ClassType •endif / / C++ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. / /* Jinbo Chen’s ElasticWave /* /* (c) Copyright 1998 / * Department of Earth Sciences /* Univeristy of Southern California /* All rights reserved. /* /* //#include<Copyright.h> #ifndef RegionMaterial_h_IS_INCLUDED #def ine RegionMaterial_h_IS_INCLUDED •include <NenuUDC.h> •include <FEM.h> •include <LinEqAdm.h> •include <TimePrm.h> •include <GridFE.h> •include <DegFreeFE.h> •include <Arrays_real. h> •include <Store4Plott ing.h> •include <FieldFormat.h> // classes to get Material property / / [ l . l ] uniform is straigt-forw ard / / Cl-2] lin e a r in terp olation / / --------------------------------------------------------------------------------- // [2] circular linear interpolation : this feature is very useful // basin simulation class RegionMaterial : public Handleld i protected: int nsd; //Dimension BooLean polar; //polar distribution Ptv(real) sPoint; //starting point Ptv(real) ePoint; //ending point Handle(ElasticMaterialProperty) sProperty; Handle (Elast icMaterialProperty) eProperty; //derivative attribute based on ElasticMaterialProperty; BooLean uniform; Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. / / — =s = = = „ = = g a = _ //sort of hack to improve efficiency a little, //its dimension is defined by GridFE. / / =====, ------------ Ptv(real) p; virtual redim(int _nsd ); public: RegionMateria(int _nsd) { redim(.nsd) >; ~RegionMaterial(); //======= : =.^ ==s ==a . J= rJ , ===s===- ===s===g===== // virtual function: allow to overwrite on sub class / / r==S ====r-T r—= i . - » . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . = = = = = = = = = virtual ElasticMaterialProperty fcgetElasticMaterialProperty(Ptvt pt); virtual void getgetElasticMaterialPropertyCPtvt pt, Mat k mat); virtual real getMassDensity(Ptvft pt) const; virtual real getDamping(Ptv k pt) const; virtual void getMaterialMatrizCPtv k pt,BooLean isotropic=DpTRUE, Mat k mat); / / = = =3===s===s===s== =s===s = = =============== // MenuSystem --------------------- virtual void scan(Is is); virtual void print(Os os); / / »===»==»==, =,-. r== = = = = = = = == s= 3= ^ = =s==3a= = „ = = - a= a //only function: set elements within this domain of grid //to have the sub domain be the region.no / / r = = = = == r ^= = S = = =s===s= = :=s= = = S S S =s= =rS = = = B = = #define ClassType MaterialRegion ♦include <Handle.h> ♦undef ClassType tendif // C++ /* Jinbo Chen’s ElasticWave */ / * * / /* (c) Copyright 1998 */ 132 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. /* Department of Earth Sciences */ /* Univeristy of Southern California */ /* All rights reserved. */ / * * / / * * / //#include<Copyright.h> •ifndef RegionPolygon.h.IS.INCLQDED •define RegionPolygon_h_IS_INCLUDED •include <NenuUDC.h> •include <Arrays_real.h> •include <ArrayGenSimplest.h> •include <Store4Plotting.b> // 2D or 3D Box region class RegionPolygon : public RegionGeometry { protected: int nsd; //Dimension int np; //number of verties ArrayGenSimplest(Ptv) pts; //Polygon verties public: RegionPolygon(); 'RegionPolygon(); virtual BooLean contains(Ptv t p) const; // MenuSystem / / = =========== ===. === == = == = == = = ==== = ==. = ========= = = virtual static void defineStatic(MenuSystem ftmenu, int level = MAIN); virtual void define(MenuSystem tmenu, int level = MAIN) virtual void scan(HenuSystem* menu); virtual void adm(HenuSystem f t menu); / / ==~===========--- ^============= = = -==== = == == = =~= = = //only function: set elements within this domain of grid // to have the sub domain be the region.no / / ==================3== ======------- r----- > •define ClassType RegionPolygon •include <Handle.h> •undef ClassType 133 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ♦endif // C++ -*• Jinbo Chen's ElasticVave (c) Copyright 1998 Department of Earth Sciences Univeristy of Southern California All rights reserved. //#indude<Copyright. h> ♦ifndef RegionBox_h_IS_INCLUDED #define RegionBox_h_IS_INCLUDED ♦include <HenuUDC.h> ♦include <FEN.h> ♦include <LinEqAda.h> ♦include <TioePrm.h> ♦include <GridFE.b> ♦include <DegFreeFE.h> ♦include <Arrays_real.h> ♦include <Store4Plotting.h> ♦include <FieldFormat.h> // 2D or 3D Box region class RegionBox: public RegionGeometry protected: Ptv(real) xMin; Ptv(real) xMax; public: RegionBox(); 'RegionBox(); virtual BooLean contains(Ptv t p) const; / / = = = = = = = = = = = = = T = = ^ - - = = = = S = = = = S = ----------------- // HenuSystea . virtual void scanCls is); virtual void print(0s os); > 134 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. •endif // C++ /* Jinbo Chen’s ElasticVave */ / * ♦/ /* (c) Copyright 1998 */ /* Department of Earth Sciences */ /* Univeristy of Southern California */ /* All rights reserved. •/ / * * / / * * / //#include<Copyright.h> •ifndef RegionElliptic.h.IS.INCLUDED #def ine RegionElliptic_h_IS_INCLUDED •include <MenuUDC.h> #include <FEN.h> #include <LinEqAdm.h> •include <TimePrm.h> •include <GridFE.h> •include <DegFreeFE.h> •include <Arrays_real.h> •include <Store4Plotting.h> •include <FieldFormat.h> // 20 or 3D Box region class RegionElliptic : public RegionGeometry • C protected: int nsd; Ptv(real) xCenter; Ptv(real) xAxis; public: RegionElliptic(); 'RegionElliptic(); virtual BooLean contains(Ptv k p) const; / / = = = = = = = = = = = = = = « r ..................._ = - - = = - - = _ - - -------------------------- // MenuSystem //====== = = « ======.- --------------- virtual static void def ineStatic(MenuSystem kmenu, int level = MAIN); virtual void def ine (NenuSy stem taenu, int level = MAIN) t def ineS tat ic (menu, level) ; } ■ virtual void scan(HenuSysteafc menu); 135 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. virtual void adm(MenuSystem t menu); //== = - - - - - - - - -,- - - - - - -~ r- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - = ■ - . - r r ■ _ //only function: set elements within this domain of grid // to have the sub domain be the region_no =S S = = S = = = , . = = = = B = = = = S=S= > ♦define ClassType RegionGeometry ♦include <Handle.h> ♦undef ClassType ♦endif // C++ -*- /* Jinbo Chen's ElasticVave * / / * * / /* (c) Copyright 1998 * / /* Department of Earth Sciences «/ /* Univeristy of Southern California */ /* All rights reserved. * / /* * / / • * / ............................ //#include<Copyright.h> ♦ifndef HaterialDomain_h_IS_INCLUDED ♦define NaterialDomain.h.IS.INCLUDED ♦include <macros.h> ♦include <genclass.h> class Is; class Os; ♦include <MenuUDC.h> ♦include <GridFE.h> ♦include <DegFreeFE.h> ♦include <Arrays.real.h> ♦include <Store4Plotting.h> class MaterialDomain : public Handleld protected: / / ======= ^============s= ==ss= ===s===ssa=ss= //the separation of property and region is to allow //different region geometry and material interoplation //says, circle region and uniform property 136 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. //or circle region but linearly increased sith depth / / ==^-=== Handle(RegionGeometry) geometry; / / = S=S= S= SS===-.=T-------------- ------------------------- //sort of hack to improve efficiency a little, //its dimension is defined by GridFE. / / ====« = =— ......................... Ptv(real) p; public: MaterialDomainO; 'MaterialDomainO; //= = = = = // setter / / = = = = virtual const RegionHaterial * regionHaterial () const { return propertyO;} virtual RegionHaterial * regionHaterial () { return propertyO ;> virtual const RegionGeometry * regionGeometry () const { return geometry() > virtual RegionGeometry k regionGeometry () - C return geometry() > // virtual function: allow to overwrite on sub class // e.g., CircleMaterialAegion / / S ==Tr==^n== S==== = ==S==s==s=====B = = = = ==r = ==== == = virtual real getHassDensity(int elm_no) const { return property.getMassDensity(p); > virtual real getDamping(int elm.no) const { return property.getDamping(p); > virtual Mat* getMaterialMatriz(int elm.no) { property.getMaterialMatriz(p); > virtual BooLean contains(int elm.no) { p = grid.getGlobalCentroid(elm.no); return geometry.contains(p); > / / = = = = = // HenuSystem / / = = = = virtual void scan(Is is); virtual void print(Os os); 137 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ♦define ClassType MaterialDomain ♦include <Handle.h> ♦undef ClassType ♦endif Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R eference List [1] J. P. Ake and A. R. Sanford. New evidence for existence and internal structure of a thin layer of magma at mid-crustal depths, New Mexico. Bulletin of the Seismological Society of A merica, 78:1335-1359, 1988. [2] Keiiti Aki. Scattering conversions P to S versus S to P . Bulletin of the Seis mological Society of America, 82:1969-1972, 1992. [3] Keiiti Aki, M. Fehler, and S. Das. Source mechanism of volcanic tremor: Fluid- driven crack model and their application to 1963 Kilauea eruption. J. Volcanol. Geotherm. Res., 2:259-287, 1977. [4] Keiiti Aki and V. Ferrazzini. A snap-shot mapping of localized heterogeneities in the lithosphere using the coda wave data from a not-so-dense seismic array, submitted to geophts. res. lett., 1997. Journal of Geophysical Research, 1998. [5] Keiiti Aki and Valerie Ferrazzini. Seismic signals from La Reunion. Technical report, 1996. [6] P. H. Bachelery, P.A. Blum, L. Chevallier J.L. Cheminee, R. Gaulon, N. Gi- rardin, C. Jaupart, X. Lalanne, J.L. LeMouel, J.C. Ruegg, and P.M. Vincent. 139 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Eruption at Piton de la Foumaise volcano on 3 February, 1981. Nature, 297:395-397, 1982. [7] S. T. Balesta, A. A. Kargopoltsev, and G. B. Grigoryan. Seismic data on magma chambers of the large Tolbachik eruption. Bull. Volcanol, 41:473-479, 1978. [8] J. J. Barton and L. R. Nackman. Scientific and Engineering C++ - An Intro duction with Advanced Techniques and Examples. Addison-Wesley, 1994. [9] A. Bonneville, J. P. Barriot, and R. Bayer. Evidence from geoid data of a hotspot origin for the southern Mascarene plateau and Mascarene islands (Indian Ocean). Journal of Geophysical Research, 93:4199-4212, 1988. [10] G. M. Brown. The layered ultrabasic rocks of Rhum, Inner Hebrides. Philos Trans. R. Soc. London, Ser. B, 240:1-53. 1956. [11] A. M. Bruaset. Tools for automatic generation of result reports. The Numerical Objects Report Series #1997:11, Numerical Objects AS, Oslo, Norway, October 6, 1997. See ftp://ftp.nobjects.com/pub/doc/N097-ll.ps.gz. [12] T. Cargill. C++ Programming Standard. Addison-Wesley, 1992. [13] F. W. Cater. Chilled contacts and volcanic phenomena associated with the Cloudy Pass batholith, Washington. U.S. Geol. Surv. Prof. Pap.. 400-B:B471- B473, 1960. 140 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [14] B. Chouet. Ground motion in the neax field of a fluid driven crack and its interpretation in the study of shallow volcanic tremor. Journal of Geophysical Research, 86:5985-6016, 1981. [15] B. Chouet. Excitation of a buried magmatic pipe: A seismic source model for volcanic tremor. Journal of Geophysical Research, 90:1881-1893, 1985. [16] J. S. Collier and M. C. Sinha. Seismic mapping of a magma chamber beneath the Valu Fa Ridge, Lau Basin. Journal of Geophysical Research, 97:14031- 14053, 1992. [17] C. E. Corry. Laccoliths; mechanics of emplacement and growth. Spec. Pap. Geol. Soc. Am, 220, 1988. [18] R. S. Crosson and D. A. Bame. A spherical source model for low-frequency volcanic earthquakes. Journal o f Geophysical Research. 90:11184-11198, 1985. [19] M. A. Dalain. The application of high-order differencing to scalar wave equation. Geophysics, 58:54-66, 1986. [20] P. T. Delaney, R.S. Fiske, A. Milklius, A. T. Okamura, and M. K. Sato. Deep magma body beneath the summit and rift zones of Kilauea volcano, Hawaii. Sciences, 247:1311-1316, 1990. [21] H. Delorme, P. Bachelery, P. A. Blum, J. L. Chemince, J. F. Delarue, J. C. Delmond, A. Him, J. C. Lepine, P. M. Vincent, and J. Zlotnicki. March 1986 eruptive episodes at Piton de la Fournaise volcano (Reunion Island). J. Vol- canol. Geotherm. Res., 36:199-208, 1989. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [22] R. S. Detrick, P. Buhl, E. Vera, J. Mutter, J. Orcutt, J. Madsen, and T. Brocher. Multichannel seismic imaging of a crustal magma chamber along the East Pacific Rise. Nature, 326:35-41, 1987. [23] J. H. Dieterich. Growth and persistence of Hawaiian rift zones. Journal of Geophysical Research, 93:4258-4270, 1988. [24] D. Dzurisin, R. Y. Koyanagi, and T. T. English. Magma supply and storage at Kilauea volcano, Hawaii, 1956-1983. J. Volcanol. Geotherm. Res., 21:177-206, 1984. [25] P. Einarsson. S-wave shadows in the Krafla caldera in NE Iceland, Evidence for a magma chamber in the crust. Bull. Volcanol, 41:187-195, 1978. [26] J. R. Evans and J. J. Zucca. Active high-resolution seismic tomography of compressional wave velocity and attenuation structure at Medicine lake vol cano, northern California Cascade Range. Journal of Geophysical Research, 93:15016-15036, 1988. [27] J. W. Ewert and D. A. Swanson (Eds). Monitoring volcanoes: Techniques and strategies used by the staff of the Cascades Volcano Observatory, 1980-90. if. S. Geol. Suru. Bull, 1966:223pp, 1992. [28] M. C. Fehler. Observation of volcanic tremors at Mount. St.Helens volcano. Journal of Geophysical Research, 88:3476-3484, 1983. [29] M. C. Fehler and B. Chouet. Operation of a digital seismic network at Mount St. Helens volcano and observation of long-period events originated under the volcano. Geophy. Res. Lett, 9:1017-1020, 1982. 142 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [30] M. G. Ferrick and W. F. St. Lawrence. Comment on “Observations of vol canic tremor at Mount St. Helens volcano” by Michael Fehler. Journal of Geophysical Research, 89:6349-6350, L984. [31] M. G. Ferrick, A. Qamar, and W. F. St. Lawrence. Source mechanism of vol canic tremor. Journal of Geophysical Research, 87:8675-8683, 1982. [32] F. Ferrucci, A. Hirn, G. De Natale, J. Virieux, and L. Mirabile. P-SV conver sions at a shallow boundary beneath Campi Flegrei caldera (Italy): Evidence for the magma chamber. Journal of Geophysical Research, 97:15351-15359, 1992. [33] R. S. Fiske, C. A. Hopson, and A. C. Waters. Geology of Mount Rainier National Park, Washington. U. S. Geol. Surv. Prof. Pap, 444:93pp, 1963. [34] J. J. Frankel. Forms and structures of intrusive basaltic rocks, in Basalts: The Polderuaart Treaties on Rocks of Basaltic Composition. Wiley-interscience. New York, 1967. [35] J. R. Grasso and P. Bachelery. Hierarchical organization as a diagnostic ap proach to volcano mechanics: Validation on Piton de la Fournaise. Geophys. Res. Letters, 22:2897-2900, 1995. [36] A. Gudmundsson. Form and dimensions of dykes in eastern Iceland. Tectono- physics, 95:295-307, 1983. [37] D. P. Hill. Earthquakes and carbon dioxide beneath Mammoth Mountain, California. Seismol. Res. Letters, 67:8-15, 1996. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [38] A. Him, J.C. Lepine, M. Sapin, and H. Delorme. Episodes of pit-crater collapse documented by seismology at Piton de la Fournaise. J . Volcanol. Geotherm. Res., 47:89-104, 1991. [39] A. Him, A. Nercessian, J.C. Lepine, and M. Sapin. Internal structure of Piton de la Fournaise volcano from seismic wave propagation and earthquake distri bution. J. Volcanol. Geotherm. Res., in press, 1995a. [40] A. Hirn, M. Sapin, J.C. Lepine, and A. Nercessian. Failure and fluid flow from earthquakes of eruptions at Piton de la Fournaise volcano. J. Volcanol. Geotherm. Res., in press, 1995b. [41] J. P. Iddings. Problems of volcanism. Yale University Press, New Haven, Conn., 1914. [42] K. Ishihara. Dynamical analysis of volcanic explosion. J. Geodynamics, 3:327- 349, 1985. [43] H. M. Iyer. Seismological detection and delineation of magma chambers mea surements beneath intraplate volcanic centers in western U.S.A., in Modeling of Volcanic Processes. Vieweg, Wiesbaden, Germany, 1988. [44] B. Joe. Geompack - a software package for the generation of meshes using geometric algorithms. Adv. Eng. Software, 13:325-331, 1991. [45] B. Joe. Three-dimensional boundary-constrained triangulations. In E. N. Houstis and J. R. Rice, editors, Artificial Intelligence, Expert System, and Sym bolic Computing, pages 215-222. Elsevier Science Publishers, 1992. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [46] B. Joe. Construction of A:-dimensional delaunay triangulation using local trans formations. SIAM J. Sci. Comput. 14, November 1993. [47] B. Joe. Tetrahedral mesh generation in polyhedral regions based on convex polyhedron decompositions. Int. J. Meth. Engng., 1993. [48] A. M. Johnson and D. D. Pollard. Mechanics of growth of some laccolithic intrusions in Henry Mountains, Utah. Tectonophysics, 18:261-309, 1973. [49] K. Kamo. Some phenomana before summit eruption at Sakurajima volcano. Bull. Volcanol. Soc. Japan, 23:53-64, 1978. [50] S. Koyanagi, K. Aki, N. Biswas, and K. Maeda. Infered attenuation from site effect-corrected T-phase recorded on the island of Hawaii. PEGEOPH, 144:1- 18, 1995. [51] A. Kubotera. Volcanic tremors at A so volcano, Physical volcanology. Elsvier Amsterdam, first edition, 1974. [52] H. P. Langtangen. Efficient finite element solution of the linear wave equation in Diffpack. World Wide Web document: Diffpack vl.4 Report Series, SINTEF & University of Oslo, 1996. U RL: http://www.nobjects.com/diffpack/reports. [53] H. P. Langtangen. An example on visualizing scalar and vector fields in Diff- pack. The Numerical Objects Report Series #1997:6, Numerical Objects AS, Oslo, Norway, October 6, 1997. See ftp://ftp.nobjects.com/pub/doc/N097~ 06.ps.gz. 145 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [54] H. P. Langtangen. Improving the efficiency of Diffpack simulators for PDEs with time independent coefficients. The Numerical Objects Report Se ries #1997:8, Numerical Objects AS, Oslo, Norway, October 6, 1997. See ftp://ftp-nobjects.com/pub/doc/NO97-08.ps.gz. [55] H. P. Langtangen. Plotting of curves in Diffpack. The Numerical Objects Report Series #1997:4, Numerical Objects AS, Oslo, Norway, October 6, 1997. See ftp://ftp.nobjects.com/pub/doc/NO97-04-ps.gz. [56] H. P. Langtangen, G. Pedersen, and V V . Shen. Finite element preprocessors in Diffpack. The Numerical Objects Report Series #1997:10, Numerical Objects AS, Oslo, Norway, October 6, 1997. See ftp://ftp.nobjects.com/pub/doc/N097- I0.ps.gz. [57] J. H. Latter. Volcanological observations at Tongariro National Park, 2, Types and classification of volcanic earthquakes, 1976-1978. N .Z Dept, of Sci. and Ind. Res. Geophysics Div., 1979. [58] W. St. Lawrence and A. Qamar. Hydraulic transients: A seismic source in volcanoes and glaciers. Sciences, 203:654-656, 1979. [59] J. A. Lehman, R. B. Smith, M. M. Schilly, and L. W. Braile. Upper crustal structure of the Yellowstone caldera from seismic delay time analyses and grav ity correlations. Journal of Geophysical Research, 87:2713-2730, 1982. [60] J. F. Lenat, P. Bachelery, A. Bonneville, and A. Hirn. The beginning of the 1985-1987 eruptive cycle at Piton de la Fournaise (La Reunion); New insights 146 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. in the magmatic and volcano-tectonic systems. J. Volcanol. Geotherm. Res., 36:209-232, 1989c. [61] J. F. Lenat, P. Bachelery, A. Bonneville, P. Tarits, J.L. Cheminee, and H. De lorme. The December 4, 1983 to February 18, 1984 eruption of Piton de la Fournaise (La Reunion, Indian Ocean): Description and interpretation. J. Volcanol. Geotherm. Res., 36:87-112, 1989b. [62] J. F. Lenat, P. Vincent, and P. Bachelery. The off-shore continuation of an active basaltic volcano: Piton de la Fournaise (Reunion Island, Indian Ocean); structural and geomorphological interpretation from sea beam mapping. J. Volcanol. Geotherm. Res., 36:1-36, 1989a. [63] J. B. Lowenstern, P. C. Wallmann, and D. D. Pollard. The W est Mageik Lake sill complex as analogue for magma transport during 1912 eruption at Valley of Ten Thousands Smokes, Alaska. Geophiscal Res. Lett., 18:1569-1572, 1991. [64] T. Minakami. Fundamental research for predicting volcanic eruption (i). Bull. Earth Res. Insti. Univ. Tokyo, 38:497-544, 1960. [65] K. Mogi. Relations between the eruptions of various volcanoes and the defor mation of the ground surface around them. Bull. Earth Res. Insti. Univ. Tokyo, 36:99-134, 1958. [66] J. Mori, H. Patia, C. Mckee, I. Itikarai, P. Lowenstein, P. de Saint Ours, and B. Talai. Seismicity associated with eruptive activity at Langila Volcano, Papua New Guinea. J. Volcanol. Geotherm. Res., 38:243-255, 1989. 147 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [67] J. B. Murray, A. D. Pullen, and S. Saunders. Ground deformation surveying of active volcanoes, in monitoring active volcanoes: Strategies, procedures, and techniques. UCL Press, pages 113-150, 1994. [68] J. B. Murray, A. D. Pullen, and S. Saunders. Real-time ground deformation monitoring, in monitoring active volcanoes: Strategies, procedures, and tech niques. UCL Press, pages 92-112, 1994. [69] K. Nakamura. Why do long rift zones develop in Hawaiian volcanoes- A possi ble role of thick oceanic sediments. Bull. Volcanol. Soc. Japan, Ser. 2:255-270, 1980. [70] R. J. O’Connell and B. Budiansky. Viscoelastic properties of fluid-saturated cracked solids. Journal of Geophysical Research, 82:5719-5735, 1977. [71] Y. Okada. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America. 75:1135-1154, 1985. [72] F. Omori. The Usu-san eruption and the earthquake and elavation phenom ena,II. Bull. Impl. Earthquake Invest. Comm., 5:101-137, 1913. [73] F. Omori. The Sakurajima eruption and earthquakes. Bull. Impl. Earthquake Invest. Comm., 8:1-525, 1918. [74] V V . H. Parsons. Volcanic centers of the Sunlight area Park County, Wyoming. Journal of Geology, 47:1-20, 1939. 148 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [75] J. P. Rancon, P. Lerebour, and T. Auge. The Grand Brule exploration drilling: new data on the deep framework of the Piton de la Fournaise vocano, P art 1. J. Volcanol. Geotherm. Res., pages 113-128, 1989. [76] P. M. Roberts, Keiiti Aki, and M. C. Fehler. A low velocity zone in the basement beneath the Valles caldera, New Mexico. Journal o f Geophysical Research, 96:21583-21596, 1991. [77] D. Rousset, A. Lesquer, A. Bonneville, and J. F. Lenat. Complete gravity study of Piton de la Furnaise volcano,Reunion Island. J. Volcanol. Geotherm. Res., 36:37-52, 1989. [78] L. R. Rowan and R. W. Clayton. The three dimensional structure of Kilauea volcano, Hawaii, from travel time tomography. Journal of Geophysical Re search, 98:4355-4375, 1993. [79] M. Ruiscetti, R. Schick, and D. Seidl. Spectral parameters of volcanic tremors at Mt. Etna. J. Volcanol. Gerotherm. Res., 2:289-298, 1977. [80] C. O. Sanders. Location and configuration of magma bodies beneath Long Valley, California, determined from anomalous earthquake signals. Journal of Geophysical Research. 89:8287-8302, 1984. [81] K. Sassa. Volcanic micro-tremors and eruption earthquakes (I). Mem. Coll. Sci. Univ. Kyoto, 18:255-293, 1934. [82] K. Sassa. Micro-seismometric study on eruption of the volcano Aso (II). Mem. Coll. Sci. Univ. Kyoto, 19:11-56, 1936. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [83] A. D. L. Sharp, P. M. Davis, and F. Gray. A low velocity zone beneath E tna and magma storage. Nature, 287:587-591, 1980. [84] M. Shima. On the second volcanic tremor at the Aso volcano. Bull. Disaster Prevent. Inst. (Jniv. Kyoto, 22:1-5, 1958. [85] M. Shima. Volcanic micro-tremor at the Aso volcano. Bull. Disaster Prevent. Inst. Univ. Kyoto, 34:2-17, 1960. [86] D. W. Steeples and H. M. Iyer. Low velocity under the Long Valley as deter mined from teleseismic events. Journal of Geophysical Research, 81:849-860, 1976. [87] C. H. Thurber. Seismic detection of the summit magma complex of Kilauea, Hawaii. Sciences. 223:165-167, 1993. [88] C. H. Thurber and A. E. Gripp. Flexure and seismicity beneath the south flank of Kiauea volcano and tectonic implications. Journal of Geophysical Research. 93:4271-4278, 1988. [89] R. I. Walcott. Flexure of the lithsphere at Hawaii. Tectonophysics, 9:435-446, 1970. [90] G. P. L. Walker. The dike complex of Koolau volcano, Oahu: Internal structure of a Hawaiian rift zone, in Volcanism in Hawaii. USGS Prof. Pap., 1350:961- 993, 1987. 150 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [91] Zeng Yuhua. Theory of scattered P- and S-wave energy in a random isotropic scattering medium. Bulletin of the Seismological Society of America, 83:1264- 1276, 1993. [92] J. J. Zucca and J. R. Evans. Active high-resolution compressional wave atten uation tomography at Newberry volcano, central Casacade Range. Journal of Geophysical Research, 97:11047-11055, 1992. 151 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Using small events to measure the evolution of the strain and stress fields before the 1992 Landers, California earthquake

PDF

A unified methodology for seismic waveform analysis and inversion

PDF

Modeling and imaging asperities on a fault plane and characterizing spatial and temporal patterns of precursory seismicity

PDF

Modeling of continuous tiltmeter data from the 1984 rifting event at Krafla Volcano, Iceland

PDF

Radium isotopes in San Pedro Bay, California: Constraint on inputs and use of nearshore distribution to compute horizontal eddy diffusion rates

PDF

Management of large earthquake data using the Antelope Relational Database and seismicity analysis of the 1999 Turkey earthquake sequences

PDF

Systematic high -resolution imaging of fault zone structures

PDF

Paleoecology and depositional paleoenvironments of Pleistocene nearshore deposits, Las Animas, Baja California Sur, Mexico

PDF

The unusual sedimentary rock record of the Early Triassic: Anachronistic facies in the western United States and southern Turkey

PDF

Marine paleoecology during the aftermath of the end-Permian mass extinction

PDF

The structure and development of Middle and Late Triassic benthic assemblages

PDF

Deterministic and stochastic modelling of the high frequency seismic wave generation and propagation in the lithosphere

PDF

Mesozoic intraplate deformation in the east Asian tectonic collage: The enigmatic northwest Ordos region, China

PDF

Simulating nonplanar fault surfaces using stochastic branching

PDF

Numerical analysis of slope stabilization concepts using piles

PDF

A tectonic model for the formation of the gridded plains on Guinevere Planitia, Venus: Implications for the thickness of the elastic lithosphere

PDF

Reduced-parameter modeling for cost estimation models

PDF

A gravity and magnetic study of the Tehachapi Mountains, California

PDF

Numerical simulations of swirling flows

PDF

Multi-view image -based rendering and modeling

Asset Metadata

Creator Chen, Jinbo (author)

Core Title Numerical modeling of long period events at the Piton de la Fournaise by the use of object-oriented scientific computing

Degree Doctor of Philosophy

Degree Program Geological Sciences

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag geophysics,OAI-PMH Harvest

Language English

Contributor Digitized by ProQuest (provenance)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-172219

Unique identifier UC11329241

Identifier 3054858.pdf (filename),usctheses-c16-172219 (legacy record id)

Legacy Identifier 3054858.pdf

Dmrecord 172219

Document Type Dissertation

Rights Chen, Jinbo

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA