Source mechanism of volcanic tremor: fluid-driven crack models and their application to the 1963 kilauea eruption

Journal of Volcanology and Geothermal Research, 2 ( 1 9 7 7 ) 2 5 9 - - 2 8 7

259

© Elsevier Scientific P u b l i s h i n g C o m p a n y , A m s t e r d a m - - P r i n t e d in T h e N e t h e r l a n d s

SOURCE MECHANISM OF VOLCANIC TREMOR: FLUID-DRIVEN CRACK MODELS AND THEIR APPLICATION TO THE 1963 KILAUEA ERUPTION

KEIITI AKI, M I K E F E H L E R

and S H A M I T A

DAS

Department of Earth and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, Mass. 02139 (U.S.A.) (Received November 1, 1976; revised version accepted March 17, 1977)

ABSTRACT

Aki, K., Fehler, M. and Das, S., 1977. Source mechanism of volcanic tremor: fluid-driven crack models and their application to the 1963 Kilauea eruption. J. Volcanol. Geotherm. Res., 2: 259--287. W e propose a model for the mechanism of m a g m a transport based on a fluid-filled tensile crack driven by the excess pressure of fluid. Such a transport mechanism can generate seismic waves by a succession of jerky crack extensions, if the fracture strength of rock varies in space, or if there is a difference between the dynamic and static values of the criticalstress intensity factor. W e also find that the opening and closing of a narrow channel connecting two fluid-filledcracks m a y be a source of seismic waves. Using the finite-difference method, we calculated the vibration of dry and fluid-filledcracks generated by: (1) a jerky extension at one end or at both ends and (2) a jerky opening of a narrow channel connecting two cracks. W e then calculated the far-fieldand near-field radiation from these vibrating cracks. The spectra show peaked structures, but interestingly, most high-frequency peaks are only present in the near-field and cannot be transmitted to the far-field.The spectral features described above are often observed for volcanic tremors and in some cases for seismic signals associated with hydraulic fracturing experiments. W e firstconsider as a model of volcanic tremor randomly occurring jerky crack extensions, and derive a formula relating the tremor amplitude to the excess pressure in the magma, the incremental area in each extension, and the frequency of extensions. These parameters are also constrained by other observations, such as the rate of m a g m a flow. Our model was tested quantitatively against observations m a d e in one of the bestdescribed case histories of volcanic tremor: the October 5--6, 1963 Kilauea flank eruption. W e found that a single, long crack extending from the summit to the eruptive site cannot explain the observations. The model of a steadily expanding crack ran into difficulties w h e n quantitative comparisons were made with observations. The extension of crack area needed to explain the amplitude of volcanic tremor should accompany a large increase in tremor period which was not observed. Our second model is a chain of cracks connected by narrow channels which open and close. The length of each crack is around 1 kin, the channel area connecting neighboring cracks is about 103 m 2, and the channel opens jerkily with the magmatic excess pressure of about 20 bars. The frequency of jerky opening of each channel is about once in 15 seconds. The channel is closed after each jerky opening, as soon as m a g m a is m o v e d through the channel.

260

INTRODUCTION Seismic signals generated by volcanic activities have been investigated for m a n y years since the pioneering work of Omori (1912). These signals are quite variable in character, ranging from those indistinguishable from tectonic earthquakes to continuous vibrations called volcanic tremor, pulsations or harmonic tremor (Minakami, 1974). The earthquake-like signals can be interpreted by means of a well-developed m e t h o d for finding the location and source mechanism of tectonic earthquakes. The m e t h o d for interpreting volcanic tremor, however, is still at a rudimentary stage because of the lack of an adequate physical model of its source mechanism. Our purpose is to propose a model of volcanic tremor which may be used in determining physical parameters for the process of magma transport when a volcano is in action. We shall consider here the mechanism of magma transport to be a crack filled with magma that is driven by the excess pressure in the magma. The idea of a fluid-filled crack for a model of magma transport is not new. For example, Weertman (1971a, b) proposed a mechanism in which a crack filled with magma nucleates at the b o t t o m of the lithosphere, gradually increases in length and volume, pinches its lower end shut, and finally rises towards the top surface by the force of buoyancy. A fluid-filled crack plays an important part in a theory proposed to explain the pattern of radial dikes and distribution of flank eruptions. Nakamura (1977) likens a volcano to an experimental site of repeated large-scale hydraulic fracturing by pressurized magma, and demonstrates its usefulness in determining the orientation of insitu tectonic stress. In the following, we shall first review seismic observations of volcanic tremor. Their similarity to seismic signals associated with hydraulic fracturing will also be mentioned. We shall then formulate a seismic source problem in terms of the parameters of a fluid-filled crack driven by internal fluid pressure We consider two possible source mechanisms: (1) a fluid-filled crack growing by jerky extension at the crack tip, and (2) two fluid-filled cracks connected by a narrow channel which opens and closes due to excess pressure in the magma. We shall also discuss application of our models to volcanic tremor observed during an eruption of Kilauea volcano. The observed amplitude of the tremor and the peak frequency of tremor spectrum, together with the magma flow estimated from subsidence of the summit, enable a rough quantitative estimate of the model parameters. OBSERVED CHARACTERISTICS OF VOLCANIC TREMOR Let us first review some important features of volcanic tremor observed at various volcanoes to see if the idea of a fluid-driven crack is adequate as its source mechanism. The most conspicuous feature of volcanic tremor is that the spectrum has

261

a dominant peak. The peak frequency varies from one volcano to another, ranging from 0.1 to 10 Hz (Kubotera, 1974). Several peaks are observed for a single volcano. For example, Sassa (1935) finds four distinct tremors at Mt. Aso, with predominant periods at 0.2, 0.4--0.6, 1, and 3.5--7 s. The Fourier spectra calculated for a few cases (e.g., Kubotera, 1974; Shimozuru et al., 1966) indicate that the spectral peaks have a width corresponding to an oscillator with Q of 5--10. The cause of the peaked structure in the tremor spectra has been a subject of speculation for many years. There are t w o main lines of thought: one attributing the spectral peaks to a source effect, such as free oscillations of a magma chamber; and the other to a path effect, such as reverberations in the stratified structure of a volcano. Evidence in favor of a path effect is that the main phase of a volcanic earthquake of shallow focus, such as Minakami's B-type earthquake (occurring within a kilometer of the summit in the case of Mt. Asama), and volcanic tremor share the same spectrum, propagation velocity and mode of particle motion (Watanabe, 1961). The most striking observation in favor of the path effect is the similarity between the seismogram of a 10-kg yield explosion in a lake and that of volcanic tremor observed by Dibble (1974) at Mt. Ruapehu. In these seismograms, the maximum amplitude probably corresponds to surface waves, which show a long oscillation with a certain predominant period that has nothing to do with the seismic source. Wada and Ono (1963) and Kikuchi (1964) showed, using independent techniques, that the volcanic tremor in the frequency range 2--5 Hz observed at Mt. Aso are primarily surface waves. On the other hand, there are two important pieces of evidence supporting the assertion that the spectral peaks of volcanic tremor are due to some source effects. One is the temporal change of peak frequency during the course of volcanic activity. For example, Kubotera (1974) shows a significant change of peak frequency of long-period tremors at Mt. Aso (type 2 of Sassa, 1935) with time. As discussed later, Shimozuru et al. (1966) observed a change in peak frequency during the October 5--6, 1963 eruption of Kilauea volcano. Another piece of evidence in support of the source effect is the presence of oscillatory motion with nearly constant frequency from the beginning of P-waves and throughout the record of an isolated event. When those events occur frequently, they merge into a continuous volcanic tremor. An example of such a case is the type-C tremor of Mt. Sakurajima, with periods in the range 0.5--1.3 s, described by Kakuta and Idegami (1970). Another clear example is the tremor classified as LPD (standing for long-period deep, with period range of 0.2--1.0 s) at the Hawaiian Volcano Observatory (R.Y. Koyanagi, personal communication), which can only be e~;plained as due to some oscillatory sources. Schick and Riuscetti (1973), who studied the volcanic tremor of Italian volcanoes, also favor the idea that the source effect is responsible for spectral peaks. They point out that a substantially similar shape of spectrum was ob-

