Jonathan M. Lees

Professor of Geophysics, Seismology, Volcanology

Crustal seismology, seismic imaging, tomography, volcano seismology, geophysical signal analysis

Seismic inversion for three-dimensional variations of velocity and attenuation are often used to delineate magma bodies in the crust and upper mantle. Problems related to spatial resolution and data noise can obscure details relevant to investigating magma chambers, and the introduction of smoothing constraints, or damping, causes blurring. Tomographic inversions for P- and S-wave velocity/attenuation are summarized including large calderas, rift zones and smaller scale subduction zone volcanoes. While results vary considerably from place to place, most anomalies are found to be in the range of ?10% perturbation, a range often controlled by the method of smoothing or regularization imposed during analysis. At many volcanoes high velocity anomalies are observed in the shallow regions below active areas where conduits, dykes or sills are expected to be present. At other locations low velocity perturbations are seen and interpreted as magma accumulation. Resolution limitations and regularization play a significant role in determining the level of perturbation observed in tomographic studies, although there may be regions where diffuse accumulations of magma do not exhibit strong anomalies and their identification will be elusive. Back to Lees Publications

I present several new packages for analyzing seismic data for time series analysis and earthquake focal mechanisms. The packages consists of modules that 1) read in seismic waveform data in various common exchange formats, 2) display data as either event or continuous recordings and 3) performs numerous standard analyses applied to earthquake and volcano monitoring. SeisR is designed as a research tool aimed at investigators who need to quickly assess large amounts of time-series as they are related to the spatial distribution of geologic structure and wave propagation. In addition to time series analysis, a spatial mapping program is included that ties waveforms and radiation patterns to geographical data-base and mapping programs. Lees, J. M., N. Symons, O. Chubarova, V. Gorelchik, and A. Ozerov (2007), Tomographic Images of Kliuchevskoi Volcano P-wave Velocity, in Volcanism and Subduction: The Kamchatka Region, edited by J. Eichelberger, E. Gordeev, M. Kasahara, P. Izbekov and J. M. Lees, pp. 293-302, American Geophysical Union, Washington, D.C.

Three-dimensional structural images of the P-wave velocity below the edifice of the great Klyuchevskoy group of volcanoes in central Kamchatka are derived via tomographic inversion. The structures show a distinct low velocity feature extending from around 20 km depth to 35 km depth, indicating evidence of magma ponding near the Moho discontinuity. The extensive low velocity feature represents, at least to some degree, the source of the large volume of magma currently erupting at the surface near the Klyuchevskoy group. Back to Lees Publications

We present a review of great earthquakes and seismicity patterns along the Alaska-Aleutian and Kamchatka-Kurile arcs as an overview of one of the longest subduction zone complexes on the planet. Seismicity patterns, double seismic zones and focal mechanism solutions are described and used to illustrate the distribution of stress in the Pacific plate as it collides with North America and Eurasia. Seismicity along the Alaska-Aleutian arc is relatively shallow as compared to the Kamchatka-Kurile arc where the plate is considerably older and thicker prior to entering the subduction zone. Tomographic inversions of the slab generally show high velocity anomalies where seismicity is high, presumably tracking the cold subducting lithosphere. Back to Lees Publications

First arrivals of seismic waves were recorded along the Kamchatka arc using broadband seismic stations deployed for one year in 1998-1999. Cross correlation methods were used from a high resolution data set for tomographic inversion of body waves. The P-wave teleseismic tomography shows evidence of slab shoaling along the northern terminus of the Kamchatka subduction zone. Tomographic anomalies corroborate trends in seismicity, geochemistry, heat flow, shear wave splitting, and surface wave inversions. Thermal ablation via contact with asthenosphere, under the proper conditions, is offered as a possible explanation of the observed shoaling of the Kamchatka slab edge. Back to Lees Publications

While wireless sensor networks can generically be used for a wide variety of applications, breakthrough innovations are most often achieved when driven by a genuine need or application, with its specific system-level and science-related requirements and objectives. Hence, our work focuses on the development of wireless sensor network system-on-chip devices and supporting software for volcano monitoring, which we call Sensor Network for Active Volcanoes (SNAV). In this paper we present preliminary results of our research and development work on intelligent sensor networks for monitoring hazardous environments especially the SNAV system-on-chip design for active volcanoes monitoring. Back to Lees Publications

We present a science-centric evaluation of a 19-day sensor network deployment at Reventador, an active volcano in Ecuador. Each of the 16 sensors continuously sampled seismic and acoustic data at 100 Hz. Nodes used an event-detection algorithm to trigger on interesting volcanic activity and initiate reliable data transfer to the base station. During the deployment, the network recorded 229 earthquakes, eruptions, and other seismoacoustic events. The science requirements of reliable data collection, accurate event detection, and high timing precision drive sensor networks in new directions for geophysical monitoring. The main contribution of this paper is an evaluation of the sensor network as a scientific instrument, holding it to the standards of existing instrumentation in terms of data fidelity (the quality and accuracy of the recorded signals) and yield (the quantity of the captured data). We describe an approach to time rectification of the acquired signals that can recover accurate timing despite failures of the underlying time synchronization protocol. In addition, we perform a detailed study of the sensor network?s data through a direct comparison to a standalone data logger, as well as an investigation of seismic and acoustic wave arrival times across the network. Back to Lees Publications

Augmenting heavy and power-hungry data collection equipment with lighter, smaller wireless sensor network nodes leads to faster, larger deployments. Arrays comprising dozens of wireless sensor nodes are now possible, allowing scientific studies that aren?t feasible with traditional instrumentation. Designing sensor networks to support volcanic studies requires addressing the high data rates and high data fidelity these studies demand. The authors? sensor-network application for volcanic data collection relies on triggered event detection and reliable data retrieval to meet bandwidth and data-quality demands. Back to Lees Publications

We studied regional QLg-1 at 1 Hz (Qo-1) around the Korean Peninsula based on broadband, vertical component seismic records of 6 IRIS Global Seismographic Network stations, and 19 Korean stations of the Korea Institute of Geoscience and Mineral Resources. Using 177 seismic events with M between 5.3 and 5.7, and depths less than 50 km, the reversed two station method was applied and 94 high quality interstation paths were selected from 869 possible pairs. These results show high Qo-1 paths around the Sea of Japan (East Sea) reflecting the typical oceanic structure, and low Qo-1 paths around the northeastern China related to inactive seismicity. Assigning these path values into 193 cells around S. Korea with a size of 1 by 1?, we observed that the regional Qo-1 decreased gradually from east to west between 2 and 1 ? 10-3. Back to Lees Publications

A detailed analysis of discrete degassing pulses, chugs, at Sangay volcano, was performed on seismic and infrasonic records to determine the physics of the conduit. Infrasonic chugging signals appear as repetitive pulses with small variations in amplitude and time lag. An automated time-domain analysis was developed to measure with high precision time intervals and amplitudes at different wave arrivals, reducing the possibility error associated with hand picking. Using this automated method, a strong positive correlation of acoustic amplitude with repose time between individual pulses on chugging signals of Sangay was found on numerous oscillating sequences. Frequency gliding of apparent harmonic frequencies generally trends from high to low frequency at Sangay, in contrast to trends at Karymsky Volcano, Russia. A new description of chugging events using wavelet transform methods, appropriate for non-stationary signals, shows subtle changes in the waveforms relate to physical processes in the volcano. A system of nonlinear feed back, based on choked flow at the vent, is postulated as the most likely source of this volcanic tremor. Back to Lees Publications

Reventador Volcano entered an eruptive phase in 2005 which included a wide variety of seismic and infrasonic activity. These are described and illustrated: volcano-tectonic, harmonic tremor, drumbeats, chugging and spasmodic tremor, long period and very long period events. The recording of this simultaneous activity on an array of three broadband, seismo-acoustic instruments provides detailed information of the state of the conduit and vent during this phase of volcanic eruption. Quasi-periodic tremor at Reventador is similar to that observed at other volcanoes and may be used as an indicator of vent aperture. Variations in the vibration modes of the volcano, frequency fluctuations and rapid temporal fluctuations suggest the influx of new material, choking of the vent and possible modification of the conduit geometry during explosions and effusion over a period of six weeks. Back to Lees Publications

In parts of Ecuador's northeastern provinces airwaves are saturated with infrasound (sub-audible acoustic energy < 20 Hz). Here we identify the locations and characterize three distinct sources of local infrasound, including two types of infrasonic sources, which are not commonplace discussed in the literature. The first of these novel sources is an intense and continuous infrasound radiator with a fixed location corresponding to the San Rafael Waterfall. The signal from the river exhibits a tremor-like envelope that is very well correlated across the 3-element infrasound network. Beyond the river, we also observe and map spatially variable sources corresponding to lightning. These transient signals have impulsive onsets, but are not well correlated across the network and are attributable to spatially-distributed source regions. Finally, we identify plentiful infrasound corresponding to Reventador's vent that is associated with unrest at the volcano. This study demonstrates the great utility of dispersed infrasound networks for distinguishing sources and improving interpretation of local types of infrasound radiators.
Back to Lees Publications