262 served at stations in different geological settings, and the spectra due to other seismic sources did not show the peaks observed in volcanic tremor. An outstanding feature of volcanic tremor reported by several observers is the strong high-frequency oscillations in the immediate vicinity of the eruption site. Shimozuru et al. (1958) operated a seismograph at the b o t t o m of the crater of Naka-dake, one of the post-caldera cones of Mt. Aso, and observed continuous oscillations with frequencies of 5--10 Hz. According to R.Y. Koyanagi (personal communication), during the initial phase of a flank eruption as well as the fountaining phase in prolonged eruptive activity in Kilauea, a strong high-frequency {% 5 Hz) tremor predominates at the eruption site. We may summarize the outstanding features of volcanic tremor as follows. (1) The spectrum of volcanic tremor has a peaked structure, which is, at least partially, due to the source effect. (2) The frequency of the spectral peak attributed to the source effect sometimes changes in the course of eruptive activity. (3) High-frequency oscillations are strong in the immediate vicinity of the eruption site. Some of the seismic signals observed during a hydraulic fracturing experiment at the Los Alamos Scientific Laboratory showed a feature similar to volcanic tremor. As reported by Potter and Dennis (1974), the spectral peak of continuous vibration caused by fluid injection changed with time (from 25 Hz to 15 Hz in 5 minutes) in such a way that the frequency decreased when the crack was presumably expanded by an increased fluid injection. Our motivation for the fluid-driven crack model as the source of volcanic tremor came from the recognition of c o m m o n features between volcanic tremor and some of the seismic signals observed during hydraulic fracturing experiments. EARLY MODELS OF VOLCANIC TREMOR Several models have been proposed to explain the peaked spectrum of volcanic tremor. Omer {1950) considered the shallow crust of Kilauea volcano as laminae composed of lava flows, ash beds and intrusives. The magma flow exerts a transverse force on one end of the layers and sets these partially free layers into vibration, in a way similar to bending a deck of playing cards clamped at one end. In Omer's model, the tremor is caused by a quasi-free vibration of solid layers forming the shallow part of the volcano. The free vibration of liquid magma has been proposed to explain the spectral peaks. For example, Shimozuru {1961) and Steinberg (1965) considered the free longitudinal vibration of a magma column in the vent as the cause of volcanic tremor. So far, the most complete analysis of seismic radiation from a liquid vibrator was made by Shima {1958) who considered the spherically symmetric, fundamental mode of vibration in a liquid sphere buried in an u n b o u n d e d elastic medium. Since the wave front is normal to

263

the liquid-solid interface in this case, the seismic radiation is most efficient. For comparison, the longitudinal vibration in a liquid column is a poor seismic radiator if the acoustic velocity in the liquid is less than seismic wave velocities. We shall discuss this point later in more detail. Kubotera (1974) calculated the frequency f0 of a spherical oscillator with the radius a, using Shima's result for various model parameters, and found that f0 is closely approximated by a simple formula: 4.4 al f0 = - - - 2~r a

(1)

where ~1 is the acoustic velocity of the liquid magma. The damping of oscillations of this liquid sphere is determined by the impedance constrast between liquid and solid. Since the wave front is normal to the interface, we expect the result is similar to the plane wave problem for which the amplitude reduction at each reflection is given by a factor (~1 p2 -alPl)/(c~2P2 + 0~lPl), where ~2, P 2 are respectively the acoustic velocity and density of the solid and al, Pl are those of the liquid. Since the reflection occurs f0 times in unit time, the amplitude decay with time can be written as:

e-¢t = ( a2p2-axP, } f°t a2P2 + c~Ipl/ or:

e = foln ( °~202÷°~lpl)_ ~2P2

~IPl

Defining Q of a resonator by: Q = ~fo/e

we obtain:

P2 -ollPl/ This simple f6rmula gives Q-values in good agreement with the numerical results obtained by Kubotera (1974) for a spherical source using Shima's formula. To explain a Q of 5--10 observed for volcanic tremors by this model would require an impedance contrast (that is, a2 P 2/a 1P 1) of 3 : 1 to 6:1. According to the laboratory experiment by Murase and McBimey (1973), the acou,~tic velocity of liquid magma is a b o u t 2 km/s. If the rock outside the magma chamber is c o m p e t e n t and has a P-wave velocity around 6 km/s, then a Q of 5 is possible by this model. If the liquid magma contains vesicles (bubbles), even a small a m o u n t of it will increase the impedance contrast drastically (Aki et al., 1976). There are, however, two unsatisfactory points with Shima's sphere model.

264

GLOSSARY

a b c C d f0 G (cv) G0(co) h k kn l l~ L

Mij n NO N1 Nc p Po

P(co) P(x) q

OF SYMBOLS

radius o f spherical oscillator b u l k m o d u l u s of fluid in c r a c k p h a s e velocity o f surface waves constant crack width f r e q u e n c y of spherical oscillator F o u r i e r t r a n s f o r m o f t o t a l seismic motion F o u r i e r t r a n s f o r m of individual signal d e p t h o f seismic s o u r c e wavenumber wavenumber of nth mode of surface wave crack half-length e i g e n - f u n c t i o n for Love waves crack l e n g t h moment tensor n u m b e r o f jerky e x t e n s i o n s o f c r a c k per u n i t t i m e stress-intensity f a c t o r stress-intensity f a c t o r d u e t o cohesive force d i s t r i b u t i o n P(x) critical s t r e s s - i n t e n s i t y f a c t o r a c o n s t a n t in f o r m u l a ( 3 5 ) f o r Rayleigh wave a m p l i t u d e u n i f o r m h y d r o s t a t i c pressure power spectrum cohesive force d i s t r i b u t i o n at c r a c k tip a c o n s t a n t in f o r m u l a ( 3 5 ) for Rayleigh wave a m p l i t u d e

Q

~fo/e

r

distance between a point on crack a n d receiver

r0

distancefrom crack center to receiver

rl, r 2 S t u v v

e i g e n - f u n c t i o n s for R a y l e i g h wave area o f o n e face o f crack surface time d i s p l a c e m e n t in x - d i r e c t i o n d i s p l a c e m e n t in y - d i r e c t i o n average d i s p l a c e m e n t d i s c o n t i n u i t y over c r a c k s u r f a c e g r o u p velocity ( a s s u m i n g a r e c t a n g u l a r s h a p e for c r a c k s u r f a c e ) l e n g t h o f t h e side of c r a c k surface p e r p e n d i c u l a r t o t h e o n e a l o n g w h i c h L is m e a s u r e d Cartesian co-ordinates

V W

x, y, z

~ ~2 7 e

0 ~,

~ p p P2 aij

c~

compressional wave velocity acoustic velocity in liquid m a g m a acoustic velocity in solid matrix shear wave velocity angle between ~I and direction to receiver characteristic time for decay of amplitude (spherical oscillator) viscosity angle between normal to crack plane and receiver direction wavelength or L a m e constant modulus of rigidity coordinates of point on crack surface coordinate axis along crack length density of m a g m a density of liquid (spherical oscillator) density of solid matrix (spherical oscillator) stress components crack surface angular frequency

265

First, it does n o t naturally offer the driving force of tremors. The model requires a source of pulsation in the liquid. Rapid degassing and phase change may be possible sources, and their seismic effect should be investigated. Secondly, the spherical geometry does not appear to be appropriate for a certain case of volcanic tremor. For example, during the initial phase of flank eruption at Kilauea, a continuous fluid connection seems to exist between the summit magma chamber and the east rift fissures (Swanson et al., 1976a). As we describe later in detail, the evidence for this includes the concurrence of summit deflation and onset of tremor at the eruption sites in the rift zone. Also, the c o n t o u r map of tremor amplitude shows an elongated pattern along the rift zone (Shimozuru et al., 1966). Our fluid-driven crack model will offer both the driving force and the geometry adequate for the magma transport mechanism during a flank eruption at Kilauea volcano. JERKY EXTENSION OF FLUID-DRIVEN TENSILE CRACK