Wireless sensor networks have the potential to greatly advance volcano monitoring. Augmenting heavy and power-hungry monitoring equipment with lighter, small wireless sensor network nodes leads to faster, larger deployments. Arrays consisting of dozens of wireless sensor nodes are now feasible, with the additional scale and resolution permitting studies not possible with existing instrumentation. Designing sensor networks for volcanic monitoring requires addressing the high data rates and high data fidelity demanded by this scientific application. We have designed a sensor network application for volcanic monitoring that relies on triggered event detection and reliable data retrieval to meet bandwidth and data quality demands. Here we describe our design and relate our experience deploying our network on Volc´an Reventador, an active volcano in northern Ecuador.
Back to Lees Publications

An experiment on the uses of shear-wave splitting as an imaging tool in fracture-controlled geothermal reservoirs was conducted at Krafla, Iceland. Fifteen days after the beginning of the seismic recording the injection was stopped for eleven days and then restarted, a sequence designed to determine whether shear-wave splitting measurements can detect the transient response of the subsurface crack system to changes in fluid pressure. It was observed that time delays between the fast and slow split shear waves changed significantly and promptly with the stoppage and resumption of injection. Large time delays occurred only during injection, decreased substantially during the stoppage phase, and increased again as injection restarted. Comparisons of these results with similar observations at the Coso geothermal field in California strongly suggest that the time delay of split shear waves can be a useful proxy to monitor fluid pressure in the cracks and changes in crack density.
Back to Lees Publications

An important objective in volcanology is the quantification of eruption intensity through the study of the elastic energy propagated through atmosphere and ground. To this end we deployed a time-synchronized, seismo-acoustic-video installation at Tungurahua Volcano (Ecuador) in November-December 2004 in an attempt to understand the relation between elastic wave radiation and eruptive manifestation. Our results indicate that plume expansion scales very poorly with both the recorded seismic and acoustic trace energy and with the seismic and acoustic signal amplitude. Heightened material flux during Tungurahua eruptions evidently appears not to be representative of elevated source accelerations, which are the primary influences on elastic energy radiation.
Back to Lees Publications

Tungurahua volcano has exhibited 5 eruptive periods since 1999. The last one began in May, 2004, reached its peak in July, and remained with minor bursts until December, 2004. Between June 30 and August 12, 2004, three temporary seismic and infrasonic stations were installed on the southwest, northwest and northeast flanks of the volcano. About 2,000 degassing signals were recorded jointly on high fidelity infrasonic and seismic instrumentation. Recorded signals, classified by waveform character, include: - 1) Explosion events: impulsive, short duration blasts on infrasound and spindle-shaped long-duration seismic signals. Amplitudes of blast signals span three orders of magnitude from 0.1 to 180 Pa (July 21, 03h32 GMT, the largest signal recorded on the closest station). - 2) Roaring events, composed of seismic and infrasound complex signals with broad frequency bands at stations close to the vent. - 3) Chugging events, composed of saw-tooth shaped infrasound signals. Collocated infrasonic and seismic instruments provided a basis for cluster analysis classification of the most conspicuous type of infrasonic signals (blast explosion events). Travel time analysis of seismic first arrivals and infrasonic waves indicates that blast explosions start with a seismic event at a shallow depth (55-218 m), followed ~1 s later by an out-flux of gas, ash and solid material through the vent, creating pressure disturbances with different patterns that do no correlate with size, location or temporal patterns.
Back to Lees Publications

This paper describes our experiences using a wireless sensor network to monitor volcanic eruptions with low-frequency acoustic sensors. We developed a wireless sensor array and deployed it in July 2004 at Volc´an Tungurahua, an active volcano in central Ecuador. The network collected infrasonic (low-frequency acoustic) signals at 102 Hz, transmitting data over a 9 km wireless link to a remote base station. During the deployment, we collected over 54 hours of continuous data which included at least 9 large explosions. Nodes were time-synchronized using a separate GPS receiver, and our data was later correlated with that acquired at a nearby wired sensor array. In addition to continuous sampling, we have developed a distributed event detector that automatically triggers data transmission when a well-correlated signal is received by multiple nodes. We evaluate this approach in terms of reduced energy and bandwidth usage, as well as accuracy of infrasonic signal detection.
Back to Lees Publications

During the summer of 2004 we recorded the seismicity at the Krafla geothermal field for forty days with an array of twenty PASSCAL L-28 4.5-Hz sensors. The Krafla field is located approximately 60 km East of Akureyri in northern Iceland. The array covered an area approximately 5 km N-S by 4 km EW. The field area is located on Holocene lava flows on the Mid-Atlantic Ridge. The array recorded approximately 5 micro-earthquakes per day at a sampling rate of 500 Hz. This high sampling rate is required to exploit newly developed theories on the frequency-dependence of shear-wave splitting (SWS). During the experiment, the injection well was stopped for ten days to study the response of the subsurface crack system to changes in water pressure. SWS is an exploration method based on the analyses of polarizations and time delays of shear waves that have been distorted by the anisotropy of the medium through which the seismic waves have propagated. Epicenters roughly align along the E-W direction, while hypocenters are shallow around the injection well and appear to be related to the on-going injection. Observations of SWS at Krafla have provided evidence for at least two major crack systems oriented approximately N-S and E-W. This last, rather unexpected direction is consistent with results from a simultaneous MT (magneto-telluric) survey. Further SWS study will lead to a more detailed understanding of the fracture locations, sizes, and orientations in the geothermal field.
Back to Lees Publications

A high resolution tomographic study, using cubic B splines parameterization and employing a systematic approach to the choosing of appropriate damping and smoothing parameters, provided a three dimensional P wave velocity map of the Loma Prieta area. 11,977 high quality raypaths from 844 aftershocks of the 1989 Loma Prieta earthquake were used in the inversion. The velocity model exhibits a low velocity feature between the San Andreas and Zayante Vergeles faults in the top 10 km of the crust. This low velocity feature is interpreted as a sedimentary unit exposed to the northwest and separated from the Salinian block by the Zayante Vergeles fault. Below 10 km no consistent change is observed between the Salinian and the Franciscan blocks. There appears to be a high correlation of aftershock activity and localized high velocity anomalies southeast of the Loma Prieta main shock. Whereas this anomaly may represent brittle rocks associated with a fault zone asperity that failed after the main shock, there is evidence to suggest it may be a body of serpentinite. The serpentinite exhibits high velocities and is potentially less competent than surrounding country rock, thus providing a sector along the fault more likely to be associated with many smaller earthquakes or creep behavior.
Back to Lees Publications

Pacific plate subduction at the Aleutian–Kamchatka juncture, or corner, could be accommodated by either a large bend or a tear in the oceanic lithosphere. In this paper, we describe a number of observations which suggest that the Pacific plate terminates abruptly at the Bering transform zone (TZ). Seismicity shoals along the subduction zone from Southern Kamchatka (600 km) to relatively shallow depths near the Kamchatka–Bering Fault intersection (100–200 km). This seismicity shoaling is accompanied by an increase in the heat flow values measured on the Pacific plate. Moreover, unusual volcanic products related to adakites are erupted on Kamchatka peninsula at the juncture. Simple thermal modeling shows that a slab torn and thinner along the northern edge of the Pacific plate would be compatible with the observations. Delayed thickening of the lithosphere due to the Meiji–Hawaiian hotspot may be responsible for the required thinning.
Back to Lees Publications

Eruptive activity at Stromboli Volcano was significantly elevated over background levels in May 2001. During 63 hours of observation, eight vents produced, on average, 17 explosions/hr with an average repose interval of three minutes. During this period the Stromboli vents exhibited consistent seismic and acoustic signatures, based on cross correlation cluster analysis. Dendrogram clustering, based on waveform cross-correlation, was used to illustrate the complexity of the near surface plumbing system of Stromboli's multi-vent edifice. Cross correlations of displacement seismograms produced by explosions at specific craters, such as the Northeast crater (NEC), form dense waveform clusters with correlation coefficients between 0.96-0.99, while displacement waveforms from other craters, such as the Southwest crater (SWC), exhibit loose clusters with correlations between 0.88-0.96. The inconsistency of SWC events, as compared to the NEC, suggests that the vent system there is more heterogeneous. Cluster linkage distance between the NEC cluster and the Central crater (CC) cluster is shorter than the linkage distance between the NEC and SWC clusters, indicating that NEC and CC are more closely related. Infrasonic observations were used to locate vent explosions confirming that the clusters of events are associated with specific vents or craters. Qualitative analyses of acoustic waveforms from approximately 500 explosions reveal that impulsive acoustic signals were associated with short, mechanically simple ground displacement responses. These events may correspond to the bursting of an individual gas slug. Similar degassing mechanisms from vents in the NEC and the CC show common characteristics in their displacement waveforms.
Back to Lees Publications