Let us consider a crack filled with fluid in an u n b o u n d e d elastic medium. For simplicity, we shall assume that the crack surface is in the x, z-plane and extends to +~o and -oo in the z-direction. The initial locations of crack tips are at x = + l as shown in Fig. 1. Neglecting gravity and other body forces, we assume that fluid in the crack and solid outside the crack are initially under a uniform hydrostatic pressure P0. With an increase of liquid pressure by Ap, a stress concentration occurs in the elastic body ahead of the crack tip. The tensile stress ayy n e a r the crack tip, say at x = l, may be written (see e.g., Barenblatt, 1962) as:

Oyy

=

AP~/l/[2(x-l)],

fory=Oandx> l

Y

i

l po+AP |%/-' I

I. I

2 ~

"1 I

I

I

I I

I I I

I I

I

I

I

I I i

I I I

I I I

I

I

I

1

Fig. 1. Fluid-filled crack e x t e n d i n g at both ends. Dynamic propagation of tip leaves a cavity behind the tip and the a m b i e n t pressure tends to pinch it shut.

(3)

266 Thus, ayy has a square-root singularity at the crack tip. The coefficient of the singular term is called the stress-intensity factor, which can be written in this case as: No = AP x/1/2

(4)

The crack will extend according to Irwin's criterion if No exceeds a certain critical value Ne. At the initiation of crack, this is equivalent to Griffith's criterion, because Ne is uniquely related to the specific surface energy at zero rupture velocity. Thus, the crack will expand if A P x / I / 2 > Nc

(5)

Suppose that the crack tips moved over a distance Al (Fig. 1). According to Barenblatt (1962), when a crack is maintained open by a flowing viscous fluid injected into it, the fluid does not fill the crack completely. There is always a fluid-free part of the crack on both sides of the wetted area. In our dynamic problem we presume that viscous liquid cannot catch up with the movement of the crack tip, so that there will be an e m p t y space left behind the crack tip over a distance A I. The ambient pressure around the crack tip then acts as the cohesive force which generates a negative stress-intensity factor and resists extension of the crack. The stress-intensity factor due to a cohesive force distribution P(x) is given by:

v~ NI -

l f P(x)

o

dx (6)

x/t2 -x2

where l is the half-crack length after extension by Al, and: P ( x ) = -Po

I-,5l
= 0

1

O
(7)

I

Substituting (7) into (6), we obtain:

- ~ - 0 . 9 Po ~ I / 2

t

(8)

The total stress-intensity factor is the sum of (4) and (8): N = No+N1 ~ APx/I[2-O.9PoX/AI/2

(9)

Thus, the crack which started propagation by the condition No > Nc will be stopped immediately after an infinitesimal increment A l. In other words, a dynamic rupture propagation will occur only quasi-statically without generating seismic waves.

267

The situation may be different, however, if the material strength changes from place to place. The crack tip may be blocked at a point with high N c. When the bond with high Nc is broken, the value of Nc over the next A l may be reduced by ANc. If this reduction in Nc satisfies the following inequality: 0.9 Po x / A l / 2 < A g c

(10)

the crack tip can make a jerky extension over Al. Putting the fractional fluctuation of Nc as e, e = A N c / N c , and Nc ~- ,~P x/ I/2, we rewrite (10) as: e21/Xl > (0.8//) (Pol,SP) 2

(11)

For a random medium, the squared fluctuation e 2 would increase with the sampling distance/x I. The above inequality shows that, for a jerky extension of a crack over Al, the rate of increase of e 2 must be greater than a critical value determined by the crack length l, ambient pressure P0 and liquid overpressure AP. Thus, for a given strength distribution, the jerky extension occurs more easily if I and AP are large and P0 is small. The average distance over which a jerky extension occurs may be determined by the correlation distance of strength fluctuation. For example, if the correlation distance is 1 m, we may put Al = 1 m. Then, for the half-crack length l ~ 1 km, P0 -~ 300 bars (lithostatic pressure at 1 km depth) and Ap _~ 20 bars, the fluctuation e of strength must be greater than 50%. Greater fluctuation is required if the crack is smaller, the depth is greater and fluid excess pressure is lower. Since we expect that material properties become more uniform with increasing depth, a jerky extension of a tensile crack will become less likely. B.V. Kostrov (personal communication, 1976) suggested an alternative mechanism for jerky extension which may operate at great depths. This mechanism requires a difference between the static and dynamic values of the critical stress-intensity factor. Pressure builds up until the critical static stressintensity factor is overcome. Propagation of the crack continues, once begun, because the critical stress-intensity factor is less for the dynamic case. SEISMIC MOTIONS DUE TO J E R K Y EXTENSION OF A TENSILE CRACK

Suppose that a jerky extension of a tensile crack occurs and the fluid with the excess pressure Ap fills the newly created crack space. We shall calculate the elastic dynamic field caused by this extension. Initially, the crack with length 2l is in equilibrium with the fluid in an excess pressure Ap. Then, an extension of the crack by Al occurs suddenly, as shown in Fig. 1. We calculate the displacement (u, v) from the equilibrium state that exists just before the crack extension. We assume that the fluid in the crack has a finite bulk modulus b. For the crack width d much smaller than seismic wave length, we may assume that the normal stress at the crack surface is equal to the pressure in the liquid. Then, putting the total fluid volume as V, the constitutive equation for liquid A p = b(A V / V ) can be translated into the boundary condition st the crack surface for the elastic motion, as:

268 oyy

= b vfi+d/2)-v(-d/2)

........ d . . . . . . . .

'

-l < x < l

(12)

The effect on seismic motion of viscosity of a fluid layer sandwiched in an elastic solid has been discussed by Aki et al. (1976). For the thickness d of the fluid layer, the shear stress will not significantly be transmitted through the layer if: (13)

d > ~/2~/¢op

where 77 and p are the viscosity and density of fluid respectively, and ¢0 is the angular frequency of seismic motion. For a viscosity of 5 X 103 P, a typical value for a basaltic lava at temperature around 1000°C (Macdonald, 1972), density of 2.8 g/cm 3, and frequency 1 Hz, the inequality (13) becomes d > 25 cm. Since the temperature would be higher for fluid magma in the subterranean channel than for lava observed on the surface, the viscosity would be lower in the channel and the above limit on d may be reduced for our problem. As discussed later, in applying our model to an actual eruption of Kilauea volcano, d is considered to be about 50 cm. Therefore, we neglect the effect of viscosity in this paper. Thus, the boundary condition for shear stress is given by: Oxy = O,

y

= O,

-l < x < l

(14)

To solve the equation of motion for an elastic medium including a crack with the boundary conditions given by (12) and (14), we used the finitedifference scheme for the in-plane elasto-dynamic problem developed by Madariaga (1976) and modified by Fehler and Aki (1977). In all the calculations reported here, Poisson's ratio was assumed to be 0.25. In the finite-difference calculation, seismic motion is generated b y the extension of the crack tip by one grid length. By this extension, the pressure at the grid point located ahead of the tip suddenly becomes P0 + AP. To eliminate the need for calculating the initial static field, we converted the problem to the one in which the time derivative is taken for displacements and stresses everywhere at all instants. This makes all the initial values for the elastic field vanish. In this modified problem, seismic motion is generated by the pressure at the crack tip grid point which is zero initially, raised to ~,P/5"~ over a short time, At, and then set to zero afterwards. The solution of this problem is integrated in time to give the solution of the original problem. We solved the problem for the following four cases: {1) dry crack (b = ~2) dry crack (b = (3} wet crack ( b / p (4} wet crack ( b / ~

0) with both ends extending, 0) with one end extending, • L i d = 5) with both ends extending, • L / d = 5) with one end extending.

The implication of the particular value of b / p . L / d chosen for the we~ crack is discussed later. The normal c o m p o n e n t particle velocity V(x~ t) as a

269