Explosions of Karymsky volcano often produce signals containing a sequence of repeating pulses recorded on acoustic and seismic sensors, known as chugging. The amplitudes of these pulses correlate with the time interval between pulses. For a given measured acoustic pressure, seismic amplitudes take on arbitrary values up to a specific, empirically determined threshold. Conversely, small seismic amplitude events yielded acoustic waves with large variations and large amplitude seismic events corresponded to large acoustics waves. These observations are not consistent with a source modeled by a resonating conduit. Rather, a model consisting of a sequence of discrete pulses explains the data and provides a framework for understanding the dynamics of degassing at the vent. The physical model for chugging involves a time varying narrowing vent where gasses are released in a series of oscillations which appear to be harmonic but instead are modeled as short term transients, or discrete pulses, suggestive of choked flow.
Back to Lees Publications

Large amplitude, secondary arrivals are modeled as scattering anomalies near the Coso, California, Geothermal field. Polarization and raytracing methods determine the orientation and location of the scattering body. Two models are proposed for the scatterer: 1) a point scatterer located anywhere in a one-dimensional, layered velocity model, 2) a dipping interface between two homogeneous half spaces. Each model is derived by non-linear, grid search inversion for the optimal solution which best predicts observed travel times. In each case the models predict a nearly vertical scatterer southwest of station S4 and Y4, each southeast of Sugarloaf Mountain, a large rhyolite dome. The interface model includes five unknown parameters describing the location and orientation of the interface in addition to the S-wave velocity of the half-space. The S-wave velocity, 3.25 km/s, agrees with independently derived one-dimensional models in this area. The large amplitude, vertical impedance contrast interface coincides with steep gradients of heat flow measured near the surface and with structural boundaries observed in surface geology. The reflector is most probably the sharp boundary between the northern part of the field where there is significant fluid flow versus the southern part where hydrothermal fluids are absent. The interface coincides with geological boundaries and faults recently observed in this region, most likely representing the hydrothermal barrier which channels hot fluids northward.
Back to Lees Publications

A model is proposed to explain temporal patterns of activity in a class of periodically exploding Strombolian-type volcanos. These patterns include major events (explosions) which follow each other every 10-30 minutes and subsequent tremor with a typical period of 1 second. This two-periodic activity is thought to be caused by two distinct mechanisms of accumulation of the elastic energy in the moving magma column: compressibility of the magma in the lower conduit and viscoelastic response of the almost solid magma plug on the top. A release of the elastic energy happens when a stick-slip dynamic phase transition in a boundary layer along the walls of the conduit occurs; this phase transition is driven by the shear stress accumulated in the boundary layer. The first-order character and intrinsic hysteresis of this phase transition explains the long periods of inactivity in the explosion cycle. Temporal characteristics of the model are found to be qualitatively similar to the acoustic and seismic signals recorded at Karymsky volcano in Kamchatka.
Back to Lees Publications

Teleseismic receiver functions (RFs) from a yearlong broadband seismological experiment in Kamchatka reveal regional variations in the Moho, anisotropy in the supra-slab mantle wedge, and, along the eastern coast, Ps converted phases from the steeply dipping slab. We analyze both radial- and transverse-component RFs in bin-averaged epicentral and backazimuthal sweeps, in order to detect Ps moveout and polarity variations diagnostic of interface depth, interface dip, and anisotropic fabric within the shallow mantle and crust. At some stations, the radial RF is overprinted by near-surface resonances, but anisotropic structure can be inferred from the transverse RF. Using forward modeling to match the observed RFs, we find Moho depth to range between 30 and 40 km across the peninsula, with a gradational crust–mantle transition beneath some stations along the eastern coast. Anisotropy beneath the Moho is required to fit the transverse RFs at most stations. Anisotropy in the lower crust is required at a minority of stations. Modeling the amplitude and backazimuthal variation of the Ps waveform suggests that an inclined axis of symmetry and 5–10% anisotropy are typical for the crust and the shallow mantle. The apparent symmetry axes of the anisotropic layers are typically trench-normal, but trench-parallel symmetry axes are found for stations APA and ESS, both at the fringes of the central Kamchatka depression. Transverse RFs from east-coast stations KRO, TUM, ZUP and PET are fit well by two anisotropic mantle layers with trench-normal symmetry axes and opposing tilts. Strong anisotropy in the supraslab mantle wedge suggests that the mantle ‘‘lithosphere’’ beneath the Kamchatka volcanic arc is actively deforming, strained either by wedge corner flow at depth or by trenchward suction of crust as the Pacific slab retreats.
Back to Lees Publications

Back to Lees Publications

In the most basic seismo-acoustic studies at volcanoes, infrasound monitoring enables differentiation between sub-surface seismicity and the seismicity associated with gas release. Under optimal conditions, complicated degassing signals can be understood, relative explosion size can be assessed, and variable seismo-acoustic energy partitioning can be interpreted. The extent to which these points may be investigated depends upon the quality of the infrasonic records (a function of background wind noise, microphone sensitivity, and microphone array geometry) and the type of activity generated by the volcano (frequency of explosions, bandwidth of the signals, and coupling efficiency of the explosion to elastic energy). To illustrate the features, benefits, and limitations of infrasonic recordings at volcanoes, we showcase acoustic and seismic records from five volcanoes characterized by explosive degassing. These five volcanoes (Erebus in Antarctica, Karymsky in Russia, and Sangay, Tungurahua, and Pichincha in Ecuador) were the focus of seismo-acoustic experiments between 1997 and 2000. Each case study provides background information about the volcanic activity, an overview of visual observations during the period of monitoring, and examples of seismo^acoustic data. We discuss the benefits and utility of the infrasound study at each respective volcano. Finally, we compare the infrasound records and eruptive activity from these volcanoes with other volcanoes that have been the focus of previous seismo-acoustic experiments.
Back to Lees Publications

This paper reviews geophysical and seismological imaging in the Coso geothermal field, located in east central California. The Coso geothermal production area covers an area approximately 6 10 km2. Although regional seismicity is addressed, as it sheds light on the magma, or heat, sources in the field, the main focus of this paper is on the main production area. Three-dimensional inversions for P- and S-wave velocity variations, distribution of attenuation, and anisotropy are presented side-by-side so that anomalies can be compared spatially in a direct manner. Velocity inversions for P- and S-waves are combined for direct determination of Poisson's ratio and indirect estimation of variations of porosity in the field. Anomalies southeast of Sugarloaf Mountain are prominent on nearly all analyses. The anomalies coincide with high levels of seismicity and with stress anomalies as determined from earthquake focal mechanism analysis and seismic anisotropy distribution. They also correlate with high heat flow in the field and the termination of geothermal production to the south. I speculate that an intrusion is present in this region that causes significant perturbation of stress in the field.
Back to Lees Publications

The temporal and spatial distribution of seismicity in the Coso and Indian Wells valley region of eastern California are discussed in this study. An analysis of fault-related seismicity in the region let us conclude that the Little Lake fault and the Airport Lake faults are the most significant seismogenic zones. The faulting pattern clearly demarcates the region as a transition between the San Andreas type strike-slip regime to the west and Basin and Range extension to the east. We present the spatial and temporal variations in seismicity immediately following significant earthquakes in nearby regions over the last fifteen years with special emphasis on larger earthquakes (M 5) over the last five years. The Ridgecrest earthquakes of 1995 show a complicated faulting pattern as the rupture changes from normal-slip to right-slip at depth. The inter-relationships between the Coso earthquakes of 1996 and 1998 are presented as a set of conjugate events. Analysis of earthquake source mechanisms shows evidence for lateral variations of faulting pattern in eastern California. Earthquake focal mechanisms are used to estimate local stress orientation within the Coso geothermal field. We have identified a boundary between a transpressional regime and a transtensional regime inside the field which correlates with observed spatial variations of heat flow and seismic attenuation, velocity and anisotropy.
Back to Lees Publications

New PASSCAL data have been acquired along the extent of the Kamchatka Peninsula to examine the interaction of the Pacific Plate and the mantle in the corner junction of the Aleutian and Kamchatka trenches. The project, called the Side Edge of Kamchatka Slab, is a collaborative effort between researchers at Yale University and the Russian Academy of Sciences Institutes in Petropavlovsk-Kamchatski, in particular the Institute of Volcanology (Alexei Ozerov) and KOMSP (Evgenii Gordeev). In Russia, Kamchatka is known as the caboose, the last car on the train. It is practically the farthest one can get from Moscow, the center of cultural life in Russia. Being nine time zones away from the financial centers creates a sense of isolation and liberation. For Russians, Kamchatka is a land of dreams and possibilities, much as Alaska is the last frontier for Americans. This gives one the sense that being in Kamchatka is like being on the edge of the world. But the feeling has more significance than the mere waywardness of the place. The majestic volcanoes, the continuous seismicity and the incredible geology are omnipresent reminders that this is a special place. The relentless subduction of the Pacific plate plunges beneath Kamchatka, 70mm per year. Kamchatka is home to the most magnificent volcanoes in the Pacific Rim (28 active volcanoes) and neighbor to the Bering strike-slip fault, marking the western end of the Aleutian-Komandorsky Islands. This land mass provides an exceptional platform for investigating the interactions of volcanism, tectonics, and mantle dynamics.
Back to Lees Publications

The junction of the Aleutian Island and the Kamchatka peninsula defines a sharp turn in the boundary of the Pacific and North American plates, terminating the subduction zones of the northwest Pacific. The regional pattern of shear-wave birefringence near the junction indicates that trench-parallel strain follows the seismogenic Benioff zone, but rotates to trench-normal elsewhere. Asthenospheric mantle is inferred to flow around and beneath the disrupted slab edge, and may influence the shallowing dip of the Benioff zone at the Aleutian junction.
Back to Lees Publications

Data are presented to show that slab melt-related volcanism (adakitic volcanism) is occurring in the Aleutian-Kamchatka junction, and in the western Aleutians where the plate boundary switches from convergent to strike-slip motion. The subducting lithosphere in these areas is relatively old (>40Ma), and is therefore not expected to melt during subduction (i.e., Defant and Drummond, 1990). It is argued that the slab is melting in these areas because the side edge of the down-going plate is warmed by asthenospheric flow on three sides. The absence of deep or intermediate-depth seismicity in these areas supports a warm slab interpretation. An anomalously large slab melt contribution to arc volcanic rocks in the Aleutian-Kamchatka junction and in the western Aleutians is interpreted to be a geochemical edge effect. The presence of this effect in these areas is consistent with the idea that the Pacific plate is being torn in an unzipping motion as it moves westward from the convergent to the strike-slip portions of the Aleutian arc.
Back to Lees Publications

A new program for exploratory data analysis in three and four dimensions is presented and described. Interactive communication between diverse datasets is stressed as the main gestalt of the Geotouch program. The primary kinds of data include point, vector, raster and wireframe data sets, in addition to specialized forms such as focal mechanisms and ellipsoidal information. The software includes methods for cutting cross sections at arbitrary angles, spinning objects in three space and animating time series of punctual data, such as hypocenter series and volcano eruptions. The program is written in POSIX compliant C using X-windows for Unix systems and has been ported to Linux. Free access, via the Internet, to executable binary, source code and documentation make this package an attractive alternative to expensive or unwieldy commercial options.
Back to Lees Publications

Frequent degassing explosions, occuring at intervals of minutes to tens of minutes, are common at many active basaltic and andesitic volcanoes worldwide. In August 1997, April 1998, and September 1998 we recorded seismic and acoustic signals generated at two andesitic volcanoes with 'Strombolian-type' activity. Despite variations in explosion frequency (5 to 15 per hour at Karymsky as opposed to 1 to 3 per hour at Sangay), the signatures of the explosions are remarkably similar at these two, diverse field sites. In all explosions, gas emission begins rapidly and is correlated with an impulsive acoustic pressure pulse. Seismic waveforms are extremely emergent and begin 1 to 2 seconds before the explosion itself. We classify explosion events at the two volcanoes as either short-duration (less than 1 minute) simple impulses or long-duration (up to 5 minutes) tremor events. Many tremor events have harmonic frequency spectra and correspond to regular 1 second acoustic pulses, often audible, that sound like chugging from a locomotive. Chugging events are intermittent, suggesting that the geometry or geochemistry of the process is variable over short time scales. We attribute the 1 Hz periodic chugs to a resonant phenomena in the upper section of conduit.
Back to Lees Publications

High resolution three-dimensional compressional and shear wave velocity models, derived from microearthquake traveltimes, are used to map the distribution of Poisson's ratio and porosity at Coso Geothermal Area, Inyo County, California. Spatial resolution of the three-dimensional Poisson's ratio and porosity distributions is estimated to be 0.5 km horizontally and 0.8 km vertically. Model uncertainties, \pm 1% in the interior and \pm 2.3\% around the edge of the model, are estimated by a jackknife method. We use perturbations of r = V_p / V_s ratio and \Psi = V_p * V_s product to derive distributions of Poisson's ratio, \sigma, and porosity, which are then used to constrain and delineate possible zones of intense heat, fracture accumulation and fluid saturation. Poisson's ratio at Coso ranges from 0.15 to 0.35 with an average of 0.224, lower than the crustal average of 0.25. High Poisson's ratios are more extensive in shallower depths (<.5 km) while lower Poisson's ratios are found in the deeper section (1.5--3.0 km) of the target area. Two major features with low Poisson's ratio are identified at geothermal production depth (1--3 km) around stations S2-S6 and S1-S3-S4. The two low \sigma features are separated by a northwest-southeast trending high \sigma belt with variable width of 1 ~ 3 km. A high \Psi body is found around S2 and S6, and extends down in depth. A circular, low \Psi belt corresponding to the high \sigma belt, is located around S2-S6 and is linked to a previously reported structure in V_s tomography. This low \Psi (highly porous) belt is probably a horizontal conduit/reservoir of geothermal fluid. A vertical, low \sigma and high \Psi channel beneath triangle S1-S3-S4 corresponds to a high attenuation, dome-like feature. We propose an upwell-and-spread magma intrusion model for the last major magmatism in the Coso region. The magmatic upwelling is centered in the S1-S3-S4 area. The model predicts potential geothermal resources to the south and west of triangle S1-S3-S4 based on local faulting patterns.
Back to Lees Publications

A new inversion method for P-wave anisotropy [{\em Wu and Lees}, 1999a] has been applied to high-precision, microseismic travel-time data collected at Coso geothermal region, California. Direction-dependent P-wave velocity, and thus its perturbation, are represented by a symmetric positive definite matrix A instead of a scalar. The resulting anisotropy distribution is used to estimate variations in crack density, stress distribution and permeability within the producing geothermal field. A circular dome-like structure is observed at the southwestern part of the geothermal region southwest of Sugarloaf Mountain. Using a linear stress-bulk modulus relation and deviatoric stress is estimated to be 3-6 MPa at geothermal production depths (1-2 km), assuming all the anisotropy is realted to stress. The stress field is compressional NNE-SSW and dilational WNW-ESE, coinciding with a previous, independent study using earthquake focal mechanisms. Following a theory on flat, elliptic cracks, residual crack density estimated from P anisotropy is about 0.0078, assuming crack aspect ratios >>1:60 and is ~0.041 when crack aspect ratios are close to 1:60. Residual crack orientation distribution is related to velocity anisotropy. Based on the anisotropic part of crack density distribution function, the anisotropic part of permeability distribution may be calculated by a statistical approach via simple parallel fluid flow along cracks.
Back to Lees Publications

Two recent earthquake sequences near the Coso geothermal field show clear evidence of faulting along conjugate planes. We present results from analyzing an earthquake sequence occurring in 1998, and compare it with a similar sequence that occurred in 1996. The two sequences followed mainshocks which occurred on November 27, 1996 and March 6, 1998. Both mainshocks ruptured approximately co-located regions of the same fault system. Following a comparison with the background seismicity of the Coso region, we have detected evidence of stress loading within the geothermal field that appears to be in response to the 1998 earthquakes. The M_L = 5.2 mainshock in the 1998 sequence occurred at 5:47 am UTC, and waslocated approximately 45 km north of the town of Ridgecrest in the Coso range. The mainshock of the 1996 sequence had a M_L magnitude of 5.3. There have been no observable surface ruptures associated with either of these sequences. Though the mainshocks for both sequences were located about 900 m apart and have nearly the same local magnitudes, the sequences differ in both their temporal and spatial characteristics. An analysis of the fault plane solutions of the mainshocks and the aftershock locations suggests that the two sequences ruptured fault planes which are perpendicular to one another. We observe a much faster temporal decay of the 1998 sequence compared to the one in 1996; moreover, while the 1996 sequence was not followed by any sizeable (i.e., M_L > 4.0) aftershocks, the 1998 sequence had four such events. From an estimate of the tectonic stressing rate on the fault that produced the 1998 sequence, we infer a repeat cycle of 135 years for an earthquake of comparable magnitude at Coso.
Back to Lees Publications