f u n c t i o n o f l o c a t i o n x on t h e c r a c k surface and t i m e t is s h o w n in threedimensional plots in Figs. 2--5. In all t h e cases, the c r a c k is initially ten grid intervals long. Because we use staggered grids f o r stress and particle v e l o c i t y in t h e f i n i t e - d i f f e r e n c e calculation, t h e dispersion is negligible d o w n t o t h e wavelength o f five grid intervals. Thus, we can get a reliable result u p t o t h e wavelength X = L / 2 , or w a v e n u m b e r k = 2~z /X = 4~ / L , w h e r e L is t h e c r a c k length. F o r a small e x t e n s i o n Al o f the crack, the g e n e r a t e d seismic m o t i o n is p r o p o r t i o n a l t o Al. T h e m a x i m u m d i s p l a c e m e n t a m p l i t u d e observed o n t h e c r a c k can be w r i t t e n as: Vmax = c o n s t a n t • A I A P / #

(15)

T h e c o n s t a n t is r o u g h l y e s t i m a t e d b y n u m e r i c a l calculation as 0.2 f o r the case o f s i m u l t a n e o u s e x t e n s i o n at b o t h ends and 0.1 f o r t h e case o f e x t e n s i o n at o n e end. F o r e x a m p l e , if an e x t e n s i o n o f 1 m occurs at an excess pressure A P ~ 20 bars, a m a x i m u m d i s p l a c e m e n t o f 20 ~ 40 # m is e x p e c t e d f o r # = 1011 d y n e s / c m 2.

t

PARTICLE VELOCITY (NORMAL COMPONENT~ AT CRACK SURFACE

~, ~

t

PARTICLE VELOCITY (NORMAL COMPONENT) AT CR

EXTENDED BY ONE GRID L_ENGTHAT BOTH ENDS.

. ~ Y ~ ~

~ ~ ' ~ " ~

TENSILE CRACK EXTENDED BY ONE GRID LENGTH AT RIGHT END.

Fig. 2. N o r m a l - c o m p o n e n t particle v e l o c i t y is p l o t t e d as a f u n c t i o n o f l o c a t i o n x o n t h e c r a c k surface a n d t i m e t f o r t h e case o f a d r y c r a c k w i t h b o t h e n d s e x t e n d i n g b y a grid length. T h e c r a c k is initially 10 grid intervals long.

Fig. 3. The same as Fig. 2, for the case of a dry crack with one end extending by a grid interval.

270

t

PARTICLE VELOCITY

~

(NORMAL COMPONENT) AT

CRACKSURFACE

PARTICLE VELOCITY (NORMAL COMPONENT) AT CRACK SURFACE

~"Y

/Z:

~f' ~ ~

" ' v~

V

FLUID FILLED TENSILE C RACK EXTENDEDBYONE GRIDL E N G T H A T BOTH E N D S

F LUID - FI LL ED TENSILE CRACK E X T E N D E D BY ONE GRID LENGTH AT R I G H T END

~ 7'

Fig, 4. The same as Fig. 2, for t h e case o f a fluid-filled crack (b/t~ . L / d = 5) w i t h b o t h ends e x t e n d i n g b y a grid interval. Fig. 5. The s a m e as Fig. 3, for the case o f a fluid-filled crack ( b / u . L / d = 5) w i t h o n e end e x t e n d i n g b y a grid interval.

FAR-FIELD RADIATION FROM AN EXTENDING TENSILE CRACK

Now, let us turn to the far-field radiation from vibrations caused by the crack extension. We use the well-known equivalent body force theorem (Maruyama, 1963; Burridge and Knopoff, 1964) to find the far-field radiation. The displacement discontinuity normal to the crack plane corresponds to a superposition of a dilatational source and a dipole (two opposite forces at a distance separated along the direction of the force j. For exa~nple, the seismic moment tensor becomes a diagonal matrix with elements XAvS, XAvS and (X + 2#)AvS, where Av is the displacement discontinuity averaged over the crack surface and S is the area of one face of the crack surface. Defining the displacement discontinuity n v(~, t) as a function of time t and location } on the crack surface, one can calculate the amplitude of farfield P-waves, uP( x, t) and S-waves, u S (x, t), in the similar manner as for an earthquake fault (Haskell, 1966; Aki, 1967). Putting the distance between a point on the crack and receiver as r = Ix - ~ [ and the angle between the crack plane normal and the direction to receiver as 0, we obtain: u P ( x , t ) = (4npa3r) -' (X + 2pcos20)

f

At;(~,t--~)dZt

(16) uS(x,t)=(4np~3r)-'psin20

f A6(},t-~-)dZ~ z

271

where a,/~ and p are P-wave velocity, S-wave velocity and density, respectively, and Z represents the crack surface. For simplicity, we shall assume that the shape of the crack is narrow, and take $~ axis long the length of crack. This simplification will not affect our conclusions seriously because we are interested mostly in a relatively low-frequency range. Further, we assume that the distance r0 from crack center to receiver satisfies the inequality:

ro > 2L2 /~ where L is the crack length and X is the wavelength. Then, we can simplify (16) to: L .P(x,t) = (41rp(x3ro)-' (X+2# cos20)W f 46 ( ~,, t - r ° - ~: c°s~ )d$l o L

uS(x,t)=(41rp~:3ro)-'psin2e . W f A6(~,,t

(17)

r o - ~1cos3, )d~l

o

where W is the linear dimension of crack surface in the direction perpendicular to $1, and ~ is the angle between $1 and the direction to receiver. Taking the Fourier transform of (17) with respect to t, we find: luP(x, 60)l = (4~p~3r0) - l ( X + 2 t t c o s 2 0 ) W

f f A~($,, t)exp(-i60t+ i60 cos 7 $,t d $ , d t \

7

(18)

l u S ( x , 60)l = (4~13r0)-1 sin 20 • W

f f AV(~x, t)exp(-i60t+i60 cos______?7~~,)d$, d t We can identify the integral at the right-hand side as the space-time Fourier transform of discontinuity in particle velocity across the crack. Defining:

~6(h, 60) = f f A~(~I, t) exp (-iwt + ik~)d~ldt

(19)

we can write: luP (x, 60)1 = luS(x,w)] = (4~r0)_lsin20.W

60co87

,60

In other words, if we plot [A~(k, w)] in a k-60 diagram, we can look at farfield spectra for all directions at a glance. Figs. 6--9 shows such diagrams for the four cases of crack vibration discussed earlier.

272

There are two surprising features in these diagrams. One is the existence of several spectral peaks. The other is that most of the high-frequency peaks are located outside the range of seismic radiation in the k-co diagram. These features make a strong contrast with the similar diagrams obtained for an earthquake-type rupture (Das, 1976) in which the crack is nucleated at a point and expands to its final size in one continuous rupture propagation. When the reversal of slip direction is inhibited by friction, the m a x i m u m amplitude is at the origin (w = 0, k = 0), decays more or less smoothly for greater l~l and ikl at roughly the same rate at both radiating and nonradiating zones. It is remarkable that these predicted features for a crack vibrating by jerky extension correspond well to the observed features of volcanic tremor. As summarized earlier, two of the three most significant observations are the existence of spectral peaks and the strong high-frequency oscillations in the immediate vicinity of eruption site. The other is the temporal change of peak frequency, which is of course expected if the crack size changes with time. In the case of a dry crack, we find a spectral peak corresponding to the lowest frequency -- fundamental mode vibration -- located in the region of seismic radiation. The frequency is about 0.23 filL in the case when both ends extend, and about 0.15 fl/L in the case when a single end extends. This peak, however, disappears in the case of a fluid-filled crack with b / i t . L / d = 5. Now, let us ask the question what would be the realistic value for the parameter b/it • L i d ? Later, in application to a Hawaiian eruption, we shall

k/2~ CRACKEXTENSIONATBOTHENDS t S/L ~

~....... . .

.

.

"~

~'_

4/L 1

2/L~~I/L

0

"u~]~

~)L

a/L 2fl/L

Fig. 6. A c o n t o u r m a p o f ( 2 0 ~ £ / . ) - '

I~(k,

~)1 defined

by equation (]9)

f o r the ease