We use a multiple-empirical Green's function method to determine source properties of small (M-0.4 to 1.3) earthquakes and P and S-wave attenuation at the Coso Geothermal field, California. Source properties of a previously-identified set of clustered events from the Coso geothermal region are first analyzed using an empirical Green's function (EGF) method. Stress drops values of at least 0.5-1 MPa are inferred for all of the events; in many cases, the corner frequency is outside the usable bandwidth and can only be constrained as being higher than ~3 MPa. P- and S-wave stress drop estimates are identical to the resolution limits of the data. These results are indistinguishable from numerous EGF studies of M2-5 earthquakes, suggesting a similarity in rupture processes that extends to events that are both tiny and induced and providing further support for Byerlee's Law. Whole-path Q estimates for P and S waves are determined using the multiple-empirical Green's function (MEGF) method of {\it Hough}, (1997), whereby spectra from clusters of colocated events at a given station are inverted for a single attenuation parameter, \kappa, with source parameters constrained from EGF analysis. The \kappa estimates, which we infer to be resolved to within 0.003 sec or better, exhibit almost as much scatter as a function of hypocentral distance as do values from previous single-spectrum studies for which much higher uncertainties in individual \kappa estimates are expected. The variability in \kappa estimates determined here therefore suggests real lateral variability in Q structure. Although the raypath coverage is too sparse to yield a complete three-dimensional attenuation tomographic image, we invert the inferred \kappa value for three-dimensional structure using a damped least-squares method and the results do reveal significant lateral variability in Q structure. The inferred attenuation variability corresponds to the heat-flow variations within the geothermal region. A central low-Q region corresponds well with the central high-heat flow region; additional detailed structure is also suggested.
Back to Lees Publications

We present results from a local earthquake tomographic imaging experiment in the greater Mount Rainier area. We inverted P-wave arrival times from local earthquakes recorded at permanent and temporary Pacific Northwest Seismograph Network seismographs between 1980 and 1996. We used a method similar to that described by Lees and Crosson [1989], modified to incorporate the parameter separation method for decoupling the hypocenter and velocity problems. In the upper 7 km of the resulting model there is good correlation between velocity anomalies and surface geology. Many focal mechanisms within the St. Helens seismic zone have nodal planes parallel to the epicentral trend as well as to a north-south trending low-velocity trough, leading us speculate that the trough represents a zone of structural weakness in which a moderate (M 6.5-7.0) earthquake could occur. In contrast, the western Rainier seismic zone (WRSZ) does not correlate in any simple way with anomaly patterns or focal mechanism fault planes, leading us to infer that the WRSZ is less likely to experience a moderate earthquake. A ~10 km-wide low-velocity anomaly occurs 5 to 18 km beneath the summit of Mount Rainier, which we interpret to be a signal of a region composed of hot, fractured rock with possible small amounts of melt or fluid. No systematic velocity pattern is observed in association with the southern Washington Cascades conductor. A mid-crustal anomaly parallels the Olympic-Wallowa lineament (OWL) as well as several other geophysical trends, indicating that the OWL may play an important role in regional tectonics.
Back to Lees Publications

High precision P and S wave travel times for 2104 microearthquakes with focus <6 km are used in a non-linear inversion to derive high-resolution three-dimensional compressional and shear velocity structures at the Coso Geothermal Area in eastern California. Block size for the inversion is 0.2 km horizontally and 0.5 km vertically and inversions are investigated in the upper 5 km of the geothermal area. Spatial resolution, calculated by synthetic modeling of a cross model at critical locations, is estimated to be 0.35 km for Vp and 0.5 km for Vs. Model uncertainties are estimated by a jackknife approach and simulation of random and associated picking errors. Low-velocity zones for both P and S waves are identified at geothermal production depths (1-3 km). A large, low Vp (-6%) zone is found at depth 2-2.5 km 2 km southwest of Sugarloaf Mountain where high attenuation has been previously reported. However, a general high-Vp zone is seen under Coso Hot Springs with a slightly low Vs zone, which is characteristic of fluid saturation. The overall distributions of Vp and Vs perturbations do not correlate. An isolated high-Vs (+9%) feature, about 2 km in diameter, is unambiguously seen 2 km due west of Sugarloaf extending from surface to depth. This feature is surrounded by a circular, low-Vs belt of ~1 km width. The surrounding belt is probably the cracked, high-porosity reservoir/conduit of geothermal fluid flow. In the 2 km southwest Sugarloaf region, we found low Vp and high Vs at geothermal production depths from 1 to 2.5 km. Combined with attenuation results, this may represent a hot, fluid-depleted center of magmatic activity.
Back to Lees Publications

A new method for inverting P-wave travel times for seismic anisotropy on a local scale is presented and tested. In this analysis, direction-dependent seismic velocity is represented by a second- or fourth-order Cartesian tensor, which is shown to be equivalent to decomposing a velocity surface using a basis set of Cartesian products of unit vectors. The new inversion method for P- and S-wave anisotropy from traveltime data is based on the tensor decomposition. The formulation is formally derived from a Taylor series expansion of a continuously extended, three-dimensional velocity function originally defined on the surface of the unit sphere. This approach allows us to solve a linear inversion instead of the standard nonlinear method. The resultant, linearized, fourth-order traveltime equation is similar to a previous, fourth-order result (Chapman and Pratt, 1992) although our representation offers a natural second-order simplification. Conventional isotropic traveltime tomography is a special case of our tensorial representation of velocities. P-wave velocity can be represented by a second-order tensor(matrix) as a first approximation, although S-wave traveltime tomography is intrinsically fourth order because of S-wave solution duality. Differences of isotropic and anisotropic parameterizations are investigated when velocity is represented by a matrix A.

The tradeoff between isotropy and anisotropy in practical tomography, which differs from the fundamental deficiency of anisotropic traveltime tomography (Mochizuki, 1997), is shown to be ~1, i.e., their effects are of the same order. We conclude that anisotropic considerations may be important in velocity inversions where ray coverage is less than optimal. On the other hand, when the ray directional coverage is complete and balanced, effects of anisotropy sum to zero and the isotropic part gives the result obtained from inverting for isotropic variations of velocity alone. Synthetic test datasets are inverted, demonstrating the effectiveness of the new inversion approach. When ray coverage is fairly complete, original anisotropy is well recovered, even with random noise introduced, although anisotropy ambiguities arise where ray coverage is limited. Random noise was found to be less important than ray directional coverage in anisotropic inversions.
Back to Lees Publications

During the summer of 1997, Karymsky Volcano produced summit explosions about six times each hour. Typical explosive episodes lasted between 30 seconds and three minutes, produced gas and ash columns several hundred meters high, and ejected some incandescent material. To better understand the physical source mechanisms responsible, we recorded hundreds of explosions with a three component broad-band seismometer and an infrasonic pressure sensor located 1650 meters from the active vent. Nearly every explosion is recorded as an emergent yet identical seismic wavelet which is followed 4.15 s later by an impulsive acoustic arrival. We interpret the signals as a near-surface gas volume burst which fractures the vent `plug,' lowers the lithostatic pressure within the magma column, and often induces further degassing. When degassing continues, it is generally manifested as either a series of regular one second `chugging' explosions, steady higher frequency `jetting', or a hybrid combination. We believe that the seismic signature for `chugs,' short duration harmonic tremor with integer overtones, is the result of repeated gas volume bursts at the vent. In contrast, seismograms for jetting are non-harmonic and contain higher frequencies. We believe that the competing degassing behaviors are influenced by the gas flux as well as the plug/conduit characteristics. We propose that a plug exists due to a viscosity gradient caused by volatile depletion in the upper conduit.
Back to Lees Publications

We have searched the Coso geothermal field for microseismicity in seismic doublets, co-located hypocenters that appear to have nearly identical waveforms. Using 1085 high quality events from 1993-1994, we identified numerous doublets, some occurring within minutes of each other. We subdivided hypocentral data into spatial clusters to reduce the computational burden and evaluated multiple cross correlation pairs, assigning scores to each pair. As an example, one spatial cluster includes 183 events yielding 96 high correlation (>0.6) paired events. To isolate potential multiplets, equivalence class analysis and cluster analysis routines were used. Among the 96 high correlation pairs 24 equivalence classes have been isolated. While most of these are doublets, 8 classes include 3 or more cluster members and one class include 16 members. Relative locations were calculated using phase shifts between corresponding events. Detailed analysis of hypocenter relocations shows elongate, vertical structure with apparent random temporal variations. The multiplets do not appear to be true repeating events; rather they are clusters of small, nearly identically oriented ruptures, perhaps representing swarms of fractures activated by fluid pressure fluctuations. Using the small volumes encompassing each multiplet, we estimate fracture densities measure between 0.02 to 0.4 1/m and are largest near injection wells.
Back to Lees Publications

Microseismicity, stress, and fracture in the Coso geothermal field are investigated using seismicity, focal mechanisms and stress analysis. Comparison of hypocenters of microearthquakes with locations of development wells indicates that microseismic activity has increased since the commencement of fluid injection and circulation. Microearthquakes in the geothermal field are proposed as indicators of shear fracturing associated with fluid injection and circulation along major pre-existing fractures. High seismicity zones are associated with the main fluid flow paths within the geothermal system. Calculated stress patterns from focal mechanisms provide direct evidence for the boundary between significantly different stress regimes within the Coso geothermal field.

Microseismicity in the Coso geothermal field is spatially but not temporally related to regional seismicity extending southeast of the field. The spatial distribution of these events defines a northwest-trending seismic-fracture zone, consistent with a previously defined northwest-striking zone. The abrupt decrease of seismicity below this fracture zone may provide seismic evidence for the existence of a brittle and ductile transition zone at shallow depth beneath the Coso geothermal field.
Back to Lees Publications

Coherency between seismic signals recorded in tight arrays is investigated by reorganizing data via waveform cluster analysis. Multitaper spectral estimates and coherency functions suggest that pair-wise coherency between stations located only a few meters apart at Pinyon Flat, CA, is low above approximately 15 Hz. This observation is examined in detail by searching for clusters of similar waveforms using single link cluster analysis on suites of stations located in the 1990 high frequency, 60 station Pinyon Flat, CA, deployment. The analysis requires that noisy line spectra, introduced by instrument noise and shallow resonances be removed prior to analysis. Spectrum reshaping is applied to remove biases associated with smearing periodic signals in the frequency domain. Cluster analysis is performed on coherency and correlation scores derived from multi-taper analysis of 60 signals in the surface and borehole arrays. Cluster analysis shows that while simple pair-wise comparisons of coherency suggests that local scattering dominates the spectrum above 15 Hz, spatial clustering indicates that there is more structure in arriving waveforms than is apparent in the simple analysis. Stations that are clustered together spatially also cluster in coherency score, at frequencies as high as 50 Hz. By re-arranging seismograms according to associated clusters, signals along north-south and east-west arms of the array are separable at frequencies ranging from 5-50 Hz. These observations have important implications for the use of high frequency seismic spectra in determination of source and path effects for small magnitude, local earthquake analyses.
Back to Lees Publications

A high resolution tomographic study, using cubic B-splines parameterization and employing a systematic approach to the choosing of appropriate damping and smoothing parameters, provided a three-dimensional P-wave velocity map of the Loma Prieta area. 11,977 high quality raypaths from 844 aftershocks of the 1989 Loma Prieta earthquake were used in the inversion. The velocity model exhibits a low-velocity feature between the San Andreas and Zayante-Vergeles faults in the top 10 km of the crust. This low-velocity feature is interpreted as a sedimentary unit exposed to the northwest and separated from the Salinian block by the Zayante-Vergeles fault. Below 10 km no consistent change is observed between the Salinian and the Franciscan blocks. There appears to be a high correlation of aftershock activity and localized high-velocity anomalies southeast of the Loma Prieta main shock. Whereas this anomaly may represent brittle rocks associated with a fault zone asperity that failed after the main shock, there is evidence to suggest it may be a body of serpentinite. The serpentinite exhibits high velocities and is potentially less competent than surrounding country rock, thus providing a sector along the fault more likely to be associated with many smaller earthquakes or creep behavior.
Back to Lees Publications

This study investigates the use of a station influence statistic to identify velocity model shortcomings in the earthquake hypocenter location problem. Two groups of microearthquake events are examined. The first is a group of 81 events from the Mount St. Helens region which occurred between November 1987 and September 1991; the second, 110 well located events from the 1992 Joshua Tree aftershock sequence. We describe a method for validating a postulated earth model. Let $lambda$ denote the hypocenter estimates that Geiger's method obtains. Systematically remove each station observation from the location problem and recompute the location estimate. Call this estimate $lambda$(i) when the i-th station is removed. For a single event define a station's influence (SI) as a weighted difference between lambda and $lambda$(i). Distributional summaries of SI statistics across events are used to identify model shortcomings: given a specified velocity model, SI distributions which are not homogeneous across stations provide evidence of model inadequacies and/or failures in the weighting scheme. We show that velocity model shortcomings detected using SI statistics for the Mount St. Helens sequence under a one-dimensional model appear to correlate with known physical anomalies; while SI distributions evaluated under a 3-dimensional model are more homogeneous and reflect a modest improvement over the 1-dimensional model. SI distributions provide evidence of model failure for the Joshua Tree sequence under a 1-dimensional model, but no evidence of failure under a 3-dimensional model. Finally, the weighting scheme's validity is verified for the Joshua Tree sequence under the 3-dimensional model.
Back to Lees Publications

Pulse width data are used to invert for attenuation structure in the Coso geothermal area, California. The dataset consists of pulse width measurements of 838 microseismic events recorded on a seismic array of 16 downhole stations between August 1993 and March 1994. The quality factor Q correlates well with surface geology and surface heat flow observations. A broad region of low Q(~30-37) is located at 0.5-1.2 km depth below Devil's Kitchen, Nicol Prospects and Coso Hot Springs. A vertical, low Q(~36) in contrast with surrounding rock of 80) region is interpreted as a channel through which hydrothermal energy is transported from depth to the surface. The location of the channel is between stations S1 and S4 and its dimension is about 1 km. At the deep end of the channel, a large, broad body of low Q is also located at 3 km depth 2-4 km to the southwest of Nicol Prospects and Devil's Kitchen. Since it lies at the bottom of the target region and beyond the scope of seismicity, further research is needed to constrain its extent. Numerical modeling with a pseudospectral method is also done to investigate the applicability of the inversion scheme to fractured regions. A linear relationship between pulse width broadening and travel time is upheld, and the proportional constants are estimated.
Back to Lees Publications

A new method for calculating boundary conditions at the free surface and along absorbing boundaries of a finite grid is presented. A finite, twice differentiable reduction function which achieves a 99% reduction over 3 wavelengths is proposed and tested. In the context of pseudospectral wave propagation, this implies a boundary layer of at least 6 grid nodes. The method is analyzed in one and two dimensions and the problems of waves impinging on corners are addressed. The reduction function recommended is gamma_R = \alpha (1+cos(pi x))^2 where \alpha is a parameter to be determined by optimization. Tests of the performance of the new method versus other common schemes are presented and analyzed. We provide a strategy for determining the optimal parameter in the reduction function. Synthetic Rayleigh waves are observed at the free surface of the simulation. Experiments with a vertical fault plane show the presence of direct, reflected, transmitted and head waves. The presence of head waves may be used to analyze velocity contrasts across fault zones.
Back to Lees Publications

Three-dimensional Q-1 variations in the aftershock region of Loma Prieta are derived by tomographic inversion. The data set consists of over 4000 aftershock recordings at 22 PASSCAL (Program for Array Seismic Studies of the Continental Lithosphere) stations deployed after the Loma Prieta mainshock of 1989. Estimates of attenuation are determined from nonlinear least squares best fits to the Fourier amplitude spectrum of P and S wave arrivals. The linear attenuation inversion is accomplished by using three-dimensional velocity variations derived previously in nonlinear velocity inversions. Low Q is observed near the surface and Q generally increases with depth. The southwest side of the San Andreas fault exhibits lower Q than does the northeast side and this feature apparently extends to approximately 7 km depth. The fault zone, as determined by the dipping plane of aftershock activity, is characterized by slightly higher Qp and lower Qs, compared to regions immediately adjacent to the fault. These correlate with high- velocity anomalies associated with seismicity at depth. The results are in agreement with earlier observations regarding the association of high-velocity anomalies, seismicity, and fault zone asperities.
Back to Lees Publications

Tomographic inversions for velocity variations in western Washington indicate a high correlation with surface geology and geophysical measurements, including gravity observations. By assuming a simple linear relationship between density and velocity (Birch's law) it is possible to calculate the gravity field predicted from the velocity perturbations obtained by local tomographic inversion. While the predicted gravity matches observations in parts of the model the overall correlation is not satisfactory. In this paper we suggest a method of constraining the tomographic inversion to fit the gravity observations simultaneously with the seismic travel time data. The method is shown to work well with synthetic data in 3 dimensions where the assumption of Birch's Law holds strictly. If the sources of the gravity anomalies are assumed to be spatially localized, integration can be carried out over a relatively small volume below the observation points and sparse matrix techniques can be applied. We have applied the constrained inversion method to western Washington using 4,387 shallow earthquakes, to depths of 40.0 km, (36,865 raypaths) covering a 150´250 km region and found that the gravitational constraints may be satisfied with minor effect on the degree of misfit to the seismic data.
Back to Lees Publications

Tomographic inversion of P-wave arrival times from aftershocks of recent 1992 Southern California earthquakes is used to produce three-dimensional images of subsurface velocity. The preliminary 1992 dataset, augmented by the 1986 M 5.9 North Palm Springs sequence, consists of 6458 high-quality events recorded by the permanent regional networkÑproviding 76, 306 raypaths for inversion. The target area consisted of a 104 ´ 104 ´ 32 km3 volume divided into 52 ´ 52 ´ 10 rectilinear blocks. Significant velocity perturbations appear to correlate with rupture properties of recent major earthquakes. Preliminary results indicate a low-velocity anomaly separates the dynamic rupture of the M 6.5 Big Bear event from the M 7.4 Landers mainshock; a similar low-velocity region separates the M 6.1 Joshua Tree sequence from the Landers rupture. High-velocity anomalies occur at or near nucleation sites of all 4 recent mainshocks (North Palm Springs- Joshua Tree-Landers-Big Bear). A high-velocity anomaly is present along the San Andreas fault between 5 and 12 km depth through San Gorgonio Pass; this high-velocity area may define an asperity where stress is concentrated and where future large earthquakes may begin.
Back to Lees Publications

An automated method for removing line spectrum elements embedded in colored spectra is presented. Since smooth spectrum estimates are desired, line spectra tend to smear out over an effective smoothing window. This introduces a bias in the spectrum estimation that can seriously degrade determination of signal-to-noise ratios, spectral deconvolution or any other operation where spectrum shape is important in analysis. Multi-taper analysis provides a simple algorithmic approach to this problem and a simple method of determining where spectral peaks are both significant and contain signal power is suggested. While the method is completely general, an illustration of the technique applied to seismic signals is provided. Examples include estimation of signal-to-noise ratio at the high frequency array at, Pinyon Flat, CA. A comparison of noise spectra line segments and signal spectra line spectra reveals similarities associated with instrument noise and shallow resonances that are stimulated by incoming seismic signals. Identification and removal of the resonances provides a better means of estimating background noise spectrum for the purposes of modeling earthquake source spectra and path effects associated with attenuation.
Back to Lees Publications

Tomographic inversion is applied to delay times from local earthquakes to image 3-D velocity variations surrounding the main rupture of the 1989 Loma Prieta earthquake. The 55x45 square km region is represented by blocks of 1 km per side laterally and by 8 layers of varying thickness to 18 km depth. High quality P-wave arrival times recorded on the USGS CALNET array from 549 crustal earthquakes with depths of 0 to 25 km were used as sources. Preliminary results several velocity variations (5-12%) that correlate with specific characteristics of the 1989 rupture. These include prominent high-velocity anomalies near the mainshock hypocenter and prominent low-velocity anomalies where the dip of the San Andreas fault appears to change significantly. The termination of prominent low velocity features existing primarily in the hanging wall to depths of 7-9 km, correlates with the top of the rupture zone. High-velocity variations along the fault dominate where aftershock activity is high. The high velocity anomaly located at depth along the fault is interpreted as imaging the asperity on which the Loma Prieta earthquake occurred.
Back to Lees Publications

We investigate the stability of tomographic analysis by comparing the results of two different methods of parameterizing the three-dimensional P-wave velocity variations in the vicinity of the 1989 Loma Prieta earthquake. A block inversion is implemented using 55x45x10 blocks of 1 km width and varying thickness to a depth of 25 km below the surface. Linear and non-linear analysis are presented. The non-linear analysis is achieved by iterating over three-dimensional raytracing and earthquake relocation relative to current three-dimensional models until solutions show only small improvements. A second parameterization is achieved by using cubic B-spline functions to span the space of the model which is rotated by 46.5û. Non-linear results are presented with several different starting models illustrating the robustness of the technique to the initial conditions. All the non-linear results produced essentially the same final model, which was structurally the same as the model obtained by linear analysis using a reasonable starting model.
Back to Lees Publications

High resolution, three-dimensional images of P-wave velocity anomalies below Mt. St. Helens, Washington, were derived using tomographic inversion. The model is a 27.5´21´20 km target volume parameterized by blocks .5 km per side. The area included 39 stations and 5454 local events leading to 35,475 rays used in the inversion. To diminish the effects of noisy data, the Laplacian was constrained to be small within horizontal layers, providing smoothing of the model. Non-linear effects were compensated for by iterating three-dimensional ray tracing (using pseudo-bending) between inversions and relocating earthquakes relative to the updated three-dimensional model. The structural differences between the linear and non-linear inversions appear to be insignificant, although the amplitudes of the anomalies are larger in the non-linear models. Results indicate a low- velocity anomaly (>7%), approximately 1 km in lateral extent, from 1.5 to 3 km depths. Between 3 and 6 km depth the anomaly appears to spread out. Below 6 km depth the low velocity feature changes to a higher velocity perturbation with lower velocity perturbations flanking around the perimeter of the volcano. The higher velocity material, which correlates with the higher seismicity at that depth, is interpreted as being a plug capping the low velocity magma chamber which begins below 9 km depth.
Back to Lees Publications

Tomographic inversion is applied to delay times from local earthquakes to image three dimensional velocity variations in the Puget Sound region of Western Washington. The 37,500 square km region is represented by nearly cubic blocks of 5 km per side. P-wave arrival time observations from 4,387 crustal earthquakes, with depths of 0 to 40 km, were used as sources producing 36,865 rays covering the target region. A conjugate gradient method (LSQR) is used to invert the large, sparse system of equations. To diminish the effects of noisy data, the Laplacian is constrained to be zero within horizontal layers, providing smoothing of the model. The resolution is estimated by calculating impulse responses at blocks of interest and estimates of standard errors are calculated by the jackknife statistical procedure. Results of the inversion are correlated with some known geologic features and independent geophysical measurements. High P-wave velocities along the eastern flank of the Olympic Peninsula are interpreted to reflect the subsurface extension of Crescent terrane. Low velocities beneath the Puget Sound further to the east are inferred to reflect thick sediment accumulations. The Crescent terrane appears to extend beneath Puget Sound, consistent with its interpretation as a major accretionary unit. In the southern Puget Sound basin, high velocity anomalies at depths of 10-20 km are interpreted as Crescent terrane and are correlated with a region of low seismicity. Near Mt. Rainier, high velocity anomalies may reflect buried plutons.
Back to Lees Publications

Three-dimensional tomographic inversion of P-wave travel-time data is used to investigate the seismic velocity structure of the crust in Thessaloniki and Chalkidiki, N. Greece. Local earthquakes recorded by two networks operating in the area are used as natural seismic sources. Two different target volumes, defined on the surface by 39o50' - 41o50'N and 21o25' - 24o20'E, and 40o10' - 41o 00'N and 22o 45' - 23o 50'E, are investigated. The first dataset is recorded by 13 stations and the second by 29. The size of the blocks used to parameterize the areas is 10 x 10 km and 3 x 3 km in the horizontal, respectively, with varying depth layering. The major seismic velocity anomalies within the crust, obtained by the tomographic inversion, are resolved with a horizontal spatial resolution of about 20 km and 7 km for the first and second target volume, respectively. Our particular interest is to illuminate velocity anomalies and more detailed characteristics of the two main Neogene- Quaternary basins in this region (Vardar-Axios and Struma-Strymon). These basins are identified as low velocity features overlying relatively higher P-wave velocity structures in the lower crust. The complex Mygdonian area reveals a similar pattern of low-velocity basin overlying higher-velocity basement. Overall the velocity patterns correlate well with the location and strike of the main geological and tectonic units in the area, as well as the basic assumptions on basin development. This highlights the utility of local tomography to illuminate structural, tectonic and rheological properties within the crust.
Back to Lees Publications

Tomographic inversion is applied to delay times from aftershocks of the 1986 ML 5.9 North Palm Springs (NPS) earthquake to image 3-D velocity variations within the northern Coachella Valley. P-wave arrival times from 1074 earthquakes, with depths ranging from 3 to 20 km, were used as sources recorded by 12 portable and 4 permanent stations. Preliminary results show well-defined high- and low-velocity anomalies (2-7%) that correlate with the rupture distribution of the 1986 mainshock. At depths less than 8 km, a low-velocity anomaly predominates between the two NE-dipping Banning and Mission Creek faults. From 8 to 12 km, where the NPS mainshock and most of the aftershocks occur, a high-velocity anomaly is observed. This high-velocity feature is interpreted as imaging the asperity responsible for the 1986 rupture; and suggests that velocity information may help to define important elements, such as asperities, that control fault rupture, and thus, may help to predict the location and size of future events.
Back to Lees Publications

Unzen Volcano, located in south west Japan, erupted on November 17, 1990 after 198 years of dormancy and has been ejecting lava since May 20, 1991. In this paper, we present three-dimensional P and S-wave velocity variations below Unzen Volcano using 22,473 P-wave, and 14,349 S-wave arrival times from 3,457 local earthquakes recorded on a local network of 12 seismic stations. The model was parameterized by 24,000 approximately 2.0 km cubic blocks, targeting a volume of 120X100X20 km. A prominent low velocity anomaly (greater than 4\% slowness perturbation) starting from 2.5 km to 5.0 km depth beneath the volcano is observed on both the P and S wave inversions. Below Unzen volcano slightly lower amplitude low velocity anomalies are further observed to a depth of 12.5 km for P and 7.5 km for S-waves. Shallow low velocity anomalies are observed below the Chijiwa Bay area (northwest of Unzen), and deeper anomalies are seen below Shimabara Bay (southwest off the Shimabara Peninsula). This northwest to southeast trending feature represents an elongated system of dykes which supply the shallow magma reservoirs.
Back to Lees Publications

A technique for improving the efficiency of shortest path raytracing (SPR) [Moser, 1991] is presented. We analyze situations where SPR fails and provide quantitative measures to assess the performance of SPR raytracing with varying numbers of nodes. Our improvements include perturbing the ray at interfaces according to Snell's Law, and a method to find correct rays efficiently in regions of low velocity contrast. This approach allows the investigator to use fewer nodes in the calculation, thereby increasing the computational efficiency. In two dimensional cross-borehole experiments we find that with our improvements, we need only use 2/3 as many nodes, saving up to 60\% in time. Savings should be even greater in three dimensions. These improvements make SPR more attractive for tomographic applications in three dimensions.
Back to Lees Publications

Tomographic inversion is applied to 17,659 P phase observations at 21 stations from 2023 earthquakes in the vicinity of Mount St. Helens to study the three-dimensional velocity structure. Block size for the inversion is 2 km horizontally and 2 km or more vertically. Locations of hypocenters are assumed known and are based on a reference one-dimensional, layered velocity structure. A conjugate gradient technique (LSQR) is used to invert the large sparse system of equations, augmented by regularization with a Laplacian roughening matrix. Resolution is estimated by computing the impulse response of the inversion for various critical locations, and uncertainties of the estimates are determined by a jackknife approach. The results of the inversion show a remarkable correlation with known geological and geophysical features. The Spirit Lake and Spud Mt. plutons are characterized by high-velocity regions extending to approximately 9 km depth. The St. Helens seismic zone, a band of diffuse seismicity extending NNW from the volcano is evident as a prominent low-velocity lineation. The change in character of the velocity anomalies south of St. Helens corresponds well with the near cessation of seismic activity there. A low-velocity anomaly beneath the crater from 6 to 16 km depths may represent modern magma accumulations.
Back to Lees Publications

Tomographic inversion is applied to delay times from local earthquakes to image three dimensional velocity variations near Parkfield, California. The 25 ´ 20 square km region is represented by nearly cubic blocks of 0.5 km per side. Arrival times of P waves from 551 local earthquakes, with depths of 0 to 15 km, were used as sources producing 3135 rays covering the target region. The data were recorded on low-noise downhole seismographs. A conjugate gradient method is used to invert the resulting sparse system of simultaneous equations. To diminish the effects of noisy data, the Laplacian of the model parameters is constrained to be small within horizontal layers, providing smoothing of the model. The resolution of the model is estimated by calculating point spread functions at blocks of interest. Estimates of standard errors of the model parameters are calculated by the jackknife statistical procedure. The results of the inversion show correlation with some of the local geological and geophysical features. Station corrections removed the long- wavelength anomaly associated with the contrast of the Salinian block southwest of the San Andreas fault versus the Franciscan to the northeast. A velocity low located a few kilometers northwest of Parkfield (depth 2.5-3.5 km), appears to lie along the gradient of the large Bouguer gravity anomaly associated with the Parkfield syncline. The south- southeastward extension of the low velocities may relate to reflections observed on the Parkfield, Consortium for Continental Reflection Profiling (COCORP) lines. We speculate on the geological meaning of these features and interpret them either as part of the local strike slip tectonics or a shallow crustal detachment. The correlation of higher-velocity features and seismic activity may indicate that earthquakes are occurring in more competent zones while aseismic slip takes place in zones of lower-velocity, less competent rocks.
Back to Lees Publications

Detailed tomographic images of the collision zone between the northern edge of the Philippine Sea and the Eurasian plates reveal a high correlation between tectonic features inferred from seismicity and P and S-wave velocity structures. The linear tomographic inversion covered a 150x150x60 km region in the south Fossa Magna centered on the northern Izu Peninsula. Thirty three 3-component stations were located in the target region and 3823 high quality earthquakes were selected from the catalogues of the NRIESDP, giving rise to 53,593 P and 50,059 S-wave phase arrivals used in the inversion. The model was parameterized by 60x60x10 rectilinear blocks measuring 2.5 km per side horizontally and 5-10 km varying thicknesses in depth. Three-dimensional perturbations from the one dimensional, NRIESDP, layered model were derived by minimizing the squared misfit of the travel time residuals. Regularization was employed to reduce the effects of noisy data by constraining the 2-dimensional Laplacian of the model, within horizontal layers, to be small. The large sparse matrix was solved using the conjugate gradient algorithm LSQR. Reduction of misfit was 50% for the P-wave inversion and 57% for the S-wave inversion.

There is a high correlation between the P and S-wave inversions. Between 5-15 km depth both 3D images indicate low velocity anomalies trending northward in the west along the Fujikawa, from Suruga Bay to northwest of Mt. Fuji. North of Izu Peninsula there is a broad high velocity zone. These results correlate with Bouguer gravity anomalies previously observed. Starting at 5 km depth, a high velocity, northeast dipping anomaly beneath Sagami Bay is inferred to be the shallow subducting PHS plate. To the north, a high velocity zone in the vicinity of the East Yamanashi seismic swarm (15-32 km depth) is clearly observed on the P-wave inversion and to a lesser extent on the S-wave inversion. Because of the high seismicity in this region, we speculate this is the zone where the PHS plate is grazing the EUR plate to the northwest, and the high velocity anomaly represents a large scale asperity. To the west, below the Mt. Fuji area, a broad low velocity is observed between 15-32 km depth. We interpret the deep, low velocity structure, which exhibits little seismicity, to be either the subduction of the volcanic arc or a region of elevated temperature associated with hot rising magma. Near the swarm area on the eastern edge of Izu Peninsula, low velocity anomalies are observed near the surface on both P and S wave inversions. Along the eastern part of Izu Peninsula a prominent low velocity anomaly is observed between 10-20 km on the S-wave image but is less pronounced on the P-wave inversion.
Back to Lees Publications

We compare two approaches for solving the large, sparse linear systems that arise in tomographic velocity inversion problems. When noise is present in the data, the system typically is inconsistent and quasi-overdetermined, and some form of regularization must be implemented to avoid the strong, undesired influence of small singular values dominating the solutions. First, a Bayesian ART (Algebraic Reconstruction Technique) algorithm is applied to the system where we solve for a set of model parameters and residuals simultaneously. Careful choice of relaxation parameters and smoothing filters insures convergence and acceptable results. This approach avoids the undesirable effects of implicit row weighting inherent in simple ART or SIRT (Simultaneous Iterative Reconstruction Technique) applications. Second, a conjugate gradient approach is implemented via algorithm LSQR where regularization is achieved by augmenting the system with additional constraint equations which minimize the roughness of the model. Specifically, applied to a 3-dimensional tomographic inversion we constrain the second derivative (the Laplacian) to be zero within horizontal layers. Comparison on synthetic data reveals that these techniques produce nearly equivalent results. By applying these methods to local earthquake data in the vicinity of Mount St. Helens, Washington, we have produced a 3-dimensional, laterally varying velocity structure in the top 40 km of the crust which correlates well with known geological and geophysical features and delineates possible accumulation of magma beneath the crater.
Back to Lees Publications

Microseismicity in the Coso geothermal field has been searched for seismic doublets, hypocenters co-located that appear to have identical waveforms. Using 1085 high quality events from 1993-1994, numerous doublets were identified, some occurring within minutes of each other. Hypocentral data were subdivided into spatial clusters to reduce the computational burden and multiple cross correlation pairs were evaluated, assigning scores to each pair. In one spatial cluster including 183 events (16471 pairs) yielded 96 high correlation (>0.6) paired events. Equivalence class analysis and cluster analysis routines were used to isolate potential multiplets. Among the 96 high correlation pairs 24 equivalence classes were isolated. While most of these are doublets, 8 classes included 3 or more cluster members and one class included 16 members. Relative locations are calculated using phase shifts between corresponding events. Detailed analysis of hypocenter relocations shows elongate, vertical structure aligned along ambient stress fields as measured by conventional focal mechanism analysis. Multiplet hypocenters are generally oriented sub-parallel to faults observed at the surface.
Back to Lees Publications

R is a free software for statical computing and graphics. It compiles and runs not only on UNIX platforms but MS Windows. The R commands are easy and offer interactive help. R is used in extensive field by implementing packages. RSEIS, the package of R, enables us to do easy graphical processing of seismic data. Here we illustrate by showing an example of seismic data processing using RSEIS.
Back to Lees Publications

Although widely used for preparation of geographic maps in the field of earth sciences, the Generic Mapping Tools(GMT) software is difficult for users to understand, and does not work well with Microsoft (MS) Window PC. By utilizing R package, GEOmap, we can do mapping work at MS window PC with commands easier than those of GMT. In addition, the R commands provide interactive help. Here we briefly introduce a few features of GEOmap, and illustrate a procedure for preparing topographic maps as an example.
Back to Lees Publications

Address:

  
Prof. Jonathan M. Lees
Department of Geological Sciences
CB #3315, 104 South Road
University of North Carolina
Chapel Hill, NC  27599-3315
(919) 962-0678
FAX (919) 966-4519