of a dry crack with both ends extending. This d i a ~ r i gives the fmr-fie|d radiation in all directions at a glance. The c r o s s - s e c t i o n o f t h e contour along k ffi (-, cos o )/c gives t h e farfield s p e c t r u m for waves w i t h velocity c radiating in a direction e f r o m the rupture propagation direction as shown in equation (20). The region [kS > c~/~ will n o t be observed at the far-field o f S-waves.

I~

#

M

m

m

m

j

I

I

i

I

rl°

-

~1

£

~I~

~I~

L\

"~ \

,\

~I~

%\

"9\

"o\

~

~I~ ~I~ ~

[

-

~I~ ~,~:~ ~f~

-~

~

~

•

m z

m

Z

o

-I

Z

S

z

m

I'n

db

d3

b~

274

5

k

CRACK

FXTENSION

AT

BOTH

ENDS

FLL

D

[

LL C D

B

[ 6 L i~.j ~

x\

o"

8L

k:w/a

c~

/

Ii:

Q~

i

f

/

q'

5 4L

L5

J--:~

16L

/

5 8L

J J

~C

5

•

/

16L

0

7a

7cz 13L

21a 26L

14a 13L

Fig. 8. The same as Fig. 6, for the case o f a fluid-filled crack extending.

35a 26L

21 a 13L

(b/# .Lid ffi 5) w i t h b o t h ends

infer the crack length 1 ~ 2 km and the width about 50 cm; the corresponding Lid is of the order of 103. The bulk modulus b of magma, on the other hand, may be comparable to the rigidity # of rock matrix in the absence of gas bubbles (Murase and McBirney, 1973). The apparent bulk modulus, however, can be drastically reduced by the presence of small amounts of bubbles (Aki et al., 1976). We feel that the value of b/# • L/d -~ 5 is the minimum estimate which may be appropriate for a very shallow crack. Only few spectral peaks are located in the range of seismic radiation for the fluid-filled crack. In the case when both ends extend, we find the spectral peak at frequencies 0.6 ~/L to 0.9 ~/L, depending on the azimuth of radiation. (We assume that a = ~f3-~.) In the case when one end extends, the peaks in the radiation range lie between 0.5 ~/L and 0.7 {UL. Besides, the pattern of spectral peak is very directional in the case when a single end extends. The prominent peak is observed only in the direction of crack propagation.

275 CRACK EXTENSION AT ONE END

kl2.w

t.0~

p.o,

1.~0'

(FLUID-FILLED)

/'°% ~I~'°}'.l~'0~

~-0~

5

2L 55 16L

/

15 8L

/

/

/

25 16L 5 4L

\

/

\.,

\

15 16L 5 8L 5 16L

K

1

a__

/

2L

J~

i

2L

w/2~,

-5 16L

8L -15 16L

4L

/

'i

I

k~.

2.._55 16L -15 8L

16L -35 f< -5_

/ ~/07

/

2L

Fig. 9. The same as Fig. 6, for the case of a fluid-filled crack (b/#.L/d = 5) with one end extending. In t h e case o f a fluid-filled crack, w h e t h e r o n e or b o t h ends are e x t e n d e d , t h e near-field is d o m i n a t e d by a strong p e a k at a f r e q u e n c y 0.9 {3/L. S i n c e f o r a small e x t e n s i o n A1, t h e a m p l i t u d e o f oscillation is p r o p o r t i o n a l t o ~ l, w e can express t h e a b s o l u t e value o f spectral d e n s i t y at a given p e a k (ko, ~ 0) in t e r m s o f m o d e l parameters. In general, t h e spectral d e n s i t y can be w r i t t e n as:

276 IA~(ko, c~o)l ~ C.LoAp.AI/p

(21)

where the value of C for each peak is a dimensionless number which can be found from our numerical results. For a dry crack with both ends extending by Al, C is a b o u t 0.5 at the frequency 0.23 [J/L. For a dry crack with one end extending by Al, C reaches about 0.25 at the frequency 0.15 6/L. As mentioned earlier, these peaks disappear in the case ef a wet crack. For a fluid-filled crack with both ends extending by Al, C is a b o u t 0.15 in the frequency range 0.6 ~/L to 0.9 [UL. For a liquid-filled crack with one end extending b y Al, C is around 0.1 for peaks in the frequency range 0.4 [~/L to 0.8 6/L. Substituting (21) into (20), we find the far-field spectral density for P- and S-waves as: lu P( x,f0)t -~ (4npa3r0) -1 iX + 2U cos 20) C . L . A P . A S P

(22)

luS(x, fo)] ~_ (4n~r0)_ 1 sin20 C . L . A P . A S //

AS = WA l is the increment of crack surface by a single jerky extension at one end. For the case of both ends extending, the actual area increase is 2AS. So far we considered only b o d y waves. Since the tremor source is shallow, it is important to know h o w the contribution of surface waves compares with that of b o d y waves. Let us find the surface wave spectrum generated by a vibrating tensile crack. The solution is well known for a vertically heterogeneous half-space (e.g., Saito, 1967). Taking a cylindrical coordinate (r, 4, z) and Cartesian coordinate ix, y, z), where x = r cos 4 and y = r sin 4, we consider a vertical crack parallel to the x-axis, located at r = 0 at depth z = d. We can write down the vertical c o m p o n e n t displacement of nth mode Rayleigh waves generated by a point force with the m o m e n t tensor corresponding to a tensile crack of the above geometry as follows:

luzn (r, 4, Z; kn, co)l -

r2(o~,kn, Z) C 2 Mxx(~O)knCOS24r,(oJ, kn, d) 4 cUI,

(dr2 t z=d I + Myy(co)knsin:4 r,(~,kn, d) + Mzz(W)\dz]

(23)

with:

Mxx(¢O) = 1._,coXW A~ (w cOS

'¢°)

I().+2.)WA~(¢°c°s4 ,~) c

Myy(a~ ) = ia)

Mzz(~) = --1,~~W A ~ ( ¢°c°s4c ' ~)

(24)

277

where rl(co, kn, z) and r2(co, kn, z) are the z-dependence of the displacement eigen-function of nth mode Rayleigh waves for the horizontal and vertical component, respectively, r~ and r~ are defined in such a way t h a t t h e y are both real and have opposite signs at z = 0. I1 is defined by:

I1 = f p(rl 2 + r~2)dz o

kn is the wavenumber of nth mode, and c is the corresponding phase velocity. The ~-component displacement of the nth mode Love waves may be written as:

lulL(r, z; kn, co) t -

l,(co, k n , z ) V 2 4 cUI1 nkr

X [Mxx(co) -Myy(co)lknsin ~ cos ~ /1(co, kn, d)

(25)

where Mxx(co) and Myy(co) are given in (24), 11(co, kn, z) is the z-dependence of the displacement eigen-function of nth mode Love waves, and:

11 = ~f pl12dz o

Substituting (21) into (23) through (25), we can express the surface wave spectrum in terms of the model parameters. Now, our final step in the model construction is to consider seismic waves from a sequence of numerous jerky extensions. RANDOM EXTENSIONS OF A FLUID-FILLED CRACK AS A MODEL OF CONTINUOUS TREMORS

Let us assume that the jerky extension described above occurs randomly at a rate of n per unit time. Within a time interval A t, nAt extensions are expected, each of which generates the far-field spectrum given in (22), (23), or (25). Since t h e y occur randomly, independently of each other, the Fourier transform G (co) of total seismic motion sampled over the time length A t will be related to the Fourier transform G0(co) of individual signal by: <[G(co)l 2 > = nAtlGo(co) l 2

(26)

where < > stands for averaging over samples of such total motion. On the other hand, the power spectrum P(co) of a random time series is related to the Fourier transform of its finite sample over time length A t by:

P(co) ~-- /2~t

(27)

Thus, we have:

P(co) -~ 7i n[Go(co)[2

(28)

278 For a rough estimate of the root-mean-square (RMS) amplitude of volcanic tremor, we shall assume that the power spectrum at the peak frequency has a box-car shape with height P(f0) and width Af. Then, the RMS amplitude of the volcanic tremor is given approximately by RMS = x/ 2P(fo)Af

(29)

Replacing Go(w) by the P- and S-wave spectra given in (22), and substituting the resultant power spectrum into the above equation, we find: R M S ( u P ) _~ (4~p~3r0)_ , ( X + 2 p c o s 2 0 ) C . L . A P . A S . ~ / n A f P

(30)

R M S ( u s ) ~- (4~13r0)-1 sin 20 C . L . A p . A S x / n A f P

(31)

The Raleigh and Love wave spectra are given respectively by: IRMS(uR)[

r2(co, kn, z) !/ 2 4CUIlcOp V n~-r I XknCOS2~ rl(w, kn, d)

(X + 2p)knsin 2 ~ r, (¢o, kn, d) + X

t C .L'AP . A S x / n ~

(32)

\ dz ]z=d and: RMS(uL) ]

11(co, kn, z) 4 cUI,¢o

2kn sin $ cos ~ l,(co, kn, d)

X C'L.AP.ASx/nAf

(33)

Thus, the far-field observation of tremor amplitude will give a constraint on the crack length L, excess pressure 5 p , frequency n and incremental area AS of a jerky extension. Additional constraints on these parameters are available from: (1) observed frequency of tremor spectral peak, (2) rate of magma transport and eruption, and (3) observed width of dike. In the next section, we shall apply our model to the actual case of eruption in Hawaii. APPLICATION TO THE OCTOBER 5--6, 1963 KILAUEA ERUPTION One of the best described case histories of volcanic tremor is given by Shimozuru et al. (1966) for the October 5--6, 1963 eruption of Kilanea volcano, Hawaii. As shown in Fig. 10, the eruption occurred along a 1 0 - k m section of the central part of the east rift zone between Napau Crater and Kalalua Cone. Napau Crater is located at a distance of a b o u t 15 km from the summit measured along the rift zone. According to Moore and Koyanagi (1969), the episode started with subsidence of Kilauea summit, as shown by

279 155Ol5 ' [

I 9°30 '

155 ° I

SUMMIT NORTH P I T / ~ ~/KILAUEA IKI , o. ~3 "4c"~'(S HALEMAUMAU ". oC~ ~9~ ~ KAO~KI i ' ~ c . ~. ~ (~L~~ K.~LALUA A 7NAPAO FOOT PRINTS

OHALE

/ /

MAKAOPUH I

/

,9°2o'-

PAHOA

o

PACIFIC

. . . . . . . .

Fig. 1 0 . M a p o f t h e c h a i n of craters a n d Kalalua C o n e a l o n g t h e east rift z o n e a n d s o m e o f

the seismograph stations cited in the text. The site of the October 5--6, 1963 eruption is located between Napau Crater and Kalalua Cone. tiltmeter records, at 0306 on October 5, followed by strong harmonic tremor near the summit and in the upper rift zone. The tremor started at 0316, 10 minutes after the initiation of summit subsidence. The eruption presumably (no direct witnesses) started at Napau Crater about 2 hours later, at 0525 when an acoustic shock was heard. This sequence of subsidence-tremoreruption strongly suggests transport of magma from the summit reservoir to the eruptive vents through conduits within the rift zone. Early eruptions were from vents in and near Napau Crater, but by 0930 a line of vents 5--6 km long east-northeast of Napau was fountaining lava 15 m high. Beginning slightly before 1400, new outbreaks occurred 5--10 km east of Napau. The lava fountains there reached as high as 45 m at 1410, but decreased to 15--20 m at 1800. The eruptive activity diminished gradually and all the vents were quiet at 0630 on October 6. The total a m o u n t of erupted magma is estimated as 7 × 106 m 3, while the subsidence at the summit is interpreted as due to 25 x 106 m 3 of magma transferred from the summit reservoir to the rift zone. The harmonic tremor began rather abruptly at all the stations near the summit as well as at Pahoa, located farther east of Kalalua Cone (Fig. 10). Comparison of tremor amplitude at various stations suggests also the underground passage of magma from the summit reservoir progressively eastward to the eruptive vents. The tremor amplitude was greatest near the summit from 0700 to 1300 on October 5, at Makaopuhi (about 3 km west of Napau Crater) from 0700 to 2300 on October 5, and at Pahoa from 1400, October 5, to 0100, October 6. A detailed study of the spatial distribution of tremor amplitude made by the use of a dense array of seismometers delineated the source area of tremor very clearly. Early in the eruption, tremor was most intense near Halemaumau (summit caldera), but after 1400, the region of greatest intensity of tremor was elongated along the upper east rift zone. Unfortunately, the station net-

280

work was limited to the area west of Napau and did n o t cover the area of vents that erupted on October 5--6, 1963. From the evidence described above, it is apparent that the volcanic tremor is caused by the movement of magma through the rift zone conduit from the summit reservoir to the eruptive vents over a distance more than 15 km. As Swanson et al. (1976b) proposed, the movement is probably driven forcefully by an excess pressure in magma rather than a passive flow into already open cracks. This is supported by geodetic data indicating contraction of the south flank in the direction normal to the rift zone during some eruptions in the east rift zone (Swanson et al., 1976b). Thus, it seems reasonable to apply our seismic source model of a growing fluid-driven tensile crack to the October 5--6, 1963 Kilauea eruption. We cannot, however, apply a simple model of a single crack growing from the summit to the eruptive vents for several reasons. First, the pattern of initial increase of tremor amplitude at Makaopuhi was quite similar to the one observed at North Pit. It does not support a single crack front progressing from the summit reservoir. Secondly, a roughly constant peak period of 0.8--0.9 s was observed at distant stations such as Kaoiki, F o o t Prints and Ohale shown in Fig. 10. For a fluid-filled crack, as discussed earlier, we expect the predominant frequency in the far-field to be about 0.5--0.9 ~/L, where is the shear wave velocity and L is the crack length. The shear wave velocity at shallow depth may vary from about 1 to 3 km/s. Taking an average of 2 km/s (corresponding rigidity about 1011 dynes/cm 2 , we find L to be 1--2 km. If we consider the jerky extension of a fluid-filled crack as a model for tremor we must have many separate cracks, each with a length of 1--2 km, rather than a single, long crack extending from the summit to the eruptive vent. A.R. McBirney, reviewing the present paper, noted that the rift zones are not continuous fissures but consist of series of en-echelon gashes (see, e.g., Swanson et al., 1976a). Now let us look at the absolute value of amplitude observed at distant stations mentioned above. From Shimozuru et al. (1966), we find that the tremor amplitude (peak to peak) averaged over the period from 0600 to 1800, October 5, is 0.50, 0.20, and 0.17 pm at Ohale, F o o t Prints and Kaoiki, respectively. Let us first compare relative contributions of body waves and surface waves using equations (30) through (33). The calculation of the surface wave spectrum requires knowledge of the eigen-function. For a rough estimate, we shall assume a homogeneous half-space and consider only Rayleigh waves. The eigen-function for a half-space with Poisson's ratio 0.25 is well known (e.g., Ewing et al., 1957): rl r2

=

e -0"8475kz

0.8475

-

0.5773

e -°'s475kz

-

e --0"3933kz

1.4679

e -°'3933kzl

(34)

where c ffi 0.9194 ~ and I1 = 1.2409 p/k. Since the d o m i n a n t period is about 0.9 s, the corresponding wavelength is

281 a b o u t 1.6 km, assuming a shear velocity of 2 km/s. Substituting the above relations into (32), we obtain: RMS (u R) = ~-~ ~

(p + q sin 2 ~ ) C: L . A p . AS X / n A f

(35)

where k is the wavelength (~ 1.6 km); values o f p and q for different focal depths, h, are given in Table 1. TABLE 1 Values ofp and q for different focal depths, h h (km)

kh

p

q

0.1 0.2 0.5 1.0

0.39 0.79 1.96 3.93

0.090 0.020 -0.061 -0.057

0.144 0.057 -0.050 -0.056

Comparing (35) with (31), we find that: RMS (uS) RMS (u R)

~

sin 20 4 ~ ( p + q sin2 $)

(36)

For the range of focal depth h from 0.1 to 1 km, and epicentral distance from 5 to 10 km, the above formula shows that the amplitude of S-waves is 0.1--0.5 that of Rayleigh waves. Since the amplitude of P-waves is less than that of S-waves by about a factor of ~/a, Rayleigh waves probably are the main waves contributing to the observed volcanic tremor. According to equation (35), then, the amplitude times x/7 must be independent of station. We find the product is 1.1, 0.6 and 0.6 # m km 1/2 for Ohale, F o o t Prints and Kaoiki, respectively. Allowing for a difference by a factor of about 3 between the peak-to-peak observed amplitude and RMS amplitude, we estimate that the average of RMS (uR) • v T i s 0.3 pm km In. Putting this into (35) and using a rough median value f o r p + q sin2~ to be 0.1, we find that: L . A P . / X S x / n A f ~- 1019 dynes cm/s

(37)

where we used the value of C (-~ 0.1) for a fluid-filled crack with one end extended. Since L is determined from the spectral peak to be 1--2 km, and h f is a b o u t 0.1 ~ 0.2 Hz from observed width of the peak, we can rewrite the above equation as: A P . AS x / n = 1 ~ 3 × 1014 dynes/s ~n

(38)

Let us n o w seek other constraints on the parameters given in the above equation. In our model, fluid magma is transported through the ends of all

282 the cracks at a rate of d - A S. n where d is the width of crack. As m e n t i o n e d earlier, during the O ct ober 5--6 eruption, 25 X 106 m 3 of magma was transported f r om the summit t o the rift zone within about 20 hours. The volume o f e r upt ed magma, on the ot he r hand, is estimated as 7 × 106 m 3. Assuming that most of the 25 × 106 m 3 passed through the upper east rift zone, which is the main source of t r e m o r observed at stations operated by Shimozuru et al. (1966), we obtain: d . / k S - n ~ 3 X l 0 s cm3/s

(39)

According t o Swanson et al. (1976b), a typical width for m any observed Hawaiian dikes is 50 cm. Since the width in our formula (39) refers t o the peripheral part of the crack, we should take a value a little less than 50 cm. Putting d ~ 30 cm, say, we obtain a r o u n d n u m b e r for: A S - n -~ 107 cm2/s

(40)

On the o t he r hand, a continuous t r e m o r with frequency near 1 Hz requires that n c a n n o t be much smaller than 1 per second. For an order of magnitude estimate o f model parameters, we shall put n -~ 1 per second. Then, bot h the constraints (38) and (40) are satisfied by the following parameters: AP ~ 20 bars } AS ~ 103 m2 l n -~ 1 per second

(41)

Thus, we n ow have a rather concrete picture of how magma is transported. Jerky extensions of cracks with lengths of 1--2 km are caused by an excess magma pressure of a b o u t 20 bars the total rate of one extension per second. If there are 10 cracks, the extension of an individual crack occurs at a rate of 1 per 10 seconds. The increment of crack area in each j erky extension is a b o u t 103 m 2. If the de pt h of a crack is 1 km, the crack tip advances by 1 m in each jerky extension. The above picture of magma transport does n o t l ook physically unreasonable in terms of quantities involved. However, we find a problem when we calculate the total increment of crack area during the episode. For AS = 103 m 2 , n -~ 1 per second, and the total duration of vigorous activity of volcanic t r e m o r 10 hours, the total crack area must be expanded by 36 km 2 during the episode. This is a considerable increase in crack size which should acco mp an y a significant increase in period of the peak t r e m o r spectrum. The d o m i n a n t period showed an increase (0.79--0.93 second) at Kaioki during the period fr om 0600 to 1800 on O c t o b e r 5, but slightly decreased from 0.84 to 0.73 second at Ohale during t he same period. Although a slight increase in period is observed at a majority of stations, the increase is n o t as large as e x p e c t e d from the model we have described.

283

A CHAIN OF CRACKS CONNECTED

BY NARROW

CHANNELS

To resolve the above inconsistency within the framework of our crack model, we must abandon the idea of a crack steadily expanding by jerky extensions. It is still possible that such expanding cracks generate volcanic tremor, but they did not play a major role in the October 5--6, 1963 eruption. We must consider a more stationary seismic source for the episode. We propose a chain of cracks connected by narrow channels shown schematically in Fig. 11 as the source of volcanic tremor. At the time of the October, 1963 eruption, these cracks might have been already filled with magma remaining from the activity in August, 1963. The channels between neighboring cracks are initially closed. When the pressure in the summit reservoir reaches a certain limit, the first channel opens and magma flows into the first crack. This in turn raises the pressure in the first crack, and the second channel opens letting magma flow into the second crack. Once the flow takes place, the pressure in the first crack will drop and the channel between the first crack and the second will be closed. It will open again when the pressure increases in the first crack due to the magma supply from the reservoir. This process will be repeated at all the channels along the chain of cracks down to the eruption site where magma finally flows to the surface. When the magma with excess pressure AP enters into a channel with area AS, the cracks on both sides of the channel are connected by fluid, and set into vibration as a single crack. The seismic source for this vibration may be modelled by a localized excess pressure AP suddenly applied over an area AS near the center of the double crack. Using the boundary conditions (12) and (14), we calculated the particle velocity at the surface of an in-plane tensile crack when a localized excess pressure is suddenly applied at its center by the same m e t h o d as used in computing the motions for expanding cracks. In the finite-difference calculation, we again used ten grid lengths for the crack and the excess pressure was applied at the two grid points in the center to take advantage of the sym-

SUMMIT ) RE

/ Fig. 11. Magma transport through a chain of cracks connected by narrow channels.

284 metry of the problem. Using (19), we then calculated the frequency wavenumber spectrum A 5 (k, w) of discontinuity in particle velocity across the crack. The results for the dry crack and. the fluid-filled crack with b/t~ • L i d = 5 are shown in Figs 12 and 13, respectively. They are quite similar to the case of crack extension at both ends as shown previously in Fig. 6 and 8, except that the fundamental spectral peak lying along the co-axis is very strong in the case of excitation at the center for both dry and fluid-filled cracks. For other peaks, the t w o cases are very much alike in both amplitude and location in the k-w diagram. In applying the steadily expanding crack model to the Kilauea eruption, we used the case of a fluid-filled crack extended at one end. In this case, the far-field peak occurred in the frequency range 0.4 [3/L to 0.8 M L and the constant C in (21) was around 0.1. In the case of a fluid-filled crack excited at the center, on the other hand, the far-field peaks occur at 0.55 filL and 1.2 {UL. The constant C is about 0.2 for both peaks. Since the excess pressure is applied at two grid points in the latter case, while only one grid point is extended in the former case, the absolute value of peak amplitude for a fixed AS is the same for the two cases. Thus, the tremor amplitude will be the same if the same AP, AS and n are used in both cases. The crack lengths estimated from the observed peak frequency are, however, slightly different. Instead of k/2~T

EXCITATION AT CENTER

5/L

/

\

F '~'-

//

"

"

"\

2L

/

~

i

,i

I'

~P'~"

T-- ~ 1 2 "It

]~\'

....

/ L ~ -

3"/2L

Fig. 12. The same as Fig. 6, for the case o f a dry crack w h e n the e x c e s s pressure is applied at t w o grid points at the center o f the crack.

285

k/2"rr 5/L

EXCITATION AT CENTER (FLUID- FILLED) ~

o,

~'/L

Fig. 13. The same as Fig. 6, for the case of a fluid-filled crack ( b / ~ . L / d = S) when the excess pressure is applied at two grid points at the center of the crack. 1 ~ 2 km obtained before, the double crack length is either 1.2 or 2.6 km. The length of individual cracks of the chain is therefore 0.6 or 1.3 km. DISCUSSION Modifying our starting model of a steadily expanding crack to resolve inconsistencies with observations, we arrived at the revised model of a chain of cracks connected with narrow channels. The length of each crack is around 1 kin, the channel area connecting neighboring cracks is about 103 m 2, and the channel opens with magrnatic excess pressure of about 20 bars. If we consider the total of fifteen cracks along the rift zone from the summit to the Napau crater, the frequency of jerky openings of each channel is about once in 15 seconds. In this preliminary analysis of volcanic tremor, we avoided the treatment of dynamics of gas and fluid inside the crack. Recently, Steinberg and Steinberg (1975) proposed a source model of volcanic tremor in which the motion of gases in the vent plays the major role. They attributed the t r e m o r to selfoscfllations generated by the transition of subsonic to supersonic gas flow, but, unfortunately, they did n o t give quantitative estimates of the tremor period or amplitude. They avoided analyzing the solid part of the tremor generating system. A complete understanding of the dynamics of such a

286

complicated system will take many years. Nevertheless, we feel that the excess magmatic pressure in magma-filled cracks is an important driving force of transport process, at least, at shallow depths. Further comparison with observations and analysis of the fluid motion involved are needed. ACKNOWLEDGEMENTS

Discussion with R.M. Potter, L. A a m o d t and J. Albright of Los Alamos Scientific Laboratories were most helpful to set the direction at the initial stage of this work. The authors greatly benefited from valuable information supplied by R.Y. Koyanagi of U.S. Geological Survey on the volcanic tremor observed in Hawaii. Thanks are also due to B. Chouet for reading the manuscript and making helpful comments. This work was supported by E R D A through contract E (11-1) 2534. REFERENCES Aki, K., 1967. Scaling law of seismic spectrum. J. Geophys. Res., 72: 1217--1231. Aki, K., Chouet, B., Fehler, M., Zandt, G., Koyanagi, R., Colp, J. and Hay, R.G., 1976. Seismic properties of a shallow m a g m a reservoir in Kilauea Iki by active and passive experiments (submitted to J. Geophys. Res.). Barenblatt, G.I., 1962. The mathematical theory of equilibrium cracks in brittle fracture. Adv. Appl. Mech., 7: 56--129. Burridge, R. and Knopoff, L., 1964. Body force equivalents for seismic dislocations. Bull. Seismol. Soc. Am., 54: 1875--1888. Clark, S.P., 1966. Handbook of Physical Constants. Geol. Soc. Am., Mere., 97: 92. Das, S., 1976. A numerical study of rupture propagation and earthquake source mechanism. Sc.D. Thesis, Massachusetts Institute of Technology, Cambridge, Mass. Dibble, R.R., 1974. Volcanic seismology and accompanying activity of Ruapehu volcano, New Zealand. In: L. Civetta, P. Gasparini, G. Luongo and A. Rapotla (Editors), Physical Volcanology. Elsevier, Amsterdam, pp. 49--84. Ewing, W.M., Jardetsky, W.S. and Press, F., 1957. Elastic Waves in Layered Media. McGrawHill, New York, N.Y. Fehler, M. and Aki, K., 1977. Numerical study of diffraction of plane elastic waves by a finite crack with application to location of a magma lens (submitted to Bull. Seismol. Soc. Am.). Haskell, N., 1966. Total energy and energy spectral density of elastic wave radiation from propagating faults. Bull. Seismol. Soc. Am., 56: 125--140. Kakuta, T. and Idegami, Z., 1970. The volcanic tremors of C type at the Sakurazima volcano. Bull. Volcanol. Soc. Japan, 15: 61--74. Kikuchi, S., 1964. On the short-period volcanic micro-tremors at Mt. Aso, H. Bull. Volcanol. Soc. Japan, 9: 9---16. Kubotera, A., 1974. Volcanic tremors at Aso volcano. In: L. Civetta, P. Gasparini, G. Luongo and A. Rapolla (Editors), Physical Volcanology. Elsevier, Amsterdam, pp. 29--48. Macdonald, G.A., 1972. Volcanoes. Prentice-Hall, Englewood Cliffs, N.J., 64 pp. Madariaga, R., 1976. Dynamics of an expanding circular fault. Bull. Seismol. Soc. Am., 66: 63~ G66. Maruyama, T., 1963. On the force equivalents of dynamic elastic dislocations with reference to the earthquake mechanism. Bull. Earthq. Res. Inst., T o k y o Univ., 41: 467--484.

287

Minakami, T., 1974. Seismology of volcanoes in Japan. In: L. Civetta, P. Gasparini, G. Luongo and A. Rapolla (Editors), Physical Volcanology. Elsevier, Amsterdam, pp. 1--28. Moore, J.G. and Koyanagi, R.Y., 1969. The October 1963 eruption of Kilauea Volcano, Hawaii. U.S. Geol. Surv., Prof. Paper 614-C. Murase, T. and McBirney, A.R., 1973. Properties of some common igneous rocks and their melts at high temperature. Bull. Geol. Soc. Am., 84: 3563--3592. Nakamura, K., 1977. Volcanoes as possible indicators of tectonic stress orientation -principle and proposal. J. Volcanol. Geotherm. Res., 2: 1--16. Omer, G.C., 1950. Volcanic tremor. Bull. Seismol. Soc. Am., 40: 175--194. Omori, F., 1912. The eruptions and earthquakes of the Asama-Yama, I, II, III. Bull. Imp. Earthq. Inv. Comm., Tokyo, 6. Potter, R.M. and Dennis, B.R., 1974. Seismic and fluid pressure response from a series of hydraulic fractures in granite. EOS, Trans. Am. Geophys. Union, 55: 431. Saito, M., 1967. Excitation of free oscillations and surface waves by a point source in a vertically heterogeneous earth. J. Geophys. Res., 72: 3689--3699. Sassa, K., 1935. Volcanic micro-tremors and eruption-earthquakes. Mere. Coll. Sci., K y o t o Imp. Univ., Set. A, 18: 255--293. Schick, R. and Riuscetti, M., 1973. An analysis of volcanic tremors at south Italian volcanoes. Z. Geophys., 39: 247--262. Shima, M., 1958. On the second volcanic micro-tremor at volcano Aso. Bull. Disaster Prevent. Inst., K y o t o Univ., 22: 1--6. Shimozuru, D., 1961. Volcanic micro-seisms - - discussion on the origin. Bull. Volcanol. Soc. Japan, 5: 154--162. Shimozuru, D., Goto, K., Nakamuta, O. and Noda, H., 1958. Observations of volcanic micro-tremors at volcano Aso. Bull. Volcanol. Soc. Japan, 3: 35--42. Shimozuru, D., Kamo, K. and Kinoshita, W.T., 1966. Volcanic tremor of Kilauea volcano, Hawaii, during July-December, 1963. Bull. Earthq. Res. Inst., Tokyo Univ., 44: 1093--1133. Steinberg, G.S., 1965. Genesis of volcanic tremors and the long-range forecasting of eruptions. Dokl. Akad. Nauk, 165: 1294--1297. Steinberg, G.S. and Steinberg, A.S., 1975. On possible causes of volcanic tremor. J. Geophys. Res., 11: 1600---1604. Swanson, D.A., Jackson, D.B., Koyanagi, R.Y. and Wright, T.L., 1976a. The February 1969 East Rift eruption of Kilauea volcano, Hawaii. U.S. Geol. Surv., Prof. Paper, 891, 30 pp. Swanson, D.A., Duffield, W.A. and Fiske, R., 1976b. Displacement of the south flank of Kilauea volcano: the result of forceful intrusion of magma into the rift zones. U.S. Geol. Surv., Prof. Paper, 963: 1--39. Wada, T. and Ono, H., 1963. Spectral study of volcanic micro-tremors. Bull. Volcanol. Soc. Japan, 8: 1--10. Watanabe, H., 1961. On the volcanic-tremor of Sakura-jima. Bull. Voleanol. Soc. Japan, 6: 29--41. Weertman, J., 1971a. Theory of water-filled crevasses in glaciers applied to vertical magma transport beneath ocean ridges. J. Geophys. Res., 76: 1171--1183. Weertman, J., 1971b. Velocity at which liquid-filled cracks move in the earth's crust or in glaciers. J. Geophys. Res., 76: 8544--8553.