University of Southern California Dissertations and Theses

Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults

(USC Thesis Other)

Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content ANALYSIS OF WAVEFORM AND CATALOG DATA OF AFTERSHOCKS FOR
PROPERTIES OF EARTHQUAKES AND FAULTS

by

Wenzheng Yang

A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSTIY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(EARTH SCIENCES)

December 2009

Copyright 2009 Wenzheng Yang

ii
Acknowledgments

Even though I am in a place far from my parents, my feeling for them is always
close. Thanks to the very cheap IP phone, I call them each weekend and tell them what
has been happening to me. They always comfort and encourage me to continue chasing
for my PhD dream, which has existed in my mind ever since I was a pupil. Without
those picture books about Thomas Edison, Albert Einstein, Marie Curie, Ernest
Rutherford, which were purchased by my parents during my childhood, I might not
have such strong determination and curiosity to do research.
Prof. Ma Li passed away two days before the Beijing Olympic games in 2008 due to
a kind of painful blood cancer. In 1996, she led me to the gate of earthquake physics
and statistics research and encouraged me to enter. She worked hard with all kinds of
passion and enthusiasm, and she treated young researchers like her own children. It was
Prof. Ma Li who introduced me to Prof. Yehuda Ben-Zion in 2001, who later became
my Ph.D. advisor.
I am deeply indebted to my advisor Prof. Yehuda Ben-Zion. He guided me in
research and trained me to be professional. I very much appreciate his generosity. I
could not remember how many times he invited us for dinner. I also deeply appreciate
Mrs. Ben-Zion, who was always kind and prepared lots of delicious and healthy food
for us at the Thanksgiving nights in the last several years. Additionally, I appreciate
Yehuda for his organization of several hiking, camping and fieldwork during these five
years.

iii
Prof. Richard Aster was my Master thesis advisor. I learned lots of computational
skills and seismic waveform analysis techniques from him. Prof. Ta-liang (Leon) Teng
is a very knowledgeable person. One day he told me a four Chinese characters
“ 业 ”, which comes from an ancient proverb by Han Yu ( 韩 , 768～824), and
its English meaning is “The expectations of life depend on diligence”. Nowadays, there
are so many things to distract me from concentrating on my work. I choose this proverb
as my motto, and try to be focused, work hard and also be consistent in both research
and life. Prof. Zhigang Peng taught me how to work with the Turkey earthquake data,
and communicated with me a lot during these five years. He always gave me wonderful
suggestions and advice on both work and life.
During my five years studying at USC, I have accumulated a list of people to
appreciate. If someone is forgotten, I hope it is forgivable. For people who work or
teach at USC: Cindy Waite, John Yu, John McRaney, Prof. Trifunac Mihailo, Prof.
Charlie Sammis, Prof. Tomas Jordan, Prof. Dahmen Karin, Prof. John Nodvik, Prof.
Thorsten Becker, Prof. Aiichiro Nakano, Prof. James Dolan, Prof. Yonggang Li, Prof.
Steven Scott, Prof. Tomas Henyey, Prof. David Okaya, Prof. Maxim Olshanii, Vardui
Ter-Simonian, Barbara Grubb, Stanley Wright, Sally Henyey and Tran Hyunh. My
fellow students offered me tremendous help in both my academic and social life: Po
Chen, Iain Bailey, Michael Lewis, Zheqiang Shi, Jeremy Zechar, Adam Fischer, Peter
Power, Ory Dor, Neta Wechsler, Tao Zhang, Yanbing Wang, Jinjun Kan, Shiqing Xu,
Mengfan Zhu, Feng Wang, Zi-Yu Wu and Thomas Goebel.
My special appreciation goes to Miss Su Li.

iv
Table of Contents
Acknowledgments ii
List of Tables vi
List of Figures vii
Abstract ix
Introduction 1
Chapter 1: Observational analysis of correlations between aftershock
productivities and regional conditions in the context of a damage
rheology model
4
Summary 4
1.1 Introduction 5
1.2 Analysis 9
1.2.1 Individual aftershock sequences of large mainshocks 11
1.2.2 Stacked aftershock sequences of moderate mainshocks 16
1.3 Discussion 23
Chapter 1 References 28
Chapter 2: Variations of strain-drops of aftershocks of the 1999 İzmit and
Düzce earthquakes around the Karadere-Düzce branch of the North
Anatolian Fault
31
Summary 31
2.1 Introduction 32
2.2 Data and Preprocessing 35
2.2.1 Data 35
2.2.2 Preprocessing 36
2.3 Method 38
2.3.1 Iterative Separation of Source, Station and Travel-Time
Spectra
39
2.3.2 Stacking of Source Terms with Similar Amplitudes 41
2.3.3 Fitting Source Spectra with Theoretical Model 43
2.4 Results 45
2.4.1 Spatial Variations of Strain-drops 46
2.4.2 Temporal Variations of Strain-drops 51
2.4.3 Station Spectral Terms 54
2.5 Discussion 57
Chapter 2 References 62

v
Chapter 3: Estimation of corner frequency ratios from P and S waves
generated by aftershocks around the Karadere-Düzce branch of the
North Anatolian Fault
65
Summary 65
3.1 Introduction 66
3.2 Data and Preprocessing 69
3.3 Joint Fitting P and S wave for Corner Frequency Ratio 73
3.4 Application and Results 76
3.5 Discussion 84
Chapter 3 References 87
Chapter 4: An algorithm for detecting clipping waveforms and suggested
correction procedures
90
Summary 90
4.1 Introduction 91
4.2 Data 93
4.3 Features of Clipped Waveforms 95
4.4 An Algorithm for Detecting Clipping Waveform 96
4.5 The Kriging and the Similar Waveform Methods 102
4.6 A Comparison of Two Correction Methods 104
4.7 Discussion and Conclusions 111
Chapter 4 References 114
Alphabetized Bibliography 116

vi
List of Tables
Table 1.1 Estimated χ and R values for the different regions

26
Table 2.1 Statistical features of the four divided sections

49
Table 3.1 Five schemes for P and S wave spectra calculations

73
Table 4.1 A summary of the observed ranges, total numbers of waveforms
and numbers of detected clipped waveforms for each of the 3
components of the ten stations

99
Table 4.2 Statistical summary of correction errors for methods comparison

110

vii
List of Figures
Figure 1.1 M ≥ 2.0 earthquakes and heat flow data in Southern California

10
Figure 1.2 A comparison of two individual aftershock sequences

12
Figure 1.3 Aftershock productivities v.s heat flow for two aftershock areas

15
Figure 1.4 A comparison of five studied regions

18
Figure 1.5 Aftershock productivities and L factor

21
Figure 1.6 Aftershock productivities, heat flow and thickness of sediment
cover

23
Figure 2.1 Aftershocks along the Karadere-Düzce region

36
Figure 2.2 Example of waveforms and spectra

37
Figure 2.3a Source, station and travel-time terms separation

41
Figure 2.3b Separated source term and stacking

42
Figure 2.4 Fiting for strain-drop

45
Figure 2.5a Spatial distribution of best-fitting strain-drop

46
Figure 2.5b Examples of source spectra fitting

48
Figure 2.6 Depth variations of log
10
strain-drops along the strike

51
Figure 2.7 Strain-drop variations before and after the Düzce mainshock

52
Figure 2.8 Temporal variation of strain-drops in four sections

54
Figure 2.9 Station spectra

55
Figure 3.1 The distribution of 8785 Aftershocks along the Karadere-Düzce
region

69
Figure 3.2 The selection wave windows for scheme A, B and C

72
Figure 3.3 Statistical features of best-fitting results for the five schemes

77

viii
Figure 3.4a Best-fitting P and S spectra of twelve examples for scheme A

78
Figure 3.4b Best-fitting P and S spectra of twelve examples for scheme B

79
Figure 3.4c Best-fitting P and S spectra of twelve examples for scheme C1

80
Figure 3.4d Best-fitting P and S spectra of twelve examples for scheme C2

81
Figure 3.4e Best-fitting P and S spectra of twelve examples for scheme C3

82
Figure 3.5 Spatial distribution of best-fitting corner frequency ratios

83
Figure 4.1 Epicentral distribution of ~26000 earthquakes

94
Figure 4.2 One clipped event recorded by station VO

95
Figure 4.3 Distribution of maximum amplitudes of 24135 waveforms

97
Figure 4.4 The variation of detected clipped waveforms to amplitude
thresholds

98
Figure 4.5 Distribution of detected clipped events

100
Figure 4.6 Magnitude and distance relationship of detected clipped
waveforms

101
Figure 4.7 Probability density of cases of consecutive clipped waveforms

102
Figure 4.8 Example of correction using similar waveform method

104
Figure 4.9 One selected waveform with its similar waveforms

106
Figure 4.10 Correction comparisons to cases of 1-4 consecutive clipped points

107
Figure 4.11 A comparison of correction errors of two methods

109
Figure 4.12 A performance comparison between two correction methods

111

ix

Abstract

We study the statistical relationship between aftershock productivity and background
geophysical features using an earthquake catalog with ~340,000 events that occurred
between 1984 and 2004 in the Southern California region. We analyze properties of
individual aftershock sequences and the average properties of stacked aftershock
sequences in five regions. The results indicate that productivities of the aftershock
sequences are inversely correlated with heat flow and existence of deep sedimentary
covers, in agreement with the damage model predictions given by Ben-Zion and
Lyakhovsky (2006). Using the observed ratios of aftershock productivities, along with
simple expressions based on the damage model, we estimate the seismic coupling
coefficients in the different regions. We estimate seismic strain-drops from a data set
with ~26,000 aftershocks that occurred along the Karadere-Düzce branch of the North
Anatolian Fault after the 1999 İzmit mainshock. We use a method that is associated
with separation of source, travel-time and station spectral terms, and stacking results at
several stages to enhance the signal-to-noise ratio. The strain-drops are obtained by
fitting iteratively the separated source spectral of 201 nearest neighboring events in
different amplitude bins to the omega square source spectra model. We analyze the
variations of the obtained strain-drops along the fault zone, with the depth extension.
Furthermore, we analyze seismic source properties using both P and S waves on
seismograms. We develop a procedure to fit for P/S corner frequency ratios together
with relative strain-drop systematically. We consider five different schemes to analyze

x
the S wave, and the obtained median values of frequency ratios from all schemes are
larger than 1.5. We describe an algorithm to detect clipped waveforms in the same
~26000 data set based on statistical characteristics of observed amplitude ranges,
properties of neighboring points in seismograms. To the detected clipped waveforms,
the intensity of waveform clipping increases with proximity to the fault and site effects.
We compare two correcting methods (the Kriging method and similar waveform
method) to several simple cases with 1-9 consecutive clipped points on artificially
clipped waveforms.

1
Introduction
Aftershocks are generally triggered by a mainshock and occur very close to its
hypocenter. Compared with other types of earthquakes, there are several advantages to
the study of aftershocks. Firstly, aftershocks can be relatively well observed, since most
of them occur along the fault zone associated with the mainshock and instruments could
be set up accordingly shortly after the occurrence of the mainshock. Secondly,
aftershocks provide a large data set of earthquakes concentrated along the fault zone.
Therefore, abundant small-scale information related to fault zone and large-scale
background geophysical characteristics could be disclosed through the analysis of
aftershock data. In my dissertation, we study large data sets of aftershocks using two
approaches in regions with different scales.
In Chapter 1, we study the statistical relationship between aftershock productivity and
background geophysical features by using a relocated earthquake catalog with ~340,000
events that occurred between 1984 and 2004 in the Southern California region (500 ×
500 km-square). Ben-Zion and Lyakhovsky (2006) developed a solution for aftershocks
decay in a damage rheology model. This solution indicates that the productivity of
aftershocks decreases with increasing temperature, and also with existence of
sedimentary rocks at seismogenic depth. We first analyze properties of individual
aftershock sequences generated by the 1992 Landers and 1987 Superstition Hills
earthquakes. The results show that the ratio of aftershock productivities in these
sequences spanning 4 orders of event magnitudes is similar to the ratio of the average
heat flow in the two regions. To perform stronger statistical tests we systematically
2
analyze the average properties of stacked aftershock sequences in five regions. The
results indicate that the productivities of the stacked sequences are inversely correlated
with the heat flow and existence of deep sedimentary covers, in agreement with the
damage model predictions. Using the observed ratios of aftershock productivities, along
with simple expressions based on the damage model, we estimate the seismic coupling
coefficient in the different regions. The employed methodology for estimating the
seismic coupling in different regions can be useful for seismic hazard studies.
In Chapter 2, we estimate seismic strain-drop distributions from a data set with
~26,000 aftershocks that occurred along the Karadere-Düzce branch (30 × 100 km
square) of the North Anatolian fault after the 1999 Izmit mainshock. The strain-drop is
a key parameter needed to understand the physics of earthquakes and infer the stress
level in the crust. We use a method that is associated with separation of source, travel-
time and station spectral terms, and stacking results at several stages to enhance the
signal-to-noise ratio. The strain-drops are obtained by fitting iteratively the separated
source spectral of 201 nearest neighboring events in different amplitude bins to the
omega square source spectra model. We analyze the variations of the obtained strain-
drops along the fault zone, with the depth extension, and with time respectively.
In Chapter 3, with the analysis infrastructure developed in Chapter 2, we extend the
analysis of seismic source properties to include both P and S waves, and develop a
method to estimate the P/S corner frequency ratios together with relative strain-drops.
We take several technique improvements to decrease individual spectra variations and
avoid inconsistencies related with practical use of standard seismic source model. We
take a 1.28 s window around P arrival for P wave. For S wave, we consider several
3
different schemes: A) S wave in a single 1.28 s window around S arrival. B) S wave in a
single 10 s window after S arrival. C) S wave in consecutive multiple (N=15) 1.28 s
windows. For scheme C, we consider three different normalizing approaches to stack
the spectra in multiple windows: C1) Normalize by travel-time difference. C2)
Normalize by amplitude ratio. C3) Without normalization. The mode values for the
obtained the 8,785 corner frequency ratios are 1.52, 1.57, 1.89, 1.47 and 1.56 for each
of the five schemes, respectively. The obtained mode values for log10 relative strain-
drops for all schemes are -4.35. This corresponds to a stress drop of 1.34 MPa assuming
a nominal rigidity of 30 GPa.
In Chapter 4, we describe a clipped waveform detection algorithm that is based on
observed ranges of amplitudes in the recorded seismograms, expected frequency-size
statistics of amplitudes associated with the Gutenberg-Richter distribution and
properties of neighboring points in seismograms. We apply this algorithm to detect
clipped waveforms automatically inside the ~ 26, 000 aftershocks data set from Chapter
2. To the detected clipped waveforms, the horizontal components have more detected
clipped waveforms than the vertical one, and the intensity of waveform clipping
increase with proximity to the fault and amplitude of site effects. We compare two
potential methods (the Kriging method and the similar waveform method) for correcting
clipped waveforms to simple cases with 1-9 consecutive clipped samples on artificially
clipped waveforms.
Through systematically and comprehensively analyzing aftershock data with different
scales and with regional comparison, using both seismicity and seismic waveforms, we
can better understand the properties of earthquakes and faults.
4
Chapter 1

Observational analysis of correlations between
aftershock productivities and regional conditions in the
context of a damage rheology model

SUMMARY
Aftershock sequences are commonly observed but their properties vary from region
to region. Ben-Zion and Lyakhovsky (2006) developed a solution for aftershocks decay
in a damage rheology model. The solution indicates that the productivity of aftershocks
decreases with increasing value of a non-dimensional material parameter R given by the
ratio of timescale for brittle deformation to timescale for viscous relaxation. The
parameter R is inversely proportional to the degree of seismic coupling and is expected
to increase primarily with increasing temperature, and also with existence of
sedimentary rocks at seismogenic depth. To test these predictions we use aftershock
sequences from several southern California regions. We first analyze properties of
individual aftershock sequences generated by the 1992 Landers and 1987 Superstition
Hills earthquakes. The results show that the ratio of aftershock productivities in these
sequences spanning 4 orders of event magnitudes is similar to the ratio of the average
heat flow in the two regions. To perform stronger statistical tests we systematically
analyze the average properties of stacked aftershock sequences in five regions. In each
5
region we consider events with magnitudes between 4.0 and 6.0 to be mainshocks. For
each mainshock, we consider events to be aftershocks if they occur in the subsequent 50
days, within a circular region that scales with the magnitude of the mainshock, and in
the magnitude range between that of the mainshock and 2 units lower. This procedure
produces 28-196 aftershock sequences in each of the five regions. We stack the
aftershock sequences in each region and analyze the properties of the stacked data. The
results indicate that the productivities of the stacked sequences are inversely correlated
with the heat flow and existence of deep sedimentary covers, in agreement with the
damage model predictions. Using the observed ratios of aftershock productivities, along
with simple expressions based on the damage model, we estimate the relative values of
the material parameter R and seismic coupling coefficient in the different regions. The
employed methodology for estimating the seismic coupling in different regions can be
useful for seismic hazard studies.

1.1 Introduction
Aftershocks are a form of brittle relaxation of the crust to rapid stress loading
generated by previous earthquakes that are referred to as mainshocks. It may thus be
expected that aftershock properties would depend on certain regional conditions such as
rock type, ambient temperature and fluid content. Such correlations have been found in
previous observational studies, but mechanical models that make these connections
explicit, and related observational analyses that employ the connections to derive
additional aspects of earthquake behavior, have been generally lacking (e.g. Ben-Zion
2008).
6
The decay rates of aftershock sequences are usually fitted by the Omori-Utsu law

ΔN /Δt = K(t+ c)
− p
, (1.1)
where N is the cumulative number of events, t is the time after the mainshock, K is the
productivity of the sequence, the exponent p has a value close to 1, and c is usually a
small fraction of a day (e.g. Omori 1894; Utsu et al. 1995). However, aftershock decay
rates can also be fitted with exponential and other functions (e.g. Mogi 1967; Otsuka
1985; Gross and Kisslinger 1994; Narteau et al. 2002). Observational studies indicate
that cold continental regions have high aftershock productivity and long sequences
associated with low effective power law decay, while hot continental regions and
oceanic lithosphere have low aftershock productivity and short sequences with fast
decay or exponential fall-off (Mogi 1967; Kisslinger and Jones 1991; Davis and
Frohlich 1991; Utsu et al. 1995; McGuire et al. 2005). However, as far as we know, a
systematic analysis of the connections between aftershock properties and heat flow
within a formalism associated with a physical model and specific quantitative
expectations has not yet been done. In the present paper we perform such an analysis in
the context of a damage rheology model (Ben-Zion and Lyakhovsky 2006) that makes
explicit connections between the ambient temperature field and other conditions that
may modify the effective viscosity of a region, the productivity of aftershock sequences,
and the degree of seismic coupling in a region.
Ben-Zion and Lyakhovsky (2006) derived an analytical solution for aftershock decay
rates in a viscoelastic continuum damage model for evolving elastic properties of rocks
sustaining irreversible brittle deformation. The damage model generalizes the strain
energy function of a solid to account for first-order macroscopic effects of existing
7
cracks (damage), and makes the elastic moduli functions of an evolving damage state
variable α representing the local crack density (e.g. Lyakhovsky et al. 1997; Hamiel et
al. 2004). An undamaged material with α = 0 is the ideal solid governed by 3D linear
elasticity, while a material with α = 1 can not support any load. Lyakhovsky and Ben-
Zion (2008) and Ben-Zion (2008) provide detailed reviews of the damage model. The
analytical solution of Ben-Zion and Lyakhovsky (2006) follows an exponential form
that depends fundamentally on a non-dimensional material parameter R given by the
ratio of timescale for brittle damage to timescale for viscous relaxation. The material
parameter R also plays a fundamental role in the degree of seismic coupling χ, which
quantifies the fraction of the stored elastic strain energy that is released in brittle
deformation (i.e., earthquakes). Ben-Zion and Lyakhovsky (2006) showed that χ is
inversely proportional to R and is given by
. (1.2)
Since the degree of seismic coupling has direct consequences for the seismic potential
of a region, an ability to estimate the value of χ from properties of aftershock sequences
can have important practical applications.
Ben-Zion and Lyakhovsky (2006) showed that for simple assumptions that are
expected to hold over relatively short time intervals (e.g. 100 days), the analytical
solution for aftershocks decay rates follows the power law relation
, (1.3a)
8
where is the initial aftershocks decay rate, φ is a small scaling constant (e.g. 10
−6
)
between the number of events and average change of rock damage, and is the initial
rock damage. Equation (1.3a) provides the following mapping between parameters of
the Omori-Utsu relation (1.1) and parameters of the damage model,
, (1.3b)
,
p =1.
In particular, the ratios of the productivity values K of aftershock sequences in different
regions, derived by fitting the aftershock decay rates to the Omori-Utsu law with an
assumed p = 1, are expected to be inversely proportional to the R values that
characterize the different regions. Since R is inversely proportional to the effective
viscosity in a region, the productivity of aftershock sequences should decrease for
lithological and ambient conditions that reduce the viscosity. These include increasing
temperature, increasing fluid content, and existence of sediments (or other rock units
that have relatively low viscosity) at seismogenic depths.
In the next section we test the forgoing expectations using observed aftershock
sequences in several sub-regions of southern California. We adopt a simple approach
associated with (a) selecting a-priori several sub-regions using available information on
the heat flow and crystalline vs. sedimentary rocks at seismogenic depth, and (b)
performing a comparative analysis of aftershock productivities in those regions using
identical values for all other parameters that characterize the entire southern California
region. Other possible controlling variables such as fluid content are not examined
9
because of the lack of systematic information that is needed for our comparative study.
The obtained productivities of aftershock sequences are found to be inversely correlated
with the heat flow and existence of deep sedimentary rocks, in agreement with the
damage model predictions. Using equations (1.2) and (1.3) and the derived aftershock
productivities, we estimate the relative degrees of seismic coupling in the different
regions. Since the regional conditions assumed by the damage model to affect the
seismic coupling and aftershock productivities evolve only over geological time, the
employed approach (and related more-sophisticated analysis procedures) may be useful
for estimates of seismic hazard.

1.2 Analysis
We analyze aftershock properties in different regions of southern California using
two separate approaches. In the first approach we compare individual aftershock
sequences of large mainshocks that occurred in regions with known strong contrast of
heat flow. To perform a more systematic analysis based on larger amount of data, we
compare in the second approach properties of stacked aftershock sequences of moderate
mainshocks in five regions that have clear differences in heat flow and sedimentary
thickness. We distinguish primarily between areas with thin sediments (e.g. the Mojave
Desert) where essentially all the seismicity occurs in dense cohesive crystalline rock,
and areas with deep sedimentary covers (e.g. Ventura Basin) where much of the
seismicity occurs within more compliant sediments. Intermediate levels of sedimentary
covers are referred to as regular or moderate and are assumed to have small effects on
the results. More detailed information on the rocks types can be obtained from the
10
SCEC velocity model (http://www.data.scec.org/3Dvelocity/), but this is not needed for
the present analysis level. Figure 1 shows a map view of the employed seismicity and
heat flow data sets. The seismic data are taken from the relocated 1984-2002 southern
California earthquake catalog of Shearer et al. (2005). The heat flow data are taken
from the USGS online heat flow database
(http://earthquake.usgs.gov/heatflow/index.html).

Figure 1.1. A map view showing epicenter locations of M ≥ 2.0 earthquakes (black dots) and heat flow
data (colored squares with values indicated in the legend) used in this work. The earthquakes are from the
relocated catalog of Shearer et al (2005) and the heat flow values are from the USGS online heat flow
database.

11
1.2.1 Individual aftershock sequences of large mainshocks
To compare properties of aftershocks of large individual mainshocks, we choose the
aftershock sequences of the 1992 M7.3 Landers earthquake and 1987 M6.6 Superstition
Hills event. As shown in Figure 1.1, the area of the Landers sequence (region B) has a
relatively low heat flow, while the area of the Superstition Hills sequence (region A)
has a relatively high heat flow. We assume that the aftershock zone of each event
consists of a 1°×1° box centered on the mainshock epicenter (Figures 1.2a and 1.2d).
We also consider a situation where the Landers aftershocks are limited to the polygon
region in Figure 1.2a. We consider events inside the employed areas to be aftershocks if
they occurred in the subsequent 50 days after the occurrence of the mainshock. This
somewhat arbitrary choice of time interval includes the most significant portions of the
aftershock sequences. We performed similar analyses with time windows ranging from
30 to 60 days and our main results remain essentially unchanged.

12

Figure 1.2. (a) Epicenters of analyzed aftershocks (circles) of the 1992 M7.3 Landers California
earthquake (red star). Events inside and outside the dashed polygon around the main rupture zone are
marked with blue and green colors, respectively. (b) Magnitude versus time of the Landers aftershocks in
(a) occurring in the first 50-day after the mainshock. (c) Numbers of Landers aftershocks per day (circles)
and fitted curves based on an integral form of the Omori-Utsu relation (equation 5b) with t
0
= 5 day. The
black and blue symbols correspond to events in the whole region and the polygon area, respectively. For
presentation purpose, we use instead of as the scale of the horizontal axis. The fitted K
values of the aftershock productivities are indicated below the fitted lines. (d) Same as (a) for the
aftershock sequence of the 1987 M6.6 Superstition Hills earthquake. (e) Magnitude versus time of the
analyzed aftershocks of the Superstition Hills earthquake. (f) Same as (c) for aftershocks of the
Superstition Hills earthquake.

In order to isolate the effect of rheological properties on the productivity of
aftershock sequences we have to account for the different magnitudes of the
mainshocks. Helmstetter et al. (2005) showed that the average productivity K of
aftershocks in southern California depends on the mainshock magnitude M as
, (1.4)
where λ = 1.05 ± 0.05 and M
min
is the magnitude cutoff. Felzer et al. (2004) obtained an
average λ value for California seismicity of about 1. To account for the size difference
13
of mainshocks on the aftershock productivities, we assume that λ = 1 and use in each
sequence aftershocks that are within the same magnitude range below that of the
mainshock. Specifically, we analyze properties of aftershock sequences that are within
4 magnitude units below those of the mainshocks. This choice allows us to include the
largest number of events that are above the completeness level of the catalog for the two
sequences. The latter is tested simply by fitting the frequency-magnitude statistics of the
data with the Gutenberg-Richter relation (e.g. Wiemer and Wyss 2000).
The temporal evolution of the employed aftershock sequences of the Landers and
Superstition Hills mainshocks are shown in Figures 1.2b and 1.2e, respectively. To
estimate aftershock productivities that are associated with the Omori-Utsu law, we use
an integral form of the Omori-Utsu law with p = 1 given by
, (1.5a)
where t
0
is the initial time. As noted earlier the constant c is typically a small fraction of
a day. The determination of c is strongly affected by the early part of the aftershock
sequence, which is associated with incomplete recording (e.g. Utsu et al. 1995; Peng et
al. 2006; Enescu et al. 2007). To reduce artifacts associated with c, we choose a
relatively large t
0
value so that c can be effectively omitted. We thus use in the analysis
, (1.5b)
where We discretize the time into successive 1 day intervals, sum the
number of aftershocks that occur inside each interval, and use a nonlinear least-squares
method (Bates and Watts 1988) to estimate the productivity parameter K.
14
Figures 1.2c and 1.2f show the fitted results (solid lines) and obtained K values for
the Landers and Superstition Hills aftershock sequences using t
0
= 5 day. The estimated
productivity values remain essentially the same with t
0
values of 4, 3, 2 and 1 day. The
deviations of the data from the fitted straight lines may be related to the fact that we fix
the exponent p to be 1, ignore possible effects of the c parameter by using large t
0
value,
and other simplifications of our analysis. However, the overall trends of the cumulative
number of aftershocks follow the fitted curves, and the employed procedure is not
expected to produce a bias that would affect our comparative study of aftershock
productivities in the different regions. Since Figure 1.2a includes aftershocks of the
1992 M6.4 Big Bear event, which may or may not be considered as belonging to the
Landers aftershock sequence, we also analyze with identical procedure only events that
occurred in the polygon region of Figure 1.2a around the main Landers rupture zone.
The resulting numbers of aftershocks (plus symbols) and associated productivity are
indicated in Figure 1.2c. The productivity value of the aftershocks in the polygon region
is smaller than that of the entire box in Figure 1.2a, but it is still larger than the
productivity value of the Superstition Hills aftershock sequence.
15

Figure 1.3. (a) Numbers of aftershocks per day (circles) and fitted curves (lines) based on the discrete
Omori-Utsu relation (equation 1) with c = 0, p = 1 and t
0
= 5 day for the Landers aftershock sequence
(blue) and the Superstition Hills aftershock sequence (red). The blue circles and solid lines are for all the
events in Figure 2a, while the plus symbols and dashed line are for the events that occur only inside the
polygon area. (b) Relations between productivity and heat flow values for aftershocks of the Landers and
Superstition Hills mainshocks. The horizontal error bars represent 95% confidence intervals of the heat
flow data. The vertical error bars represent the standard deviations of the productivities.

Figure 1.3a displays the observed aftershock sequences of the Landers and
Superstition Hills mainshocks, along with lines corresponding to the discrete form of
the Omori-Utsu law (equation 1.1) with the parameters associated with the fits of
Figures 1.2c and 1.2f (c = 0, p = 1, t
0
= 5). The aftershock productivity in the Landers
area is seen to be about 2.4 (or 1.5 using only the polygon region) times higher than it is
in the Superstition Hills area. To compare the obtained aftershock productivities with
the background heat flow, we use the USGS heat flow data (Figure 1.1) inside the
aftershock zones. We take the mean of the data in each region as a representative value
16
of the heat flow, and the 95% confident interval of the data to be the error range. The
relationship between heat flow values and aftershock productivities for the different
aftershock zones is displayed in Figure 1.3b. The results show inverse correlation
between the ratios of heat flow and aftershock productivity values in the Landers and
Superstition Hills aftershock zones. It is clear, however, that analysis associated with a
small number of sequences (here only 2) may be affected strongly by statistical
fluctuations rather than giving an appropriate representation of the relations between
aftershocks and environmental properties. We therefore expand our analysis in the next
section to include additional regions, with tens of aftershock sequences in each region
and consistent definition of the aftershocks area of each mainshock.

1.2.2 Stacked aftershock sequences of moderate mainshocks
To reach stable estimations of aftershock productivities that are based on a larger
amount of data, we analyze properties of stacked aftershock sequences in five
representative southern California regions (Figure 1.1). The regions are selected, as
before, based on well known environmental conditions that are assumed by the damage
model of Ben-Zion and Lyakhovsky (2006) to affect the aftershock properties. The
stacking of data smoothes out fluctuations associated with individual aftershock
sequences and highlights the common features of all sequences inside each region.
The analysis is done as follows. We select aftershock sequences of mainshock
magnitudes between 4.0 and 6.0 in each of the five regions of Figure 1.1. The areas of
the regions are listed in Table 1.1 and decrease in the following order: A, B, D, C and E.
The Imperial Valley and Coso regions (A and E) have relatively high heat flow levels,
17
the largest and smallest areas, and regular to moderate sedimentary covers. The Landers
and Hector-Mine area (B) has low heat flow, relatively large area, and very thin
sedimentary cover (order 100-200 m) that is above the seismicity. The San Bernardino
region (C), which includes both the San Bernardino valley and mountains, has low
average heat flow, moderate area, and moderate average sedimentary cover. Finally, the
Ventura Basin region (D) has low heat flow, moderate area, and thick sedimentary
cover (order 10 km) that extends nearly to the bottom of the seismogenic zone.
18

Figure 1.4. Results associated with analysis of stacked aftershock sequences in the five regions of Figure
1.1. Each row gives results for the region with the indicated name. The left panels show magnitude versus
time of aftershocks with 2.0 ≤ M ≤ 6.0. The middle panels display the number of events per 0.5 day for all
aftershock sequences in a given region (colored curves with plus symbols) and the corresponding stacked
aftershock sequences (bold black points). The number of employed sequences in each region is shown at
the top right. The right panels show the stacked number of aftershocks per 0.5 day (circles) and fitted
curves (lines) based on the integral Omori-Utsu relation (equation 1.5b) with t
0
= 5 day. The fitted K
value of the aftershocks productivity in each region is indicated at the bottom right.
19

The left panels in Figure 1.4 show magnitude versus time of earthquakes with 2.0 ≤
M ≤ 6.0 that occurred in each region between 1984 and 2002. We consider events to be
aftershocks of moderate mainshocks with 4.0 ≤ M ≤ 6.0, if they occur in the subsequent
50 days within a circular region that scales with the magnitude of the mainshock, and in
the magnitude range between that of the mainshock and 2 units below. The use of
circular regions to represent aftershock areas is appropriate for moderate mainshocks
with M ≤ 6.0, which generally do not saturate the seismogenic zone (e.g. Ben-Zion
2008), and it allows us to perform a systematic analysis with aftershock areas that scale
with the sizes of the mainshocks.
We estimate the rupture radius (r) of each mainshock using the following scaling
relation (e.g. Ben-Zion 2003) for a circular crack with a uniform stress drop,
, (1.6)
where P
0
is the scalar seismic potency of the event, Δε
s
is the static strain drop assumed
here to be 10
−4
and c = 2.283. The P
0
value of each event is calculated from the
empirical potency-magnitude scaling relation of Ben-Zion and Zhu (2002),

log
10
P
0
= 0.06M
L
2
+ 0.98M
L
− 4.87 , (1.7)
with M
L
being the local magnitude of the event and P
0
is in km
2
times cm. Figure 6 of
Ben-Zion (2008) shows r values as a function of M
L
calculated with the same
procedure. Since aftershocks extend beyond the mainshock rupture area, we multiply
the estimated rupture radius by a factor L that varies between 1 and 50 and use the result
to be the radius of a circular aftershock zone. In the analysis below we divide the
productivities of the stacked aftershock sequences by the assumed L factor. The results
20
show that L = 20 provides normalized productivities that remain similar when larger L
values are used (Figure 1.5). Dynamically triggered aftershocks may extend to larger
distances (e.g. Felzer and Brodsky 2006; Hill and Prejean 2007), but the considered
areas with L = 20 include the bulk of classical aftershocks and reduce mixing
aftershocks with background seismicity. Since the stacked data have more events than
individual aftershock sequences, we use here discrete time intervals of 0.5 day (half the
size that was used in section 1.2.1). This procedure produces 28-196 aftershock
sequences in each of the five regions (middle panels in Figure 1.4). The stacked data
associated with all sequences in a given region (solid circles in the middle panels)
highlight the common characteristics of the aftershock properties in the region.
21

Figure 1.5. (a) Estimated productivity values of stacked aftershocks versus the L factor that scales the size
of the considered aftershock areas for the different regions of Figure 1.1. (b) Estimated productivity
values for the different regions normalized by the used L factor. The results are stable for 15 < L < 50.
The value L = 20 (dashed vertical line) is used for the results in Figures 1.4 and 1.6 and Table 1.1.

22
We fit the stacked aftershock sequences with the Omori-Utsu law using the same
method as discussed in section 1.2.1 with c = 0, p = 1 and t
0
= 5. The fitted results are
shown in the right panels of Figure 1.4. The regional aftershock productivities decrease
in the following order: the Landers and Hector-Mine area, the San Bernardino region,
the Ventura Basin, the Imperial Valley, and the Coso region. To verify that the
comparison of the K values for the different regions is not strongly affected by the
different numbers of aftershock sequences used in the different regions, we perform
similar analysis employing for each region 25 randomly selected sequences from the
available data. We calculate the parameters of stacked aftershock sequences using 500
realizations of 25 randomly selected sequences. The mean results differ from those of
Figure 1.4 only by 0.4-1.3% and the order of the productivity values remains the same.
Figure 1.6a shows the stacked aftershock sequences and lines corresponding to the
Omori-Utsu curves in discrete form for the five different regions using the parameters
associated with the fits obtained in Figure 1.4. The relations between the estimated heat
flow values and normalized aftershock productivities (K values from Figure 1.4 divided
by L = 20) for the five aftershock regions are shown in Figure 1.6b. The observed
correlations between aftershock productivity, heat flow and seismogenic zone that
includes sedimentary rocks are compatible overall with the damage model predictions
of Ben-Zion and Lyakhovsky (2006).
23

Figure 1.6. (a) Numbers of stacked aftershocks per 0.5 day (circles) and fitted curves (lines) based on the
discrete Omori-Utsu relation (equation 1.1) with c = 0, p = 1 and t
0
= 5 day for the five selected regions
(colored circles and curves as indicated by the legend). (b) Relations between aftershock productivity and
heat flow values in the five examined regions. The horizontal error bars represent 95% confidence
intervals of the heat flow data. The vertical error bars represent the standard deviations of the
productivities.

1.3 Discussion
The results from both the analysis of large aftershock sequences in two different
regions and analysis of many stacked sequences in five different regions indicate clearly
(Figures 1.3b and 1.6b) that there is a general inverse relationship between the mean
heat flow and aftershocks productivity in a region. These results are compatible with the
previous observational findings of Mogi (1967), Kisslinger and Jones (1991) and other
studies summarized by Utsu et al. (1995). Among the five examined regions, the
Landers and Hector-Mine area, the Imperial Valley and the Coso region have seismicity
24
that is generally below the sediments. The corresponding points for these three regions
in the productivity versus heat flow plot of Figure 1.6b fall nearly on a straight line.
This indicates that for seismogenic zones with similar lithology (e.g. crystalline rocks),
the aftershock productivity is primarily determined by the heat flow (or more correctly
the temperature field). The Landers and Hector-Mine area, the San Bernardino region
and the Ventura Basin have similar heat flow levels but considerably different thickness
of sedimentary covers. The aftershock productivities are lower in the latter two regions,
where the sedimentary covers extend partially (in the San Bernardino region) or almost
fully (in the Ventura Basin) into the seismogenic zone. The above observations are
consistent with the expectations based on the damage rheology model of Ben-Zion and
Lyakhovsky (2006).
As mentioned in the introduction, the seismic coupling coefficient χ is expected to be
inversely proportional to the material parameter R (equation 1.2), and to be directly
proportional to the aftershocks productivity in a region (equations 1.3a and 1.3b). Since
the maximum event size and seismic hazard in a region increase with increasing seismic
coupling, it is important to develop methodologies for estimating the χ values in
different regions. Based on the damage rheology model, the largest value of χ in the
study area is expected to be in the Landers and Hector-Mine region, and the lowest
values in the Imperial Valley and the Coso region. This is consistent with the facts that
the largest earthquake in the analyzed 1984-2002 catalog (the M7.3 1992 Landers
earthquake) is in region B, and the smallest mainshocks are in regions A and E
(although the relative short duration of the catalog compared to recurrence times of
large earthquakes precludes drawing conclusions from this alone).
25
A direct estimate of the χ value that characterizes a given fault system requires a
comparison between the seismic release and geodetically observed deformation in the
region. Lohman and McGuire (2007) compared geodetic data with earthquake
properties within the Salton Trough area of southern California (green rectangle in
region A of Figure 1.1) and estimated the χ value for their study area to be 0.2. Their
results on earthquake sequences in the area are consistent with the expected correlations
of Ben-Zion and Lyakhovsky (2006). Most seismically active regions (especially
outside the US and Japan) do not have geodetic networks that operate over sufficiently
long time (e.g. > 10 yr) and length (e.g. 10 km) scales to allow direct estimates of the
seismic coupling coefficient. An ability to estimate χ and the material parameter R,
through analysis of aftershocks productivity of the type done in this paper, can therefore
be useful for seismic hazard assessment and various other studies. This can be done as
follows. From equations (1.2) and (1.3b) the ratios of the R and χ values in different
regions can be written as,
, (1.8)
.
Combining these two relations leads to
(1.9)
In the absence of any direct information on the seismic coupling (or R), one may
assume a χ value for one region and then derive corresponding relative values using the
observed ratios of the aftershock productivities for other regions. As mentioned, based
26
on the damage model results of Ben-Zion and Lyakhovsky (2006) we expect the
Landers and Hector-Mine area to have relatively low R value and relatively high χ
value. Assuming that area B of Figure 1.1 is characterized by χ = 0.85, we can derive
from the observed normalized K values of Figure 1.6 and equations (1.9) and (1.8) the χ
and R values for all five study areas. The results are summarized in the left χ and R
entries of Table 1. 1.
Region Area (km
2
) K/L χ R
A 10428 1.8 0.70, 0.66 0.42, 0.51
B 7274 4.3 0.85
©
, 0.82 0.18, 0.21
C 5548 2.9 0.79, 0.76 0.26, 0.31
D 5740 2.3 0.75, 0.71 0.33, 0.40
E 3510 0.54 0.42, 0.37 1.4, 1.7
Green box in A 233 0.23 0.23, 0.20* 3.3, 4.0

Table 1.1 Estimated χ and R values for the different regions of Figure 1.1 based on aftershock
productivities K in areas that scale with the mainshock magnitude and numerical factor L = 20. The χ
value denoted by
©
is an educated guess for the Landers and Hector-Mine area (region B). The other left
entries of χ and R are based on this value. The χ value denoted by * is from Lohman and McGuire (2007).
All other right entries of χ and R are based on this value.

See text for more explanation.

In our case, we can also use the estimated value χ = 0.2 of Lohman and McGuire
(2007), along with estimated K value for the same area, and then derive from those
results χ and R values for all other regions. Applying the stacking approach of section
1.2.2 to the green box in area A of Figure 1.1, we obtain a normalized K = 0.23 from
data associated with that region. The χ and R values for the different regions based on
χ = 0.2 and K = 0.23 for the green box in Figure 1.1 are summarized in the right χ and
R entries of Table 1.1. The obtained χ = 0.82 for area B of Figure 1.1 is close to our
assumed 0.85 value.
27
Observations of earthquake sequences can be done relatively easily. The damage
model results of Ben-Zion and Lyakhovsky (2006) and forgoing procedure provide an
effective tool for deriving from properties of aftershock sequences fundamental
parameters that characterize the rheology and seismic potential of a region. Our analysis
employs a simple procedure associated with pre-selection of regions, based on
supplementary knowledge on clear differences in assumed controlling regional
conditions, and mapping the differences of aftershocks rates in given magnitude ranges
to the productivity parameter. More sophisticated analyses (e.g. using smoothly-varying
adjustable selection of regions, assimilation of more detailed information on heat flow,
lithology and other relevant information, joint inversion of multiple parameters, etc.)
may lead to better estimates of the seismic coupling coefficient. This is left for future
work.

28
Chapter 1 References

Bates, D.M. & Watts, D.G., 1988. Nonlinear Regression Analysis and Its Applications,
Wiley.

Ben-Zion, Y., 2003. Key Formulas in Earthquake Seismology, International handbook
of earthquake and engineering seismology, Part B, 1857-1875, Academic Press.

Ben-Zion, Y., 2008. Collective Behavior of Earthquakes and Faults: Continuum-
Discrete Transitions, Progressive Evolutionary Changes and Different Dynamic
Regimes, Rev. Geophysics, 46, RG4006, doi:10.1029/2008RG000260.

Ben-Zion, Y. & Lyakhovsky, V., 2006. Analysis of aftershocks in a lithospheric model
with seismogenic zone governed by damage rheology, Geophys. J. Int., 165, 197-
210.

Ben-Zion, Y. & Zhu, L., 2002. Potency-magnitude scaling relations for Southern
California earthquakes with 1.0 < M
L
< 7.0, Geophys. J. Int.,148, F1–F5.

Davis, S.D. & Frohlich, C., 1991. Single-link cluster analysis of earthquake aftershocks:
decay laws and regional variations, J. geophys. Res., 96, 6335–6350.

Enescu, B., Mori, J. & Miyazawa, M., 2007. Quantifying early aftershock activity of the
2004 mid-Niigata Prefecture earthquake (Mw6.6), J. Geophys. Res., 112, B04310,
doi:10.1029/2006JB004629.

Felzer, K.R., Abercrombie, R.E. & Ekstrom, G., 2004. A common origin for
aftershocks, foreshocks and multiplets, Bull. Seismol. Soc. Am., 94, 88-98.

Felzer, K.R. & Brodsky, E. E., 2006. Decay of aftershock density with distance
indicates triggering by dynamic stress, Nature, 441, 735-738.

Gross, S.J. & Kisslinger, C., 1994. Tests of models of aftershock rate decay. Bull.
Seismolo. Soc. Am., 84, 1571-1579.

Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y. & Lockner, D., 2004. A Visco-
elastic damage model with applications to stable and unstable fracturing, Geophys.
J. Int., 159, 1155–1165, doi: 10.1111/j.1365-246X.2004.02452.x.

Helmstetter, A., Kagan, Y.Y. & Jackson, D.D., 2005. Importance of small earthquakes
for stress transfers and earthquake triggering, J. Geophys. Res., 110, B05S08,
doi:10.1029/2004JB003286.

29
Hill, D.P. & Prejean, S.G., 2007. Dynamic Triggering, in Earthquake Seismology -
Treatise on Geophysics ed. H. Kanamori and G. Schubert, V4, 257-291, Elsevier,
Amsterdam.

Kisslinger, C. & Jones, L.M., 1991. Properties of aftershock sequences in southern
California. J. Geophys. Res., 96, 11,947-11,958.

Lohman, R.B. & McGuire, J.J., 2007. Earthquake swarms driven by aseismic creep in
the Salton Trough, California. J. Geophys. Res., 112, B04405,
doi:10.1029/2006JB004596.

Lyakhovsky, V., Ben-Zion, Y. & Agnon, A., 1997. Distributed Damage, Faulting, and
Friction, J. Geophys. Res., 102, 27635-27649.

Lyakhovsky, V. & Ben-Zion Y., 2008. Scaling relations of earthquakes and aseismic
deformation in a damage rheology model, Geophys. J. Int., 172, 651-662, doi:
10.1111/j.1365-246X.2007.03652.x.

McGuire, J.J., Boettcher, M.S. & Jordan, T.H., 2005. Foreshock sequences and short-
term earthquake predictability on East Pacific Rise Transform Faults, Nature, 434,
457–461.

Mogi, K., 1967. Earthquakes and fractures. Tectonophysics, 5, 35-55.

Nauteau, C., Shebalin, P. & Holschneider, M., 2002. Temporal limits of the power law
aftershock decay rate, J. Geophys. Res., 107, B12, 2359,
doi:10.1029/2002JB001868.

Omori, F., 1894. On the aftershocks of earthquakes, J. Coll. Sci., Imp. Univ. Tokyo, 7,
111-200.

Otsuka, M., 1985. Studies on aftershock sequences-Part 1. Physical interpretation of
Omori formula, Sci. Rep. Shimabara Earthq. Volcano Obs., 12, 11-20 (in Japanese).

Peng, Z., Vidale, J.E. & Houston, H., 2006. Anomalous early aftershock decay rates of
the 2004 M6 Parkfield earthquake, Geophys. Res. Lett., 33, L17307,
doi:10.1029/2006GL026744.

Shearer, P., Hauksson E. & Lin, G., 2005. Southern California hypocenter relocation
with waveform cross-correlation, Part 2: Results using source-specific station terms
and cluster analysis, Bull. Seismol. Soc. Am., 95, 904-915, doi:
10.1785/0120040168.

Utsu, Y., Ogata, Y. & Matsu’uara, R.S., 1995. The centenary of the Omori formula for
a decay law of aftershock activity, J. Phys. Earth, 43, 1–33.
30

Wiemer, S. & Wyss, M., 2000. Minimum Magnitude of Completeness in Earthquake
Catalogs: Examples from Alaska, the western United States and Japan, Bull.
Seismol. Soc. Am., 90, 4, 859-869.

31
Chapter 2

Variations of strain-drops of aftershocks of the 1999
İzmit and Düzce earthquakes around the Karadere-
Düzce branch of the North Anatolian Fault

SUMMARY
We estimate the strain-drops of 7498 aftershocks of the 1999 İzmit and Düzce
earthquakes using P waveforms recorded by a local seismic array along the Karadere-
Düzce branch of the North Anatolian fault in the 6 months following the İzmit
mainshock, and evaluate the site effects of the various recording stations. The method is
associated with separation of source, travel-time and station spectral terms and stacking
results at several stages to enhance the signal-to-noise ratio. The strain-drops are
obtained by fitting iteratively the separated source spectra of 201 nearest neighboring
events in different amplitude bins to the ω
-2
source spectral model. The obtained strain-
drops generally increase with depth between 3 and 10 km and remain approximately
constant for the deeper section. A local shallow patch of seismicity north of the
Karadere segment has relatively low strain-drop values. Along the relatively straight
Karadere segment the strain-drops are lower than along the geometrically complex
bounding regions. In some sections the range of strain-drop values decay with time
from the İzmit and Düzce mainshocks, while in others the values fluctuate within about
32
constant ranges. The observed spatial variations of strain-drops can be explained by
increasing normal stress with depth along with the degree of geometrical complexity of
different fault sections and the expected slip deficit at different depth sections. The
seismic energy in the separated station spectra decreases overall with distance from the
rupture zone and the spectra of various stations have 3 peaks at 6, 14 and 25 Hz. The
spectral peak at 6 Hz is also observed in trapped waves studies and may be related to
the damaged fault zone layer.

2.1 Introduction
Observed seismograms result from a convolution of source, propagation path, and
site-instrument effects. Many useful source, path, and site parameters can be obtained
from recorded seismic waveforms after properly separating the various effects. One of
the most important source parameters that can be computed through seismological
methods is the reduction of the average elastic strain (e) or stress (s) over the earthquake
rupture area. The strain- or stress-drops can be estimated from the shape of the
earthquake source spectrum (i.e. corner frequency f
c
) by assuming various source
models (e.g. Brune 1970; Madariaga 1976; Aki & Richards 2002; Ben-Zion 2003).
These parameters provide fundamental information on the physics of earthquake rupture
and spatio-temporal evolutions of stress-strain conditions around earthquake faults.
Many methods have been developed in the past to isolate source signals from the
observed seismograms. The most commonly used one is referred to as the empirical
Green function (EGF) method (e.g. Berckhemer 1962; Mueller 1985; Hough 1997). If
waveforms of two neighboring events with large magnitude difference are recorded by a
33
common station, the source signal of the smaller event could be approximately treated
as a delta function. Since the radiations from both events share very similar path and
site responses, the relative source function of the larger event can be obtained by
deconvolving its seismic record with that of the smaller event. However, the application
of the EGF method is somewhat limited because it requires two events with close
locations and with large size difference. In addition, the deconvolution operator may
introduce some artificial signals if the noise level of the smaller event is high. Another
group of methods take advantage of the fact that multiple events are recorded by
multiple stations, and separate the relative source, path and station terms by solving a
series of linear equations (e.g., Andrews 1986; Boatwright et al. 1991). However, the
obtained results may not be stable due to the existence of random noise. Warren &
Shearer (2002) introduced a stacking technique to enhance the common features of
signals and decrease the effect of random noise. Such method has been applied to obtain
reliable source properties of small earthquakes at many regions (Prieto et al. 2004,
Shearer et al. 2006; Allmann & Shearer 2007).
The above developments notwithstanding, basic aspects of source and sites effects
remain unclear, primarily because of the limited direct information and considerable
uncertainties at both ends of the source-observation system. Several sources of
uncertainties contribute to the large scatter and discrepancies of the available results.
First, estimations of the corner frequency, stress/strain-drop, rupture velocity and
seismic energy depend on the high-frequency part of the seismic spectra, where
correcting for attenuation and other path effects may be difficult. Second, previous
studies have been done either in one location with a single or a few stations (e.g.
34
Abercrombie 1995), at a large region with a mixture of tectonic environments (e.g.
Shearer et al. 2006), or using combinations of data associated with different studies
(e.g. Ide & Beroza 2001). Systematic analysis of large seismic waveform data recorded
by multiple stations in nearby regions can provide better constraints on earthquake
source properties and site effects.
In this study we apply an iterative stacking technique to separate source, propagation
and site effects in a large seismic data set recorded along and around the Karadere-
Düzce branch of the North Anatolian Fault (NAF) during the 1999 İzmit and Düzce
earthquake sequences in Turkey. The earthquakes (Figure 2.1) are highly concentrated
in space (within ~ 10-km-wide 100-km-long region around the Karadere-Düzce branch
of the NAF) and time (within ~ 6 months). Moreover, the data were recorded over a
period that includes both the early post-mainshock and final pre-mainshock stages of
large earthquake cycles (in relation to the İzmit and Düzce events, respectively). We
focus on analysis of P waveforms and estimate strain-drops of about 7500 earthquakes
and site effects of stations with various distances from the rupture zones of the İzmit
and Düzce mainshocks. The derived strain-drops are generally higher along fault
sections with larger geometrical complexity and they increase with depth. The seismic
energy associated with site effects includes spectral peaks around 6, 14 and 25 Hz, and
is largest at a station that was installed within the rupture zone of the İzmit earthquake.

35
2. 2 Data and Preprocessing
2.2.1 Data
A temporary PASSCAL network of 10 short-period seismic stations was deployed
along and around the Karadere-Düzce branch of NAF a week after the August 17, 1999,
Mw7.4, İzmit, Turkey, earthquake (Figure 2.1). All stations had REFTEK recorders and
three-component L22 velocity sensors with a sampling frequency of 100 Hz (Seeber et
al. 2000; Ben-Zion et al. 2003). This network operated for 6 months and recorded
approximately 26000 events, including the November 12, 1999, Mw7.1, Düzce
earthquake. In this study, we use a subset of 18556 events that are within 30 km of the
network (Peng & Ben-Zion 2004). The employed events have magnitudes in the range
between -0.1 and 3.9. The horizontal location errors are generally less than 1 km around
the center of the network and 1-2 km near the margins. The vertical location errors are
somewhat larger. Additional details on the seismic experiment and data set are given by
Seeber et al. (2000) and Ben-Zion et al. (2003).
36

Figure 2.1. Hypocentral distribution of 18556 earthquakes (green points) along the Karadere-Düzce
branch of the North Anatolian Fault. The strain-drops of 7498 events (red points) are estimated in this
study. The gray lines denote observed rupture traces of the İzmit and Düzce mainshocks and the black
triangles are the ten seismic stations. Focal mechanisms of both the İzmit and Düzce mainshocks point to
their epicenters. The events inside the black circle are used with 201 nearest-neighbors to illustrate the
separation of the source, station and travel-time terms (Figure 3a), to stack the separated source term
(Figure 3b) and to obtain the best fitting strain-drop (Figure 4). The inset map shows the location of the
study area (black rectangle) in a large-scale regional map.

2.2.2 Preprocessing
We analyze the P-waves recorded on the vertical component of the short-period
instruments. The preprocessing steps are as follows. We remove the mean and trend
from each seismogram, correct the instrumental response, and convert the waveforms
from velocity to displacement records. We also remove traces that are clearly clipped.
These analysis steps are done using Matlab functions in the CORAL package for
seismic waveform analysis (Creager 1997). We take the signal window to be 0.28 s
before and 1.0 s after the picked P arrival, and take a background noise window to be
37
with the same length before the picked P arrival. The displacement spectra for the noise
and signal windows are computed using a multitaper algorithm (Thomson 1982; Park et
al. 1987). We require the signal to noise ratio (SNR) of the obtained spectra in the
frequency band 2.35-20 Hz to be larger than 3.0. We also require for further analysis the
minimum number of qualified observed waveforms for an event to be larger than 3.
After the data preprocessing, we have 48,580 P-wave spectra from 7,498 events.
Examples of waveforms and corresponding signal and noise spectra for one event at
several stations are shown in Figure 2.2.

Figure 2.2. Example results based on a Jan 20, 2000, M
L
2.0 earthquake recorded by four stations (BV,
VO, FP, and FI). The left panels show the velocity waveforms around the P-wave window. The station
name, component and epicentral distance are marked on the left. The right panels show the log
10

displacement spectra for both signal (solid line) and noise (dashed line).
38
2.3 Method
We use a method that generally follows the procedure of Shearer et al. (2006) to
estimate potency values and strain-drops of earthquakes in a large regional data set. The
method involves spectral stacking at several steps to increase the stability of the results
and is associated with three major analysis stages. In the first stage we separate
iteratively the source, travel-time and station spectra, using for each earthquake the
observed data of its N nearest neighbor events (here N = 200). In the second stage we
bin the separated source terms in different amplitude ranges and stack the source spectra
in each bin. In the third stage we fit the stacked source spectra of each bin with a
theoretical source model (we assume

ω
−2
spectral decay) and obtain best-fitting
earthquake strain-drops.
Our method differs from the studies of Shearer et al. (2006) and related works (e.g.
Prieto et al. 2004; Allmann & Shearer 2007) primarily by using in the analysis strain-
based quantities (potency and strain-drops) rather than the more traditional stress-based
quantities (moments and stress drops). This makes the results independent of assumed
rigidity values in the source regions. As discussed by Ben-Zion (2003, 2008), the
rigidity values vary rapidly and are hence poorly defined in the space-time windows
associated with earthquakes. Moreover, the rigidity values in the source volumes do not
affect the seismic radiation in the surrounding elastic solid and can thus be assigned
arbitrarily. For these reasons we prefer to use parameters associated with strain-based
quantities. Additional details on the various analysis procedures are given below.

39
2.3.1 Iterative Separation of Source, Station and Travel-Time Spectra
In the log-frequency domain, the observed displacement spectrum can be expressed
as
ij j ij i ij
r + s + t + e = d , (2.1)
where d
ij
is the observed displacement spectrum of the i-th event recorded by the j-th
station, e
i
is the source term of the i-th event, s
j
is the station term of the j-th station, t
ij
is
the travel-time term for the event-station pair, and r
ij
is a residual term. Assuming that
the attenuation is associated with a constant value, t
ij
can be discretized into a sequence
of k intervals that cover the range of all possible travel time terms (Shearer et al. 2006).
In this procedure, we consider spectral data in a frequency range between 2.35 and 45
Hz that is wider than the range used in the subsequent analysis steps. We discretize the
travel-time term t
ij
at 1.0 s increments, and set the range of acceptable residuals r
ij
to be
within [-1.0 - 1.0] nm/Hz. Each event is analyzed jointly with its nearest N = 200
events. The iterative separation of the different terms of Eq. 1 for each group of N + 1
events is done as follows:
1. We set the initial source, travel-time and residual terms at the first iteration to be
zero. We solve for the initial station term as the stacked value of the observed
displacement spectra,

s
j
1
=
1
m
d
ij
i=1
m
∑
, (2.2)
where m ≤ N is the total number of events recorded by the j-th station.
2. We solve for the source term at the p-th iteration (p = 2,3,…) using
40

e
i
p ( )
=
1
n
d
ij
− s
j
p−1 ( )
− t
k ( ) ij
p−1 ( )
− r
ij
p−1 ( )
⎡
⎣
⎢
⎤
⎦
⎥
j=1
n
∑
, (2.3)
where n is the total number of stations which recorded the i-th event and the superscript
k denotes the appropriate discrete travel time for various source-station pairs ij.
3. We solve for the travel-time term at the p-th iteration (p = 2,3,…) using
( ) ( ) ( ) ( )
[ ]
∑
− −
− − −
1 1
1
p
ij
p
j
p
i ij
k
p
k
r s e d
N
= t , (2.4)
where N
k
is the total number of event-station pairs within the k-th travel-time interval (k
= 1,2,…).
4. We solve for the station term at the p-th iteration (p = 2,3,…) as

s
j
p ( )
=
1
m
d
ij
− e
i
p ( )
− t
k ( ) ij
p ( )
− r
ij
p−1 ( )
⎡
⎣
⎢
⎤
⎦
⎥
i=1
m
∑
. (2.5)
5. We solve for the residual term at the p-th iteration (p = 2,3,…) as
( ) ( ) ( )
( )
( ) p
ij k
p
j
p
i ij
p
ij
t s e d = r − − − . (2.6)
If the mean value of
ij
r is outside the acceptable residual boundaries, the corresponding
waveform will be dropped from subsequent iterations. If the corresponding event has
less than three waveforms, it will be dropped from further analysis.
6. We calculate the following sums of the absolute difference values:
( ) ( )
∑
−
−
i
p
i
p
i
e e
1
,
( ) ( )
∑
−
−
j
p
j
p
j
s s
1
,
( ) ( )
∑
−
−
i
p
ij k
p
ij k
t t
1
) ( ) (
.
The iteration stops when each sum value is less than a
very small number (e.g. 1×10
-4
) and otherwise it returns to step 2.
As in the study of Shearer et al. (2006), the analysis steps discussed above usually
converge and terminate in 2-3 iterations. An illustration of the separated source, station,
41
travel-time and residual terms of 201 nearest neighboring events, with locations inside
the black circle of Figure 2.1, are shown in Figure 2.3a.

Figure 2.3a. Illustration of the source, station and travel-time terms separation for 201 nearest
neighboring events.

2.3.2 Stacking of Source Terms with Similar Amplitudes
The separated 201 source spectra obtained by the procedure of section 2.3.1 can have
large fluctuations (left panel in Figure 2.3b) that are at least partially due to errors from
random noise. To capture the common features of the neighborhoods of earthquake
sources that are analyzed jointly, we follow Shearer et al. (2006) and stack the source
spectra over 0.2 logarithmic amplitude bins. The shapes of the obtained stacked source
42
spectra at the different amplitude bins are relatively smooth (right panel in Figure 2.3b).
As expected from the Gutenberg-Richter frequency-size statistics of earthquakes, the
spectral shapes for bins at lower amplitudes are stacked from a larger number of
individual spectra than those at higher amplitude bins. This explains why the stacked
source spectra for the larger events still have some shape irregularities at high
frequency. However, larger events generally have higher SNR and their corner
frequencies are better contained within the frequency range (2.35 - 20 Hz) analyzed in
our study.

Figure 2.3b. Illustration of separated source spectra and the corresponding stacking. The left panel shows
the separated source spectra for 201 nearest neighboring events. The right panel shows the stacked source
spectra with 0.2 logarithmic amplitude increments.
43
2.3.3 Fitting Source Spectra with Theoretical Model
In this analysis stage the binned source spectra are fitted with a theoretical source
model, and a near-source common EGF is simultaneously removed from the results. We
use the theoretical source model of Brune (1970)

u f
( )
=
Ω
0
1+ f / f
c
( )
n
, (2.7)
where u(f) is the displacement spectra, f
c
is the corner frequency and Ω
0
is the zero
frequency asymptote (proportional to the scalar seismic potency P
0
and scalar seismic
moment M
0
= m P
0
with m being a nominal rigidity at the source). Since the binned
source spectra are relatively flat at the low frequency range (Figure 2.3b), we use the
mean amplitude of the first four points (2.35 - 3.53 Hz) of the source displacement
spectra to estimate Ω
0
.
The corner frequency f
c
can be related to the scalar seismic moment M
0
and seismic
potency
0
P by assuming a circular rupture model (Madariaga 1976; Abercrombie 1995;
Prieto et al. 2004; Shearer et al. 2006; Allmann & Shearer 2007) with a constant rupture
velocity of 0.9 times the shear wave velocity β,

f
c
=
0.42β
( M
0
/Δσ )
1/3
=
0.42β
(P
0
/Δε)
1/3
, (2.8)
where Δσ and Δε are the stress- and strain-drops of the earthquake, respectively.
Since the zero frequency asymptote Ω
0
obtained from the separated source spectrum
is a relative value, it does not provide information on the corresponding absolute
seismic potency. To calibrate the average potency in each amplitude bin with the
corresponding average local magnitude, we use the empirical potency-magnitude
44
equation of Ben-Zion & Zhu (2002)

P
0
= 10
0.0612M
L
2
+ 0.988M
L
− 4.87
.
(2.9)
Using Eqs (2.7)-(2.9), we fit the source spectra with a procedure involving only
variations of the assumed strain-drops for the groups of analyzed events (Figure 2.4).
This is done as follows:
1. We set a range of possible values of strain-drops 3×10
-6
≤ Δε ≤ 1×10
-2
(Figure
2.4a). Assuming a nominal rigidity 30 GPa, this corresponds to stress drops in the range
0.09-300 MPa. The range of allowed log
10
strain-drops is divided to intervals of 0.0035,
giving a total number of 1000 strain-drops values. For each possible strain-drop value,
we compute the theoretical source spectra for each bin using the ω
-2
model (n = 2 in Eq.
7).
2. We calculate the difference between the theoretical source spectrum for a given
strain-drop value and stacked source spectra at each amplitude bin (Figure 2.4b). We
then stack the differences of all amplitude bins to obtain a single common EGF (Figure
2.4c).
3. We compute the Root Mean Square (RMS) misfit, which is the sum of differences
between the theoretical source spectra and the EGF-corrected source spectra for all
amplitude bins over the employed frequency band (2.35 - 20 Hz). The value of Δε that
provides the minimum RMS misfit (Figures 2.4a and 2.4d) is selected as a best-fitting
strain-drop for the group of events.
45

Figure 2.4. Illustration of the strain-drop fitting for an event and its 200 nearest-neighbors. (a) The fitting
RMS curve in the log10 strain-drops ranges between -5.0 and -2.0. (b) The separated source spectra are
binned with 0.2 increments in magnitude (solid color lines). The theoretical source spectra with the strain-
drop as unknown parameter are shown in blue dashed lines. The black dashed line marks the upper limit
of the fitting frequency. (c) The difference between the theoretical and separated source spectra in each
bin (color curves) and the stacked difference (bold green curve) giving the common EGF. (d) The EGF-
corrected source spectra (red curves) and the theoretical source spectra (blue dashed curves). The best-
fitting strain-drop corresponds to the minimum fitting RMS value in (a).

2.4 Results
For each of the 7,498 events that satisfy the discussed quality criteria (section 2.2),
we use the above procedure to derive potency and strain-drop values with both constant
and depth-variable velocity models. About 95% of the derived log
10
strain-drops are in
the range from -4.3 to -2.6 (corresponding to 1.5-75 MPa stress-drops assuming a
46
nominal rigidity of 30 GPa), with a mean log
10
value of -3.5 (9.5 MPa). Details of the
obtained results are presented below.

2.4.1 Spatial Variations of Strain-drops

Figure 2.5a. A map of the best-fitting log
10
strain-drops for 7498 aftershocks of the 1999 İzmit and Düzce
earthquakes along the Karadere-Düzce branch of the North Anatolian Fault. The results are colored in
equal increments of log
10
Δε with color bar at the upper-left corner. Four sub-regions are separated by
dashed lines with labels I, II, III and IV at the up-left corners. Two sub regions in section II with higher
strain-drops are marked with circles and labeled as IIA and IIB. The black triangles represent the ten
stations. The seismicity projected along line AB of N71°E (parallel to the fault strike of the Karadere
segment) is shown in Figures 2.6 and 2.7.

Figure 2.5a shows the best-fitting strain-drops for the analyzed 7498 events (red dots
in Figure 2.1). The results are derived by assuming a constant shear wave velocity b =
3.5 km/s (i.e. a constant P-wave velocity α = 6 km/s, and α /b = √3). Figure 2.5b
illustrates the goodness of obtained source spectra and strain-drop fitting using ten M
L
≈
2.1 events at different locations inside the study region (marked as cyan circles in
47
Figure 5a). For each event, we plot the EGF-corrected stacked source spectra of
earthquakes with different spectral amplitude bins belonging to the N = 201
neighborhood of the selected events (panels in Figure 2.5b), along with the theoretical
source spectra of the best-fitting strain-drop (dashed blue lines) and the common EGF
(dashed green line). As in Figure 2.4b, the stacked source spectra are relatively smooth
and generally match well the theoretical curves associated with the best-fitting strain-
drops. As shown in Figure 2.5a, the strain-drops vary smoothly for nearby events, as
expected from the smearing effects of the 201 nearest neighboring approach. However,
there are clear variations of strain-drops along the different fault traces in the study area.
Based on the overall observed patterns of the strain-drops and aftershock locations, we
divide the examined region into four sections and discuss the patterns in each section
separately.

48

Figure 2.5b. The EGF-corrected stacked source spectra (red curves) for each of the ten M
L
≈ 2.1 events
marked in Figure 2.5a with its 200 nearest-neighbors and the theoretical source spectra (blue dash curves)
associated with the best-fitting strain-drop. In each panel, the number at the top-right corner corresponds
to the appropriate number within the cyan circle in Figure 2.5a, and the best-fitting log
10
strain-drop is
marked at the bottom-left corner. For a better comparison of the spectral curves, the various spectral
amplitudes are shifted according to the median event magnitude in each panel. Other symbols and
notations are the same as in Figure 2.4.

Section I is between longitudes 30º20´ and 30º40´, where the strike of the surface
trace of the NAF is along the E-W direction. Most aftershocks in this section were
generated by the İzmit mainshock and they are associated generally with depth larger
than 5 km (Figure 2.6). The strain-drops in this section are relatively high. Section II is
between longitudes 30º40´ and 31º, where the fault strike of the Karadere segment is
N71
o
E. Most aftershocks were generated by the İzmit mainshock and they form a linear
band (line AB), which is parallel to the surface trace of the Karadere segment with a 2-3
km offset to the northwest as the Karadere segment steeply dips to the NW. Additional
49
aftershocks occurred on the SE side of station CH and NE side of station GE, which are
far from the surface ruptures of the İzmit mainshock. Numerous aftershocks in this
section occurred in the top few km of the crust (Figure 2.6), and both the strain-drops
and magnitude of the events are relatively low. Two small areas in this section (IIA and
IIB in Figure 2.5a) have clusters of events off the main fault with depths between 5-12
km and relatively high strain-drops (IIA and IIB in Figure 2.6). Section III is defined as
the region between longitudes 31º and 31º12´, and is associated with an overlap
between the eastern end of the İzmit aftershock zone and the western part of the Düzce
aftershock zone. The events in this section are generally deeper than 5 km (Figure 2.6)
and they are associated with relatively high strain-drops. Section IV is between
longitudes 31º12´ and 31º20´, where the seismicity was primarily triggered by the
Düzce mainshock (Figure 2.7). This section has the highest average strain-drops among
all sections. Table 1 summarizes statistical features of the strain-drops in the different
sections.
Section I II III IV
Mean log
10
Δε -3.29±0.31 -3.85±0.32 -3.40±0.25 -3.04±0.09
Event Number 1528 3289 1931 750
Mean magnitude 1.73±0.56 1.37±0.63 1.98±0.61 2.29±0.59
Table 2.1. Statistical features of the four divided sections.

To illustrate further the spatial variations of earthquake strain-drops along the
Karadere-Düzce branch of the NAF, we project the events onto a vertical plane
intersecting the surface at line AB of Figure 2.5a, which is parallel to the strike of the
fault trace in the central Karadere segment (section II). The distribution of the strain-
drops along the vertical plane is shown in Figure 2.6. To examine the effects of the
assumed constant β = 3.5 km/s (red line in the left panel of Figure 2.6), we also compute
50
strain-drops with a 1-D velocity model for this region given by black line in the left
panel of Figure 2.6 (Seeber et al. 2000). The average values of the strain-drops versus
depth derived using both velocity models are plotted in the right panel of Figure 2.6.
For better spatial resolution, we examine the average depth variation of strain-drops of
all events (solid lines) as well as the events occurring only inside section II (dashed
lines). The results obtained with the two different velocity models are very similar
except for the shallow section. With the 1-D velocity model, the average strain-drops
(black lines) in the top 4 km have higher values than those obtained with the constant
velocity model (red lines). With the constant velocity model, the average log
10
strain-
drop of all events increases with depth from -4.2 around 3 km to -3.3 at 10 km depth
(corresponding to nominal stress-drops of 1.9 and 15.0 MPa, respectively). The strain-
drops inside section II have weak depth dependency with the constant velocity model,
and nearly no depth dependency with the 1-D velocity model.
The top panel of Figure 2.6 shows variations of the average strain-drops with distance
along the line AB of Figure 2.5a, using results obtained with both velocity models. For
better spatial resolution, we present average strain-drops for all events (solid lines) as
well as for events in the depth range 5-15 km (dashed lines) where most aftershocks
occurred. As mentioned earlier, the strain-drops inside section II are lower compared to
those in the other sections. Most of the earthquakes with depth less than 5 km occurred
in the top-right portion of section II (NE of station GE in Figure 2.1), and most of them
have relatively small strain-drop values. The differences somewhat decrease for results
based on the 1-D velocity model (black lines), but the overall pattern remains the same.
51

Figure 2.6. Depth variations of log
10
strain-drops along the strike of the seismic lineation in section II
(line AB in Figure 5a). The results are colored in equal increments of log
10
Δε with color bar at the
bottom. The panel at the left side shows the constant (red line) and 1D velocity (black line) velocity
models. Four different fault sections are separated by vertical dashed lines and the horizontal dashed lines
between 5 and 15 km contain the majority of the events. The two off-fault high strain-drop clusters in
section II (IIA and IIB) are marked with dashed rectangles. The panel on the top shows the average log
10
Δε vs. horizontal coordinate for all events (solid line), as well as for events with depth between 5 and 15
km (dashed lines) for the constant (red) and 1D (black) velocity models. The right panel shows the
average log
10
Δε vs. depth in section II (dashed line) and in the whole aftershock zones (solid line) for the
constant (red) and 1D (black) velocity models.

2.4.2 Temporal Variations of Strain-drops
The availability in our data set of both foreshocks and aftershocks of the 1999 Düzce
mainshock, including many events near the hypocenter of the Düzce mainshock,
motivates us to search for possible informative temporal evolution of the derived strain-
drops. Figure 2.7 shows the obtained strain-drop distributions of earthquakes occurring
before and after the Düzce mainshock. Although the locations and rates of subsets of
the seismicity change dramatically (e.g. Peng & Ben-Zion 2005), the overall patterns of
52
strain-drops appear to be relatively stable before and after the occurrence of the Düzce
mainshock. To further quantify the temporal changes of strain-drops, we study the
evolution of log
10
strain-drops versus time inside each section (Figure 2.8). The median
values of strain-drops (red dashed lines) change after the Düzce mainshock in sections I
and II, although the differences are within the range of the standard deviation (grey
lines), and remain about the same in section III and IV.

Figure 2.7. Map views of strain-drops (color scale) of events before (a) and after (b) the Düzce mainshock
(red star). The dashed lines separate the sections with labels at the left-bottom corner.

One way of delineating the temporal evolution of the shapes of data points associated
with groups of events is to plot the 95% bounds of the distributions (e.g. Rolandone et
al. 2004). However, this method is suited for describing samples that have a large
53
number of events that obey the Gaussian distribution. In our case, the temporal
distribution of strain-drops does not satisfy either of these two conditions. Instead we
use the mean of the highest and lowest three strain-drop values in successive time
intervals of 8 days (black lines) to delineate the temporal evolution of the strain-drop
envelopes. The shapes remain stable when we use mean values based on two and four
points instead of three. With this simple method, we observe different temporal
characteristics in the strain-drop envelopes at the four different sections. The envelopes
for sections I and IV exhibit decays in the range of strain-drops with time from the İzmit
and Düzce mainshocks. The envelope in section II fluctuates within about constant
range, while the envelope in section III fluctuates at a constant level before the Düzce
mainshock and decays thereafter. Since the range of values is expected to increase with
increasing number of events (which lead to a more complete sampling of the available
distributions), the stability of the strain-drops in section II is somewhat surprising.
54

Figure 2.8. Temporal variation of strain-drops of earthquakes inside section I (a), II (b), III (c) and IV (d).
The red dashed lines show the median of log
10
strain-drops before and after the occurrence of the Düzce
earthquake. The solid blue lines in the middle and the gray lines above and below it are the median value
of log
10
Δε and corresponding ranges of standard deviations in 8 days intervals. The solid black lines are
used to delineate the temporal variation of the strain-drops in each section. The vertical dashed lines in
each subplot mark the time of the Düzce mainshock.

2.4.3 Station Spectral Terms
So far we focused on the source spectra and ignored the other terms of Eq. (2.1).
Since the seismic network was deployed close to the İzmit rupture zone (Figure 2.1), the
data provide important opportunity for studying station terms as a function of the
distance from the rupture zone. In particular, it is interesting to examine whether the
damaged low velocity fault zone layer that produces seismic trapped waves (Ben-Zion
et al. 2003) and non-linear wave propagation effects (Karabulut & Bouchon 2007; Wu
55
et al. 2009) in the region has a clear signature on the station spectra. Figure 9 presents a
summary of the station spectra obtained by stacking the station terms derived for each
event assuming a constant quality factor Q = 300. We divide the ten seismic stations to
three general categories based on their distances (top inset of Figure 2.9) from the
observed rupture traces: within the rupture zones (stations VO, LS and MO; red
symbols in Figure 2.9), near the rupture zones (stations FI and FP; blue symbols), and
away from rupture zone (stations BV, BU, CH, GE and WF; green symbols).

Figure 2.9. Stacked station spectra terms. The vertical dashed lines mark peak frequencies that are shown
in various stations. The inset top-right panel shows the distribution of stations (triangles) in relation to
observed rupture traces (black dashed line). The lower inset gives the integral of the site spectral energy
for each station. The employed colors for stations and spectra (legend at the top-left corner) are based on
distances from the rupture traces: red - within the rupture zone, blue - within 500 m from the rupture
zone, green - more than 1.0 km from the rupture zone.

The station spectral shapes are overall relatively smooth, reflecting the fact that each
56
shape is stacked from 7498 separated station terms. We observe three clear peaks in the
spectra of various stations around 6, 14 and 25 Hz. The peak frequency around 6 Hz
appears most clearly for stations within and near the rupture zones, and it generally
decays with distance from the rupture zone. This peak frequency is also observed in the
trapped waves study of Ben-Zion et al. (2003) using shear waveforms in the same
region, as well as the trapped S waves studies near the rupture zone of the 1992 Landers
earthquake (Peng et al. 2003) and the San Jacinto fault in California (Lewis et al. 2005).
These studies also show that the spectral amplitude around this peak is inversely
correlated with the distance from the rupture trace. The other two peaks around 14 and
25 Hz are seen clearly at a mixture of stations. Stations VO and FI that are within or
near the fault zone and off-fault station GE have a clear peak frequency near 14 Hz, and
their highest spectral amplitudes are associated with that peak. The 25 Hz peak
frequency is present at stations VO, FI, FP and LS which are within or near the rupture
zone. The large spectral peak at 14 Hz in the records of stations VO and FI may reflect
a damaged fault zone layer that is observed in the P wave analysis of this work and is
probably narrower than the trapping structure of the shear waves. The origin of the 14
Hz spectral peak at station GE, and the fact that fault zone station MO does not have
spectral peaks at 14 and 25 Hz remain at present unclear. Further observations and
analysis are needed to better explain these results.
The “site energy” contained in the different station spectral terms can be estimated
by df f fu
∫
2
)] ( 2 [ π , giving the integrated square of the station displacement spectra over
the examined frequency range. The inset in the low-right portion of Figure 2.9 gives the
obtained log
10
squared root values of the station spectral terms integrated over the
57
frequency range 2.35 - 30 Hz. The data of station VO that is within the highly damaged
rupture zone have the largest site energy, while the seismic records at the off-fault
stations BU, CH and BV have the smallest site energy. The order of the site energies of
the other stations does not show a clear relation to the distance from the fault and is
probably associated with various local site effects.

2.5 Discussion
A major general goal of earthquake seismology is to establish connections between
properties of fault zone structures, properties of earthquake sources on the faults and
generated seismic radiation at various stations. In this work we derived source, travel-
time and station spectra for groups of earthquakes from observed P-waves of 7498
events in the aftershock zones of the 1999 İzmit and Düzce earthquakes along the
Karadere-Düzce branch of the NAF. The analysis is based on the procedure of Shearer
et al. (2006) with small modifications and employed data-stacking multiple times. At
the iterative separation of the source, station and travel-time terms (section 2.3.1), the
observed source spectra of 201 nearest neighboring events, recorded by various stations
and shared within given epicentral distance intervals, are stacked. In section 2.3.2, the
separated source terms are stacked within given amplitude bins. At the strain-drop
fitting stage (section 2.3.3), the differences between binned separated source spectra and
theoretical source spectra are stacked to form a common EGF. Through stacking,
random noise in each individual spectrum is largely suppressed and the common
features of the source, travel-time and station spectra are highlighted. However, the
stacking may also suppress genuine internal variations of properties and it limits the
58
resolution of the results to common properties of the used groups of events.
We fit earthquake strain-drops from the separated source spectra with the ω
-2
source
model. The mean of the log
10
of strain-drop for the data is -3.5 (corresponding to a
mean stress drop value of 9.5 MPa). The spatial distributions of seismicity and derived
strain-drops can be used, along with the overall geometrical properties of the faults, to
divide the study area into 4 sections (Figures 2.5-2.7). The distributions of strain-drops
show temporal evolutions that probably reflect primarily the temporal changes in the
numbers of events (Figure 2.8). The ranges of the derived values are larger in the early
time of the aftershock sequences and decay gradually with time. In some sections there
are small step-like changes in the average strain-drop values across the time of the
Düzce mainshock, but the results are within the standard deviations and may not be
statistically significant.
Using a constant velocity model, Shearer et al. (2006) found that the median stress-
drops of earthquakes in southern California increases with depth from 0.6 MPa near the
surface to 2.2 MPa at 8 km. Our results from the İzmit-Düzce aftershock zones with a
constant velocity model show a similar pattern, with larger changes of values, in which
the mean log
10
strain-drops of all events increases from -4.2 (1.9 MPa) at 0-4 km depth
to -3.3 (15.0 MPa) at 10 km depth. We note that the strain-drops over the depth range 0-
4 km are associated primarily with events in region IIB (NW of station GE), but the
overall tendency of strain-drops to increase with depth between 4 and 10 km is a
general characteristic of the different regions. Fletcher & McGarr (2006) inferred from
analysis of slip models that stress drops increase from 10 MPa at 6 km to 43 MPa at
16.5 km for the Northridge earthquake, and from 7 MPa at 1.5 km to 20 MPa at 6.2 km
59
for the Landers earthquake. Allmann & Shearer (2007) showed that the depth-
dependency of stress drops in the Parkfield section of the San Andres Fault disappears
when a more accurate depth-dependent velocity model replaces the constant velocity
model.
To test if the depth variations of our derived strain-drops are controlled by the
velocity model, we recomputed the strain-drops with a 1-D velocity model for the study
area (black line in left panel of Figure 2.2.8). While the derived log
10
strain-drop values
change somewhat, the overall pattern of increasing strain-drops with depth remains. The
derived values of strain-drops are likely to change further if we use a 3D velocity
model. In particular, additional reductions of the shallow S-wave velocities would
increase the derived shallow strain-drops and reduce the overall change of strain-drops
with depth. However, the employed 1-D velocity model probably captures the overall
depth variations of velocities, and the obtained increasing strain-drops with depth likely
reflect, at least partially, the increasing confining pressure with depth and the reduction
of the stored elastic strain toward the surface and the brittle-ductile transition.
The properties of earthquakes along the relatively straight Karadere segment in
section II are different in many aspects from those of earthquakes in the other sections
(Figures 2.5 – 2.8). Seismotectonically, sections I and III are inside the extensional
Akyazi and Düzce basin, respectively. Both basins belong to a series of active pull-apart
structures from Düzce to the Sea of Marmara (e.g. Sengor et al 2005). The dominant
focal mechanisms in sections I and III are normal, and the dominant focal mechanisms
along the Karadere segment in section II are strike-slip faults (Bohnhoff et al. 2006).
The rupture zone of the İzmit earthquake changes its strike from nearly EW in the
60
Akyazi basin to N71°E along the Karadere segment. The rupture zone of the Düzce
earthquake changes back to nearly EW direction in the Düzce basin. We note that
abundant aftershocks tend to concentrate near stations CH and BV where the fault
rupture changes its strike, probably due to reactivations of secondary normal faults by
stress concentrations near abrupt changes in the geometry of the mainshock ruptures.
The low strain-drops inside section II probably reflect the geometrical simplicity (and
hence reduced stress concentrations) along that section. Compared with the clear
lineation of microseismicity inside section II, the earthquake distributions in sections I,
III and IV are considerably more distributed implying strong geometrical
heterogeneities. Additional evidence of fault heterogeneity could be derived from the
range of event sizes, which should become broader as the fault structures become more
heterogeneous (Ben-Zion 1996; Zöller et al. 2005). This is consistent with the fact that
the maximum magnitude is the smallest in section II among all sections.
The analysis of the stacked station spectra indicates strong variations of site effects
and three frequencies with clear peak amplitudes that are observed at various stations
(Figure 2.9). The amplitude of the station spectra are generally anti-correlated with the
distance to the rupture traces up to a frequency of about 6 Hz, where all station spectra
have either a peak or an inflexion point. Some of the stations have peak spectral
amplitudes around 14 and 25 Hz. The peaks at 6 and 14 Hz may be related to fault zone
trapped waves in the S and P portions of the seismograms, respectively. The integral of
the station spectra over the frequency range 2.35-30 Hz is largest for the fault zone
station VO that records clear trapped waves and temporal changes of properties (Ben-
Zion et al. 2003; Peng & Ben-Zion 2006; Wu et al. 2009) and smallest for the off-fault
61
stations BU, CH and BV. The derived patterns of source and site terms should be
substantiated by additional results based on S waveforms and other types of analysis.
This will be done in future work.

62
Chapter 2 References

Abercrombie, R., 1995. Earthquake source scaling relationships from −1 to 5 M
L
using
seismograms recorded at 2.5 km depth, J. Geophys. Res., 100, 24,015–24,036.

Allmann, B.P. & Shearer, P.M., 2007. Spatial and temporal stress drop variations in
small earthquakes near Parkfield, California, J. Geophys. Res., 112, B04305,
doi:10.1029/2006JB004395.

Aki, K. & Richards, P.G., 2002. Quantitative Seismology, University Science Books.

Andrews, D.J., 1986. Objective determination of source parameters and similarity of
earthquakes of different size, in Earthquake Source Mechanics (S. Das, J.
Boatwright, and C.H. Scholz, editors), American Geophysical Union Monograph, 37,
259–267.

Ben-Zion, Y., 1996. Stress, slip, and earthquakes in models of complex single-fault
systems incorporating brittle and creep deformations, J. Geophys. Res., 101, 5677-
5706.

Ben-Zion, Y., 2003. Appendix 2, Key Formulas in Earthquake Seismology, in
International Handbook of Earthquake and Engineering Seismology, eds. W. HK
Lee, H. Kanamori, P. C. Jennings, and C. Kisslinger, Part B, 1857-1875, Academic
Press.

Ben-Zion, Y., 2008. Collective Behavior of Earthquakes and Faults: Continuum-
Discrete Transitions, Evolutionary Changes and Corresponding Dynamic Regimes,
Rev. Geophysics, 46, RG4006, doi:10.1029/2008RG000260.

Ben-Zion, Y., Armbruster, J.G., Ozer, N., Micheal, A.J., Baris, S., Aktar, M., Peng, Z.,
Okaya, D. & Seeber, L., 2003. A shallow fault-zone structure illuminated by trapped
waves in the Karadere-Düzce branch of the North Anatolian Fault, western Turkey,
Geophys. J. Int., 152, 699–717.

Ben-Zion, Y., & Zhu, L., 2002. Potency-magnitude Scaling Relations for Southern
California Earthquakes with 1.0 < M_L < 7.0, Geophys. J. Int., 148, F1–F5.

Berckhemer, H., 1962. Die Ausdehnung der Bruchfläche im Erdbebenherd und ihr
Einfluss auf das seismische Wellenspektrum, Gerlands Beitr. Geophys., 71, 5–26.

Boatwright, J., Fletcher, J.B., & Fumal, T.E., 1991. A general inversion scheme for
source, site, and propagation characteristics using multiply recorded sets of
moderate-sized earthquakes, Bull. Seismol. Soc. Am., 81, 1754–1782.

63
Bohnhoff, M., Grosser, H. & Dresen, G., 2006. Strain partitioning and stress rotation at
the North Anatolian Fault after the 1999 İzmit Mw=7.4 earthquake, Geophys. J. Int.,
166, 373-385.

Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from
earthquakes, J. Geophys. Res., 75, 4997–5009.

Creager, K.C., 1997, Coral, Seismol. Res. Lett., 68, 269–271.

Fletcher J.B. & McGarr, A., 2006. Distribution of stress drop, stiffness, and fracture
energy over earthquake rupture zones, J. Geophys. Res., 111, B03312,
doi:10.1029/2004JB003396.

Hough, S.E., 1997. Empirical Green's function analysis: Taking the next step, J.
Geophys. Res., 102, 5369–5384.

Ide, S., & Beroza G.C., 2001. Does Apparent Stress Vary with Earthquake Size?,
Geophys. Res. Lett., 28(17), 3349–3352.

Karabulut, H. & Bouchon, M., 2007. Spatial variability and non-linearity of strong
ground motion near a fault Geophys. J. Int., 170, doi:10.1111/j.1365-
246X.2007.03406.x

Lewis, M.A., Peng Z., Ben-Zion Y., & Vernon F.L., 2005. Shallow Seismic Trapping
Structure in the San Jacinto fault zone near Anza, California , Geophys. J. Int., 162,
867.881, doi:10.1111/j.1365-246X.2005.02684.x, 2005.

Park, J., Lindberg, C.R. & Vernon, F.L., 1987. Multitaper spectral analysis of high
frequency seismograms, J. Geophys. Res., 92, 12,675–12,648.

Peng, Z., Ben-Zion Y., Michael A.J. & Zhu L., 2003. Quantitative analysis of seismic
trapped waves in the rupture zone of the 1992 Landers, California earthquake:
Evidence for a shallow trapping structure, Geophys. J. Int., 155, 1021–1041.

Peng, Z. & Ben-Zion, Y., 2004. Systematic analysis of crustal anisotropy along the
Karadere-Düzce branch of the north Anatolian fault, Geophys. J. Int., 159, 253–274.

Peng, Z. & Ben-Zion Y., 2005. Spatio-temporal variations of crustal anisotropy from
similar events in aftershocks of the 1999 M7.4 İzmit and M7.1 Düzce, Turkey,
earthquake sequences, Geophys. J. Int., 160, 1027–1043.

Peng, Z. & Ben-Zion, Y., 2006. Temporal changes of shallow seismic velocity around
the Karadere-Düzce branch of the north Anatolian fault and strong ground motion,
Pure and Appl. Geophys. 163, 567–600.
64

Prieto, G., Shearer, P.M., Vernon, F.L. & Kilb, D., 2004. Earthquake source scaling and
self–similarity estimation from stacking P and S spectra, J. Geophys. Res., 109,
B08310, doi:10.1029/2004JB003084.

Madariaga, R., 1976. Dynamics of an expanding circular fault, Bull. Seismol. Soc. Am.,
66, 639–666.

Mueller, C.S., 1985. Source pulse enhancement by deconvolution of an empirical
Green's function, Geophys. Res. Lett., 12, 33–36.

Rolandone, F., Bürgmann R. & Nadeau R.M., 2004. The evolution of the seismic-
aseismic transition during the earthquake cycle: Constraints from the time-
dependent depth distribution of aftershocks, Geophys. Res. Lett., 31, L23610,
doi:10.1029/2004GL021379.

Seeber, L., Armbruster, J.G., Ozer, N., Aktar, M., Baris, S., Okaya, D., Ben-Zion, Y. &
Field, N., 2000. The 1999 Earthquake Sequence Along the North Anotolian
Transform at the Juncture Between the Two Main Ruptures: Published in Special
Volume The 1999 İzmit and Düzce Earthquakes: Preliminary Results, Eds. Barka,
A., Kozaci, O., Akyuz, S., and Altunel, S., ITU, Istanbul.

Sengor, A.M.C., Tuysuz, O., Imren, C., Sakınc¸ M., Eyidogan, H., Gorur, N., Le
Pichon, X. & Rangin, C., 2005. The North Anatolian Fault: A New Look, Annu.
Rev. Earth Planet. Sci., 33, 37–112, doi: 10.1146/annurev.earth.32.101802.120415.

Shearer, P., Prieto, G. & Hauksson, E., 2006. Comprehensive analysis of earthquake
source spectra in southern California, J. Geophys. Res., 111, B06303,
doi:10.1029/2005JB003979.

Thomson, D.J., 1982. Spectrum estimation and harmonic analysis, Proceeding of the
IEEE., 70, 1055–1096.

Warren, L.M. & Shearer, P.M., 2002. Mapping lateral variations in upper mantle
attenuation structure by stacking P and PP spectra, J. Geophys. Res., 107(B12),
2342, doi: 10.1029/2001JB001195.

Wu, C., Peng Z. & Ben-Zion Y., 2009. Non-linearity and temporal changes of fault
zone site response associated with strong ground motion, Geophys. J. Int., 176,
265-278, doi: 10.1111/j.1365-246X.2008.04005.x.

Zöller, G., Holschneider, M. & Ben-Zion, Y., 2005. The role of heterogeneities as a
tuning parameter of earthquake dynamics, Pure Appl. Geophys., 162, 1027-1049.
65
Chapter 3

Estimation of corner frequency ratios from P and S
waves generated by aftershocks around the Karadere-
Düzce branch of the North Anatolian Fault

SUMMARY

We propose a method to systematically estimate the corner frequency ratios of P and
S waves, and apply the method to obtain corner frequency ratios generated by ~9,000
aftershocks of the 1999 İzmit and Düzce earthquakes recorded at 10 seismic stations
along the Karadere-Düzce branch of the North Anatolian Fault. The analysis is
associated with separation of source, travel-time and station spectra terms and stacking
results at several stages to enhance the signal-to-noise ratio. We analyze source spectra
of P waves selected from a 1.28-second window around the P arrival, and examine
results associated with S wave signals selected using the following 5 schemes: (A) A
single 1.28-second time window for S wave. (B) A single 10-second time window for S
wave. (C1-C3) Stacked S wave spectra in N (N = 15) 1.28 second consecutive S wave
windows, where the spectral amplitudes are stacked according to travel-time
normalization (C1), amplitude normalization (C2), and without any normalization (C3).
The corner frequency ratios and relative strain-drops are obtained by fitting iteratively
the separated P and S source spectra of 201 nearest neighboring events in different
66
amplitude bins to the omega-2 source spectral model. The mode values of the obtained
P/S corner frequency ratios are 1.52, 1.57, 1.89, 1.47 and 1.56 for each of the five
schemes, respectively. The obtained median values for log
10
relative strain-drops for all
schemes are -4.35, which corresponds to a stress drop of 1.34 MPa assuming a nominal
rigidity of 30 GPa. Schemes A and B have larger fitting errors compared with scheme
C1-C3. Scheme C1 over-estimates the corner frequency ratio compared with the other
four schemes, and the obtained corner frequency ratios given by scheme C1 and C3
have higher standard deviations compared with the other schemes. The variations of the
P/S corner frequency ratios obtained for different groups of events may be related to
changes in average azimuths between the event locations and stations, variable focal
mechanisms spatial changes of fault heterogeneities, and rupture directivity effects.

3.1 Introduction
Earthquake source properties and associated rupture process could be inferred from P
and S waves recorded on the seismogram. Numerous investigations have been
concerned with acquisition earthquake source properties from seismic data and
theoretical simulation. One particular source parameter is the corner frequency, which is
defined as the frequency where the high and low frequency trends intersect on the shape
of displacement seismic source spectrum (e.g. Molnar et al., 1973). Physically, the
corner frequency is uniquely related to the fault area (Brune, 1970). Due to the different
energy release modes, the corner frequencies for P and S waves are different, and their
ratio is referred as corner frequency ratio, or corner frequency shift (Hanks, 1981).
Different theoretical models were constructed with relating corner frequency to
67
different properties of the earthquake source and predict different corner frequency
ratios (Haskell, 1964; Brune, 1970; Savage, 1972; Dahlen, 1974; Madariaga, 1976; Sato
& Hirasawa, 1973). The predicted P/S corner frequency ratios range from less than 1.0
(Savage, 1972; Dahlen, 1974) to more than 1.0 (Sato & Hirasawa, 1973; Madariaga,
1977). For the models that predict average P/S corner frequency ratios above 1.0, the
predicted average P/S corner frequency ratios range from 1.3 (Sato & Hirasawa, 1973)
to 1.5 (Madariaga, 1977). The model given by Madariaga (1977) predicted that P/S
corner frequency ratios are less than 1.0 with azimuth angle below 30 degree, while the
model given by Sato and Hirasawa (1973) predicted that corner frequency ratio is
always above 1.0. Numerous studies had been done to estimate the P/S corner
frequency ratios from individual to hundreds of natural earthquakes (e.g. Furuya, 1969;
Wyss and Hanks, 1972; Trifunac, 1972; Molnar et al, 1973; Abercrombie, 1995; Prieto
et al, 2004), but some debates were concerned on the significance of results with few
observed data and uncertainty associated with these methods due to the existence of
attenuation and site effect (Hanks, 1981, 1982; Burdick, 1982, Langston, 1982).
Therefore, systematically separating both P and S source spectra from observed seismic
waveforms, and calculating the corner frequency ratios of earthquakes in a large data set
would provide important and reliable constraint on various model predictions and
benefit our understanding on the features of seismic source physics.
By using waveform data within a single [-0.28 - 1.00] second P wave window from a
group of neighboring events (e.g. N = 200), Shearer et al (2006) proposed a method
which includes a series of procedures to separate seismic source spectra, and thereby fit
for stress drop to the Brune type ω
-2
source model. Through stacking similar terms
68
iteratively, Shearer’s method can efficiently decrease random noise and isolate the
source spectra. Such a method together with its prototypes had been used to obtain
seismic source parameters and study source scaling relationship in different regions at
different scales (Warren and Shearer, 2002; Prieto et al, 2004; Shearer et al, 2006;
Allmann & Shearer 2007, 2009; Yang et al, 2009).
With a slower velocity, the recorded signal of S wave is inevitable mixed with the P
wave coda on seismogram, especially at the part of high frequency (Prieto et al, 2006).
It might not be able to distinguish the source spectra associated with S wave from the
preexisted P wave source spectra if we take the same 1.28-second time window for S
wave as it is for the P wave. Since there is no preexisting procedure to process S wave,
we try to systematically analyze the S wave according to several different schemes in
this study. The first scheme (A) is to simply consider P and S wave in the same [-0.28 -
1.00] second time window related to each arrival respectively. Aki (1982) identified that
the S wave coda is essential the scattering S wave. For the second scheme (B), we
consider S wave in a [-0.28 - 9.72] second time window related to the S arrival. Another
way is to take multiple time windows that contain the target signal, and stack the spectra
among them to increase the significance of target signal (e.g. Imanishi and Ellsworth,
2006). Therefore, for the third scheme (C), we take multiple time windows (N = 15)
following the S arrival, and stack the calculated spectra among these windows
according to three different sub-schemes: Normalized by travel-time difference (C1),
amplitude ratio (C2), and without normalization (C3).

69

3.2 Data and Preprocessing

Figure 3.1. Hypocentral distribution of ~26000 earthquakes (green points) along the Karadere-Düzce
branch of the North Anatolian Fault. The corner frequency ratios of 8,785 events (red points) are
estimated in this study. The gray lines denote observed rupture traces of the İzmit and Düzce mainshocks
and the black triangles are the ten seismic stations. Focal mechanisms of both the İzmit and Düzce
mainshocks point to their epicenters. The inset map shows the location of the study area (black rectangle)
in a large-scale regional map. Twelve events (numbers in circles) are selected to illustrate the spectra
fitting.

A temporary PASSCAL network of 10 short-period seismic stations was deployed
along and around the Karadere-Düzce branch of NAF a week after the August 17, 1999,
Mw7.4, İzmit, Turkey, earthquake (Figure 3.1). All stations had REFTEK recorders and
three-component L22 velocity sensors with a sampling frequency of 100 Hz (Seeber et
al. 2000; Ben-Zion et al. 2003). This network operated for 6 months and recorded
approximately 26,000 events, including the November 12, 1999, Mw7.1, Düzce
earthquake. In this study, we use 8,785 events with magnitudes in the range between 0.0
70
and 3.8. The horizontal location errors are generally less than 1 km around the center of
the network and 1-2 km near the margins. The vertical location errors are somewhat
larger. Additional details on the seismic experiment and data set are given by Seeber et
al. (2000) and Ben-Zion et al. (2003).
We analyze both P and S waves recorded on all three components of the short-period
instruments. The preprocessing steps are as follows: We apply a simple algorithm
(Yang and Ben-Zion, 2009) to detect clipped waveforms, and remove the associated
clipped events. We remove the mean and trend from each seismogram, correct for the
instrumental response, and convert the waveforms from velocity to displacement
records. These analysis steps are done using Matlab functions in the CORAL package
for seismic waveform analysis (Creager 1997). The displacement spectrum in each
waveform window is computed using a multitaper algorithm (Thomson 1982; Park et
al. 1987). We take the P wave window to be 0.28 s before and 1.0 s after the picked P
arrival, and take a background noise window to be with the same length before the
picked P arrival. We require the obtained P wave to noise spectra ratio on the vertical
component in the frequency band [4.0 - 30.0] Hz to be larger than 3.0. We also require
for further analysis the minimum number of qualified observed waveforms for an event
to be larger than 3. We calculate the total spectra for signals from all components
through vector summation. After data preprocessing, a total of 8,785 aftershocks with
47,172 P/S wave spectra fulfill the required criteria.
For the calculation of S wave spectra, we consider three different methods, namely,
scheme A, B and C. In scheme A, we take waveform in 0.28 s before and 1.0 s after the
picked S arrival time window. In scheme B, we take waveform in 0.28 s before and 9.72
71
s after the picked S arrival time window. In scheme C, we take S waveform in
successive (N =15) multiple 1.28-second windows with the first window in 0.28 s
before and 1.0 s after the picked S arrival (Figure 2). Examples of waveforms and
corresponding signal and noise time window are shown in Figure 3.2. Additionally, for
scheme C, we take spectrum in a 1.28 s window ahead of P arrival to be the noise
spectrum, and we keep the observed spectra in those S windows with SNR larger than 3
on each component. For the satisfied spectra among the 15 multiple S wave windows,
we consider three different ways for normalization: C1) Normalize the S wave spectra
amplitudes according the travel-time difference between each S wave window and the
first S wave window. C2) Normalize the S wave spectra amplitudes according the
absolute maximum amplitude ratio between each S wave window and the first S wave
window. C3) Without normalization. Then after, we stack the processed observed
spectra on each component. All five schemes are summarized in Table 3.1.
72

Figure 3.2. The selection of noise, P wave and S wave windows for schemes A (a), scheme B (b) and
scheme C (c). In each panel, green color for noise, red color for P wave and blue color for S wave. In
panel (c), 15 consecutive 1.28-second windows are selected for S wave.
73

Scheme
P wave
(Relative to P
arrival)
S wave
(Relative to S
arrival)
Normalization
A [-0.28 – 1.00] s No
B [-0.28 – 9.72] s No
C1 Travel-time
C2 Amplitude
C3
[-0.28 – 1.00] s
15 consecutive
[-0.28 – 1.00] s
No
Table 3.1. Five schemes for P and S wave spectra calculations

3.3 Joint Fitting P and S wave for Corner Frequency Ratio
To each event and its 200 neighbors in the 8,785 events data set, we generally
follow the procedures of Shearer et al (2006) to separate seismic source, travel-time and
station terms iteratively, and the detailed analysis steps were documented in Yang et al
(2009). To the separated P and S source spectra, we stack them over 0.2 amplitude bins
respectively. Due to the Gutenberg-Richter distribution and the limited detecting
capability of network, the numbers of spectra are not distributed evenly among all
amplitude bins. Thus, to decrease the variation of individual spectra for source spectra
fitting, we require the minimum number of spectra in each amplitude bin to be 10.
Given the binned separated P and S source spectra for each scheme, we fit them
with the Brune-type ω
-2
source model (Brune, 1970):
Ω( f)=
Ω
0
1+( f / f
c
)
2
, (3.1)
74
where Ω
0
is the displacement amplitude at zero frequency approximation, and f
c
is the
corner frequency. The corner frequency is related with stress drop (Δσ ) and moment
(

M
0
) or strain-drop (Δε ) and potency (

P
0
), shear wave velocity (β ) and coefficient
(m) given by different models with the following formula:
f
c
=mβ
Δε
P
0
⎛
⎝
⎜
⎞
⎠
⎟
1/3
=mβ
Δσ
M
0
⎛
⎝
⎜
⎞
⎠
⎟
1/3
(3.2)
The P/S corner frequency ratio is the ratio of coefficients between P and S wave spectra
through fitting with Eq 3.1 and 3.2. For example, Madariaga (1976) obtained the
average coefficient is 0.42 for P wave and 0.28 for S waves theoretically, assuming a
constant rupture velocity of 0.9 times the shear wave velocity.
In our analysis, we fit for P and S source spectra in [4 - 30] Hz frequency band, and
the employed spectra decay equation is:
Ω( f )=
Ω
f
0
[1+ ( f
0
/ f
c
)
2
]
1+ ( f / f
c
)
2
. (3.3)
Where

f
0
= 4Hz, and

Ω
f
0
is the source displacement amplitude at

f
0
.
We assume a constant rupture velocity with 0.9 of the shear wave velocity, and use
the following two equations to represent corner frequencies for P and S waves:
f
cp
= 0.42β
Δε'
Ω
f
0 p
⎛
⎝
⎜
⎞
⎠
⎟
1/3
, (3.4)
and
f
cs
=
0.42β
R
fc
Δε'
Ω
f
0s
⎛
⎝
⎜
⎞
⎠
⎟
1/3
. (3.5)
75
Where

Ω
f
0 p
and

Ω
f
0 S
are the displacement amplitudes at frequency

f
0
for P and S
waves respectively.

R
fc
is the P/S corner frequency ratio. is the relative strain-drop,
and absolute strain-drop could be obtained by calibrating it with catalog magnitude.
Through using Eq (3.4) and (3.5), we can avoid including error associated magnitude
calibration.
With Eq (3.3) - (3.5), we jointly fit P and S source spectra to estimate relative strain-
drop and corner frequency ratio with the following grid-search steps:
1) Give a range of possible log
10
strain-drops values (e.g. 100 values linearly
distributed in [-6.0 - -2.0]). For each strain-drop, do the following:
a) Calculate the theoretical P wave source spectra.
b) Give a range of possible corner frequency ratio values (e.g. 50 values linearly
distributed in [0.1 - 6.0]), and calculate the theoretical S wave source spectra.
c) Find the differences between binned source spectra and theoretical spectra, stack
them to be the common empirical Green’s functions (EGFs), and remove the
common EGFs for both P and S source spectra respectively. Then calculate the
average fitting RMS among all spectra in the [4 – 30] Hz frequency band as
fitting error.
2). Assign the strain-drop and corner frequency ratio that are associated with the
minimum fitting error to be the best fitted strain-drop and corner frequency ratio for this
event.

76
3.4 Application and Results
For each of the 8,785 events, we use the above procedures to fit for P and S source
spectra and obtain the best-fitting corner frequency ratio and relative strain-drop for
each event to each of the five schemes. We use a 1D velocity model for
Akyazi/Karadere-Düzce region: Vp (km/s) = [2.90 4.70 5.30 6.06 6.88 8.06], depth
(km) = [0.0 1.0 7.0 11.0 15.0 35.0] (Bulut et al, 2007). We set the range of corner
frequency ratio to be in [0.1 - 6.0], and the range of log
10
relative strain-drop to be in [-
6.0 - -2.0].
The statistical features of obtained corner frequency ratios, relative strain-drops and
fitting errors for the five schemes are shown in Figure 3.3. The distributions of corner
frequency ratios for all schemes are of the skewed distribution (panels a, d, g, i and m in
Figure 3.3), and the associated mode values are 1.52, 1.57, 1.89, 1.47 and 1.56 for each
of the five schemes, respectively. There are less than 3% of all events that could be not
be fitted well with obtained corner frequency ratios reaching the upper boundary of the
given parameter range. The obtained relative strain-drops in log10 scale share the same
distribution for all schemes and have the same median values (-4.35). The distributions
of log
10
fitting errors follow the Gaussian distribution for all schemes, and schemes with
stacked observed spectra among multiple S wave windows (C1-3) have smaller fitting
errors compared to schemes A and B that have observed spectra in a single S wave
window.
The mode values of corner frequency ratios for scheme A, B, C2 and C3 are very
close, and these values are also close to the theoretical average corner frequency ratios
(=1.5) given by Madariaga (1976). The percentages of corner frequency ratios with
77
values less than 1.0 in all schemes are smaller than it is in Madariaga (1976), which
may be due to that most events in our data set located in a limited azimuth range to the
network. Scheme C1 and C3 have larger standard deviations of corner frequency ratios
compared with scheme A, B and C2.

Figure 3.3. The statistical features of best-fitting corner frequency ratios, log
10
relative strain-drops and
log
10
fitting error for the five schemes (the name of scheme is at the left side). In side each panel, the
mode/median value for each distribution is marked by red line and labeled at up-right corner.

To illustrate the P and S spectra fitting, we randomly select twelve events with spatial
distributed equally in the aftershock zone (circles with labels in Figure 3.1) among the
8,785 events, and plot the spectra fitting for both P and S wave to each scheme (Figure
3.4a-e).

78

Figure 3.4a. The best-fitting P and S spectra to twelve events and their 200 neighbors for scheme A. Blue
curves are EGF removed source spectra, and the red curves are theoretical source spectra. The fittings for
P and S wave spectra are given in two successive panels with “P-“ and “S-“ marked at the up-right
corners, and the event is labeled at the up-left corner in each P source spectra fitting panel to match with
labels in Figure 3.1. The best-fitting corner frequency ratio is marked at the left-bottom corner in each
panel for the S source spectra fitting.

79

Figure 3.4b. The best-fitting P and S spectra to twelve events and their 200 neighbors for scheme B.
Labels and legends are the same as in Figure 3.4a.
80

Figure 3.4c. The best-fitting P and S spectra to twelve events and their 200 neighbors for scheme C1.
Labels and legends are the same as in Figure 3.4a.
81

Figure 3.4d. The best-fitting P and S spectra to twelve events and their 200 neighbors for scheme C2.
Labels and legends are the same as in Figure 3.4a.

82

Figure 3.4e. The best-fitting P and S spectra to twelve events and their 200 neighbors for scheme C3.
Labels and legends are the same as in Figure 3.4a.

83

Figure 3.5. Best-fitting corner frequency ratios for scheme A (a), B (b), C1 (c), C2 (d) and C3 (e). Corner
frequency ratios are colored with scale at the right bottom. Stations are marked by triangles. For the effect
of visualization, the up-limit of colorbar is set to be 3.0.

84
Figure 3.5 shows the spatial distribution of corner frequency ratios of the 8,785
events for the five schemes. The corner frequency ratios distribution patterns of for all
five schemes generally match with each other. Due to the using of longer S wave
window and the stacking of S wave spectra over 15 consecutive 1.28-second windows,
more events in scheme B, C1 and C3 have corner frequency ratios with high values (Rcf
≥ 3), and distributed in broader area. Scheme C2 is less affected by the using of multiple
windows and has a corner frequency ratios distribution pattern most similar to it is of
scheme A among all schemes that use S coda waves. There are two common areas that
events with corner frequency ratio less than 1.0 are mostly concentrated in among all
schemes. One minor common area is at the east side of the aftershock region (Longitude
< 30.5). One major common area locates at the intersection between the strike of the
karadere segment and the strike of the Düzce aftershock zone. It is interesting to note
that beyond the intersection along the strike of the Düzce aftershock zone to the NW,
events with large corner frequency ratios (Rcf>=3) concentrate.

3.5 Discussion

We discuss a procedure to estimate P/S corner frequency ratio and relative strain-
drop. The method could be used to estimate seismic strain-drops and corner frequency
ratios systematically in a large data set.
We first separate P and S source spectra by generally using the method of Shearer et
al (2006). For the selection of waveform window for S wave, we test with five different
schemes. To the separated 201 source spectra, we require the minimum number of
85
spectra in each amplitude bin to be 10, and as a result, the stacked source spectra get
smoother and the associated fitting error decreases.
The strain-drop is a source parameter and keeps the same for P and S wave on
seismogram. By jointly fitting P and S wave spectra, we can estimate stable relative
strain-drop and solve for corner frequency ratio systematically. By using grid-search
steps, we obtain a pair of best-fitting relative strain-drop and corner frequency ratio that
minimize the source spectra fitting error. There exist some outliers (corner frequency
ratio equaling the up-limit of a given range) with a small percentage (1~3%).
Instead of calibrating seismic potencies with catalog magnitudes for the source model
fitting, we use the relative spectra amplitudes at low frequency directly to represent the
seismic potencies. With such an improvement, we can decrease the fitting error
dramatically, and avoid the inconsistence between Eq (3.1) and (3.2) in practical.
Scheme A takes the same single 1.28-second window around the phase arrivals for P
and S waves. Since the early part of S wave is under the shadow of P wave coda, we
were expected that the mode value of the obtained corner frequency ratios to be close to
1.0. However, Scheme A has a mode value to be 1.52. This value is close to the mode
values given by other schemes employing S coda wave, and is also close to the average
corner frequency ratio predicted by Madariaga (1976). Scheme B, C1, C2 and C3 use
longer S wave time windows with different approaches, and all give similar mode
values of corner frequency ratios except that scheme C1 has a larger mode value.
Among all schemes, scheme C1, C2 and C3 give smaller fitting errors, due to the
stacking of observed spectra among multiple windows. Scheme C2 has a smaller
standard deviation of corner frequency ratios compared with scheme C1 and C3.
86
The corner frequency ratios obtained from all schemes have values less than 1.0 with
different percentages. Those corner frequency ratios theoretically given by Madariaga
(1976) also have values less than one for events with azimuth less than 30 degree. The
spatial concentration of large and small values of corner frequency ratios in neighboring
areas with different azimuths could be found in Figure 3.5. However, due to the spatial
limitations of network layout and the concentration of events along the rupture zones,
we could not deduce a relationship between corner frequency ratio and azimuth.
Besides the changes in average azimuths between the event locations and stations, the
variations of corner frequency ratios might also be related to variable focal mechanisms,
spatial changes of fault heterogeneities, and rupture directivity effects.

87
Chapter 3 References

Abercrombie, R.E., 1995. Earthquake source scaling relationships from -1 to 5 ML
using seismograms recorded at 2.5-km depth, J. Geophys. Res., 100, 24015-24036.

Allmann, B.P. & Shearer, P.M., 2007. Spatial and temporal stress drop variations in
small earthquakes near Parkfield, California, J. Geophys. Res., 112, B04305,
doi:10.1029/2006JB004395.

Allmann, B.P. & Shearer, P.M., 2009. Global variations of stress drop for moderate to
large earthquakes, J. Geophys. Res., 114. B01310, doi:10.1029/2008JB005821.

Aki, K., 1982. Scattering and Attenuation, Bull. Seism. Soc. Am., 72, 319-330.

Ben-Zion, Y., Armbruster, J.G., Ozer, N., Micheal, A.J., Baris, S., Aktar, M., Peng, Z.,
Okaya, D. & Seeber, L., 2003. A shallow fault-zone structure illuminated by trapped
waves in the Karadere-Düzce branch of the North Anatolian Fault, western Turkey,
Geophys. J. Int., 152, 699–717.

Bohnhoff, M., Grosser, H. & Dresen, G., 2006. Strain partitioning and stress rotation at
the North Anatolian Fault after the 1999 İzmit Mw=7.4 earthquake, Geophys. J. Int.,
166, 373-385.

Boatwright, J. & Fletcher J.B., 1984. The partition of radiated energy between P and S
waves, Bull. Seism. Soc. Am., 74, 361-376.

Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from
earthquakes, J. Geophys. Res., 75, 4997–5009.

Bulut, F., Bohnhoff, M., Aktar, M. & Dresen, G., 2007. Characterization of aftershock-
fault plane orientations of the 1999 Izmit (Turkey) earthquake using high-resolution
aftershock locations, Geophys. Res. Lett., 34, L20306, doi:10.1029/2007GL031154.

Burdick, L.J., 1982. Comments on "the corner frequency shift, earthquake source
models, and Q," by T. C. Hanks, Bull. Seism. Soc. Am., 72, 1419-1426.

Castro, R.R., Anderson, J.G. & Brune, J.N., 1991. Origin of high P/S spectral ratios
from the Guerrero accelerogram array, Bull. Seism. Soc. Am., 81, 2268-2288.

Creager, K.C., 1997. Coral, Seismol. Res. Lett., 68, 269–271.

Dahlen, F.A., 1974. On the ratio of P-wave to S-wave corner frequencies for shallow
earthquake sources, Bull. Seism. Soc. Am., 64, 1159-1180.

Furuya, I., 1969. Predominant period and magnitude, J. Phys. Earth 17, 119-126.
88

Hanks, T.C., 1981. The corner frequency shift, earthquake source models, and Q, Bull.
Seism. Soc. Am., 71, 597-612.

Hanks, T.C., 1982. Reply to "comments on 'the corner frequency shift, earthquake
source models, and Q'", Bull. Seism. Soc. Am., 72, 1433-1444.

Hanks, T.C., and Wyss M., 1972. The use of body-wave spectra in the determination of
seismic source parameters, Bull. Seism. Soc. Am., 62, 561-589.

Imanishi, K. & Ellsworth, W.L., 2006. Source scaling relationships of microearthquakes
at Parkfield, CA, determined using the SAFOD pilot hole array. In Abercrombie,
R.E., McGarr, A., Di Toro, G., and Kanamori, H. (Eds.), Earthquakes: Radiated
Energy and the Physics of Faulting, Geophys. Monogr. 170, Washington D.C.
(American. Geophysical Union), 81–80.

Langston, C.A., 1982. Comments on "the corner frequency shift, earthquake source
models, and Q," by T. C. Hanks, Bull. Seism. Soc. Am., 72, 1427-1432.

Park, J., Lindberg, C.R. & Vernon, F.L., 1987. Multitaper spectral analysis of high
frequency seismograms, J. Geophys. Res., 92, 12,675–12,648.

Peng, Z. & Ben-Zion, Y., 2004. Systematic analysis of crustal anisotropy along the
Karadere-Düzce branch of the north Anatolian fault, Geophys. J. Int., 159, 253–274.

Prieto, G., Shearer, P.M., Vernon, F.L. & Kilb, D., 2004. Earthquake source scaling and
self–similarity estimation from stacking P and S spectra, J. Geophys. Res., 109,
B08310, doi:10.1029/2004JB003084.

Prieto, G.A., Thomson, D.J., Vernon, F.L. & Shearer P.M., 2006. How much does P-
wave coda bias S-wave spectral estimates?, Seismol. Res. Lett., 77, 167.

Madariaga, R., 1976. Dynamics of an expanding circular fault, Bull. Seismol. Soc. Am.,
66, 639–666.

Mayeda, K., Hofstetter, A., O'Boyle, J.L. & Walter, W.R., 2003. Stable and
Transportable Regional Magnitudes Based on Coda-Derived Moment-Rate Spectra,
Bull. Seism. Soc. Am., 93, 224-239.

Molnar, P., Tucker, B.E. & Brune, J.N., 1973. Corner frequencies of P and S waves and
models of earthquake sources, Bull. Seism. Soc. Am., 63, 2091-2104.

Savage, J.C., 1972. Relation of Corner Frequency to Fault Dimensions, J. Geophys.
Res., 77, 3788–3795.

89
Sato, T. & Hirasawa, T., 1973. Body wave spectra from propagating shear cracks, J.
Phys. Earth, 21, 415-431.

Seeber, L., Armbruster, J.G., Ozer, N., Aktar, M., Baris, S., Okaya, D., Ben-Zion, Y. &
Field, N., 2000. The 1999 Earthquake Sequence Along the North Anotolian
Transform at the Juncture Between the Two Main Ruptures: Published in Special
Volume The 1999 İzmit and Düzce Earthquakes: Preliminary Results, Eds. Barka,
A., Kozaci, O., Akyuz, S., and Altunel, S., ITU, Istanbul.

Shearer, P., Prieto, G. & Hauksson, E., 2006. Comprehensive analysis of earthquake
source spectra in southern California, J. Geophys. Res., 111, B06303,
doi:10.1029/2005JB003979.

Thomson, D.J., 1982. Spectrum estimation and harmonic analysis, Proceeding of the
IEEE., 70, 1055–1096.

Trifunac, M.D., 1972. Stress estimates for the San Fernando, California earthquake of
February 9, 1971: Main event and thirteen aftershocks, Bull. Seism. Soc. Am., 62,
721-750.

Warren, L.M. & Shearer, P.M., 2002. Mapping lateral variations in upper mantle
attenuation structure by stacking P and PP spectra, J. Geophys. Res., 107(B12),
2342, doi: 10.1029/2001JB001195.

Wyss, M. & Hanks, T.C., 1972. The source parameters of the San Fernando earthquake
inferred from teleseismic body waves, Bull. Seism. Soc. Am., 62, 591-602.

Yang, W. & Ben-Zion, Y., 2009. An algorithm for detecting clipped waveforms and
suggested correction procedures. Submit to Seismol. Res. Lett.,

Yang, W., Peng, Z. & Ben-Zion, Y., 2009. Variations of strain drops in aftershocks of
the 1999 Izmit and Düzce earthquakes along the Karadere-Düzce branch of the
North Anatolian fault, Geophys. J. Int., 177, 235-246.

90
Chapter 4

An algorithm for detecting clipped waveforms and
suggested correction procedures

SUMMARY
We discuss a detection algorithm and two possible methods for correcting clipped
waveforms. The procedures are applied to seismograms of ~26,000 earthquakes,
recorded with sampling frequency of 100 Hz by a 10-station PASSCAL array, where
values of clipped points are set to zero. We use a multi-step detection algorithm that is
based on the observed ranges of amplitudes in the recorded seismograms, expected
frequency-size statistics of amplitudes associated with the Gutenberg-Richter
distribution and properties of neighboring points in seismograms. As expected, the
horizontal components have more detected clipped waveforms than the vertical one, and
the intensity of waveform clipping increases with proximity to the fault and amplitude
of site effects. We compare a linear interpolation method with corrections of clipped
points using scaled-up versions of unclipped similar waveforms. With the sampling rate
and other properties of our data, the results may be summarized as follows. In cases
with 1-2 consecutive clipped points, the interpolation method generally produces
corrections with smaller errors. In cases with 3-5 consecutive clipped points, the similar
waveform method performs better if similar waveform with cross-correlation
91
coefficient (CCC) values larger than, respectively, 0.96, 0.91 and 0.88 are available. In
cases with 6 or more consecutive clipped samples, the similar waveform performs
generally better. However, similar waveforms with CCC values larger than 0.85 are
needed to produce corrected waveforms that may lead to error in earthquake magnitude
less than about 0.3.

4.1 Introduction
Waveform clipping due to instrumental limitation is a common problem in seismic
data recorded by local and regional networks. With high-amplitude data missing,
clipped waveforms, which usually have short epicentral distances and high signal-to-
noise ratios, have to be removed from the estimation of magnitudes and other
earthquake properties. An efficient automatic algorithm for detecting clipped
waveforms can help to eliminate artifacts associated with analysis of clipped
seismograms. In addition, a reliable algorithm for correcting some clipped waveforms
can increase the availability of data closer to the source location and data for larger
magnitude events.
The severity of waveform clipping at a given station is expected to increase with
earthquake size. Since the earthquake frequency generally decreases with size following
the Gutenberg-Richter relation, a large number of clipped waveforms are likely to have
only a small number of clipped points. In such cases, it may be possible to provide
reasonably accurate corrections for the clipped waveforms.
Interpolation methods are often used to reconstruct missing data. By using a spline
interpolation method, Karabulut and Bouchon (2007) corrected the peak values of
92
clipped waveform recorded by strong motion accelerometers during the 1999 Mw7.1
Düzce earthquake in Turkey. Subsets of earthquakes recorded by local or regional
network often have very similar waveforms with high cross-correlation coefficient
(CCC) values (e.g. Aster and Scott, 1993; Nadeau et al., 1994). Another possible way of
correcting clipped waveforms may be by scaling up portions of unclipped similar
waveforms with sufficiently high CCC values.
In the following, we analyze statistical features of clipped waveforms in a large data
set of 26,080 aftershocks recorded by a ten-station near-fault seismic network in 6
months following the Mw7.4, İzmit mainshock (Figure 4.1). Clipped waveforms can
generally be classified, based on the response of instruments to amplitude saturation,
into two types: Flat-Top (FT) where the value of clipped points are set to be the upper
limit value, and Back-to-Zero (BZ) where the value of clipped points are set to zero. In
our data set, the clipped waveforms are of the BZ type.
We propose an algorithm for automatic detection of clipped waveforms and discuss
the possibility of correcting clipped seismograms. About 3.6% of the recorded
waveforms in the data set are found to be clipped. We consider waveform corrections
based on a linear interpolation with the Kriging method and using unclipped similar
waveforms of smaller events. A comparison of the performance of these two methods,
using artificially clipped waveforms within a 10-second time window and with a 100
Hz sampling frequency, indicates that in cases of 1-2 consecutive clipped samples the
Kriging method performs better. In cases of 3-5 consecutive clipped samples, the
Kriging method performs better unless similar waveforms with CCC values larger than,
respective, 0.96, 0.91 and 0.88 are available. In case with 6 or more consecutive clipped
93
samples, the Kriging method has large correction error while the similar waveform
method general performs better. However, corrections using similar waveforms with
CCC<0.85 are likely to produce overly large errors for accurate magnitude
determinations.

4.2 Data
One week after the occurrence of the 1999 Mw7.4 İzmit earthquakes, a PASSCAL
seismic network with ten short-period seismic stations was deployed along the
Karadere-Düzce branch of the North Anatolian Fault (NAF). All stations had REFTEK
recorders and three component L22 velocity sensors with a sampling frequency of 100
Hz (Seeber et al. 2000; Ben-Zion et al. 2003). This network operated for 6 months and
recorded 26,080 events with magnitudes between -0.1 and 7.2 and hypocentral distances
less than 100 km (Figure 4.1).

94

Figure 4.1. Epicentral distribution of ~26,000 earthquakes along the Karadere-Düzce branch of the North
Anatolian Faults. The earthquake symbols are colored by depth (bottom left legend) and scaled by
magnitude (bottom right legend). Stations are marked by triangles with names underneath. The black and
gray curves are surface ruptures associated with the İzmit and Düzce mainshocks, respectively.

Waveforms recorded by the various instruments and components were clipped to
different degrees. Figure 4.2 shows an example set of three component clipped
waveforms recorded by station VO located within the Izmit rupture zone and generated
by an earthquake with M
L
=3.3.
95

Figure 4.2. Three components clipped waveforms recorded by the fault zone station VO. The amplitudes
of clipped samples were set to be zero by the instrument.

4.3 Features of Clipped Waveforms
Through visually inspecting clipped waveforms in the data set, we can summarize
the features of clipped waveforms as follows: 1) Clipped waveforms contain points with
amplitude close to the upper range of the recorded data. 2) Clipped points are in
waveform portions with amplitude close to the peak value in the seismogram. 3)
Clipped points are located within waveform sections having large amplitude fluctuation.
4) All clipped data points have zero amplitude values. 5) Data on both sides of a clipped
point have the same sign or one neighboring point is also zero. We define “the upper
96
range of the recorded data” in feature (1) to be the Observed Range of the relevant
instrument.
Features (1) and (2) are intuitive, and feature (3) limits further the location of
clipped sample points inside the waveform. Features (4) and (5) represent the BZ type
of clipped waveforms. If there are two or more consecutive clipped points, at least one
side of the clipped point has zero amplitude. It is difficult to judge whether a zero value
sample with neighboring samples of different signs is clipped or not. However, the high
sampling frequency (100 Hz) in our data set could rule out most such cases.

4.4 An Algorithm for Detecting Clipped Waveform
Based on the features of clipped waveforms listed in section 4.3, we propose an
algorithm for clipped waveform detection with the following steps: 1) Treat waveforms
with maximum amplitudes larger than a given threshold (e.g. 60% of the Observed
Range) as potentially clipped waveforms. 2) For each potential clipped waveform,
assume that points with zero amplitude, between the first and last samples with absolute
values larger than 0.5 of the waveform peak amplitude, are tentatively clipped. 3)
Require that clipped samples are bounded on both sides by points with the same sign or
another zero. 4) Require that the maximum and minimum values of +/- 10 samples on
either side of the candidate clipped points are larger than 0.8 of the waveform peak
amplitude.
Figure 4.3 presents the statistics of the recorded maximum amplitudes on the
vertical component of fault zone station VO. In agreement with the Gutenberg-Richter
statistics, the number of observed waveforms generally decreases with increasing
97
maximum amplitude. Because of clipping, however, the population increases
dramatically beyond a certain amplitude (e.g. 6×10
4
in Figure 4.3). A reasonable choice
for a clipping threshold is somewhat below the end of the power-law range in the
frequency-size statistics of the recorded maximum amplitudes.

Figure 4.3. The distribution of maximum amplitude within 24,135 waveforms recorded by station VO on
the vertical component. The vertical dashed lines represent 98%, 95%, 90%, 85%, …, 45%, 40% of the
Observed Range.

To implement the clipped waveform detection algorithm, we choose a series of
thresholds from 98% to 40% of the Observed Range on a given component of a given
station (dashed lines in Figure 4.3). If the algorithm works properly, the numbers of
detected clipped waveforms will first increase and then become saturated with the
decreasing of the assumed threshold. Figure 4.4 demonstrates that such a pattern is
obtained for waveforms recorded by the vertical component of station VO. Application
98
of this algorithm to the data recorded by all stations indicates that as the threshold drops
below approximately 65%, the numbers of detected clipped waveforms saturate for all
stations and instrument components.
Figure 4.4. The variation of detected clipped waveforms with different amplitude thresholds, using data
from 24,135 waveforms recorded by station VO on the vertical component. Triangles are candidate
clipped waveforms with maximum amplitudes above threshold. Squares are detected clipped waveforms
that satisfy all the algorithm criteria.
99

Observed Range
(counts) Station Component
Min Max
Total
waveforms
Clipped
waveforms
E -83220.7 83218.2 281
N -83220.7 83218.2 233

BU
Z -83220.7 83218.2
14496
163
E -83220.7 83218.2 625
N -83220.7 83218.2 485

BV
Z -83220.7 83218.2
19980
530
E -83220.7 83218.2 375
N -83220.7 83218.2 437

CH
Z -83220.7 83218.2
16922
351
E -83220.7 83218.2 1360
N -83220.7 83218.2 809

FI
Z -83220.7 83218.2
21046
579
E -83220.7 83218.2 1178
N -83220.7 83218.2 1341

FP
Z -83220.7 83218.2
21394
717
E -4918765.5 4815583.0 0
N -5204982.5 5131400.5 0

GE
Z -4781097.5 4635998.5
21569
0
E -83220.7 83218.2 1028
N -83220.7 83218.2 1427

LS
Z -83220.7 83218.2
24407
521
E -332920.2 332889.7 39
N -332920.2 332889.8 28

MO
Z -332920.2 332889.8
19793
11
E -83220.7 83218.2 2680
N -83220.7 83218.2 3784

VO
Z -83220.7 83218.2
24135
1345
E -83220.7 83218.2 971
N -83220.7 83218.2 633

WF
Z -83220.7 83218.2
22100
541
Table 4.1. A summary of the observed ranges, total numbers of waveforms and numbers of detected
clipped waveforms for all stations and components.

100
Table 4.1 summarizes basic features of the recorded data and results of applying the
detection algorithm to the data set. Among the 26,080 aftershocks, 4,177 generated at
least one clipped waveform (Figure 4.5). The total number of detected clipped
waveforms is 22,472 among all 617,526 waveforms. The intensity of clipping generally
decreases with distance from the fault, and the only station that apparently did not
record any clipped waveform is station GE. In the other stations, the horizontal
components generally have more clipped waveforms than the vertical component,
although all components have similar Observed Ranges at each station. For stations
with the same Observed Range, the percentage of clipped waveforms varies from 1% on
the vertical component of station BU to 16% on the EW component of station VO.

Figure 4.5. Distribution of ~4,200 events with detected clipped waveforms (red circles) among the
~26,000 aftershocks (green circles) recorded by the network. The size of each circle scales with event
magnitude as shown by the legend.

Yang et al (2009) found that the observed waveform spectral energy associated
with the site effects of eight stations decreases following the order: VO, FP, FI, WF, LS,
101
BV, CH and BU. Figure 4.6 shows relations between event magnitudes and hypocentral
distances producing detected clipped waveforms at each of these eight stations. A
comparison of the results with the amplitude of the site spectra at the various stations
(Yang et al., 2009) illustrates, as expected, that the intensity of waveform clipping is
correlated strongly with the amplitude of the local site effects.

Figure 4.6. The magnitude and distance relations of detected clipped waveforms on the vertical
components of the various stations (denoted by 2 letters at the top-left corners). The number of detected
clipped waveforms decreases following the order: VO, FP, FI, WF, LS, BV, CH and BU.

One approach to quantitatively compare the severity of clipping is to count the
number of consecutive clipped points on a waveform. Figure 4.7 gives the probability
density distribution of clipped waveforms with different numbers of consecutive clipped
points among all detected clipped waveforms. The dominant cases of clipped
waveforms are seen to be associated with single and two consecutive clipped points
102
(45%, 53% of all clipped waveforms respectively). There are 1.2% cases of clipped
waveforms with three consecutive clipped samples and less than 0.8% cases with four
or more consecutive clipped samples.

Figure 4.7. Probability density of clipped waveforms according to the maximum numbers of consecutive
clipped points on each clipped waveform.

4.5 The Kriging and Similar Waveform Correction Methods
In this section we discuss two methods that could be used to correct less severely
clipped waveforms. Interpolation is commonly used to replace missing data points (e.g.
Kokaram et al. 1995; Karabulut and Bouchon, 2007). One of the most used
interpolation techniques in geo-statistics is the Kriging method. Below we implement
the Kriging interpolation with the DACE package (Sacks et al, 1989).
103
Various studies demonstrated that in many cases the seismicity contains subsets of
events with very similar locations and focal mechanisms. These events, referred to as
repeating earthquakes, generate very similar waveforms with high CCC values, which
in some cases match almost wiggle by wiggle (e.g., Poupinet et al., 1984; Nadeau et al.,
1994; Schaff et al., 1998; Peng and Ben-Zion, 2006). In cases where sufficiently similar
unclipped waveforms exist, scaled-up portions of similar unclipped waveforms may be
used to replace missing samples in clipped waveforms.
The correction of clipped points with the similar waveforms method is done with
the following steps:
1) For every clipped point of the BZ type in each clipped waveform, we first replace
its amplitude with the maximum observed value for that instrument component, and use
the same sign as that of one of its neighboring unclipped points with larger absolute
amplitude (Figure 4.8a). We then search for similar unclipped waveforms in the data set
and use in the next step the unclipped waveform with the highest CCC value (Figure
4.8b).
2) The similar unclipped waveform is aligned with the clipped waveform and
scaled-up to match the amplitudes of the data points with the clipped waveform (Figure
4.8c). The employed scaling-up factor is the median value of the amplitude ratios
between all the unclipped points of the clipped seismogram and the corresponding
points of the unclipped similar waveform. The values of clipped points are then
replaced with the scaled-up versions of the corresponding points on the unclipped
similar waveform (Figure 4.8d).
104

Figure 4.8. Example of correcting clipped waveform using the similar waveform method. (a) Replace the
amplitudes of clipped points with the maximum observed value. (b) Find a similar unclipped waveform
with the highest CCC (here CCC=0.96). (c) Project the similar unclipped waveform (red-dashed) onto the
clipped waveform (blue). (d) A comparison of the clipped (blue) and corrected (red) seismograms.

4.6 A Comparison of Two Correction Methods
In the following we compare the performance of the two correction methods using
artificially clipped waveforms associated with cases of 1-9 consecutive clipped samples
(Figure 4.7), which represent 99.9% of all clipped waveforms detected in our data set.
To prepare artificially clipped waveform, we use unclipped waveforms and clip samples
105
with the largest amplitude and in some cases also 1-8 neighboring points. We use
earthquake clusters with high CCC values identified by Peng and Ben-Zion (2005,
2006). The systematic tests employ 10-second waveforms with all the main seismic
phases, recorded on the horizontal component that usually contain the largest
amplitudes. Similar waveforms of one cluster recorded by station VO on the EW
component are displayed in Figure 4.9. Among all unclipped similar waveforms, we
select the waveform with the largest amplitude to be target waveform. On the selected
target waveform, we artificially clip a number of points with amplitude larger than a
given threshold.
106

Figure 4.9. (a) A selected large amplitude waveform (top trace) and its similar waveforms arranged
following the CCC values above traces (b) A stacking of all normalized similar waveforms in (a).

To compare cases with single clipped sample, we consider the sample with the
largest amplitude on the select unclipped waveform to be a clipped sample (solid blue
square in Figure 4.10a). We calculate CCC values between the artificially clipped
waveform and each of its similar waveforms in the 10-second time window excluding
the clipped sample, and then correct the amplitude of the clipped sample from similar
waveforms with different CCC values (open green squares in Figure 4.10a-i). We also
107
apply the Kriging method to interpolate for the amplitude of the clipped sample using
data with 17 samples at each side adjacent to the clipped part (red curve in Figure
4.10a-i). To compare cases with two consecutive clipped samples, we consider the
sample with largest amplitude and one of its neighboring points that has larger
amplitude to be clipped samples (solid blue squares in Figure 4.10b). For three and
more consecutive clipped samples, we consider the sample with the largest amplitude
and its neighboring samples as clipped samples (solid blue squares in Figures 4.10c-i).
In all cases the amplitudes of the artificially clipped samples are larger than the
amplitudes of unclipped samples along the waveform.

Figure 4.10. Corrections with the two methods for cases of 1-9 consecutive clipped points. The black
solid circles are original unclipped samples. The blue solid squares are artificially clipped samples. The
red curves show the correction with the Kriging method for each case. The green squares are corrections
with the similar waveform method.
108

Figure 4.10 shows a comparison of the two methods using example artificially
clipped waveform with 1-9 consecutive clipped points. With more consecutive points
clipped, additional information associated with the waveform structure is lost, the
interpolated data becomes further removed from the original amplitudes of the
artificially clipped samples, and the corrected interpolated results become less accurate.
With the similar waveform method, the CCC values between a clipped waveform and
its similar waveforms are nearly the same for cases with different number of
consecutive clipped points.
We define the error associated with the corrections to be:

err = max( log
10
A
org
− log
10
A
cor
) , (4.1)
where is the original amplitude of the selected clipped point, and is the
corrected amplitude. For case with two or more consecutive clipped points, the error in
Eq (4.1) is defined as the largest correction error since such error is associated with the
peak position along the consecutive clipped samples. The employed log scale in Eq
(4.1) produces a correspondence between the correction errors and uncertainty on the
local Richter magnitude scale.
We apply both methods to correct a group of 1,600 artificially clipped waveforms
with 1-9 consecutive clipped points, and calculate the corresponding correction errors
using Eq (4.1). We compare both methods using the median values of the correction
errors in each case. The median values of the errors associated with the similar
109
waveform method are calculated in every 0.01 CCC bin (red curves in Figure 4.11). For
each method, we also calculate the 97.5% confidence level for reference.

Figure 4.11. A comparison between errors of the two correction methods for cases with different numbers
of consecutive clipped points (name at the up-left corner). In each panel, the blue dots are the correction
errors of 1,600 selected waveforms from their similar waveforms. The red and gray curves are the median
values and 97.5% confidence levels of the correction errors for the similar waveform method in every
0.01 CCC value interval respectively. The green solid and dashed lines are respectively, the median value
and 97.5% confidence level of the correction errors for the Kriging method. For cases of 6-9 consecutive
clipped points, the median correction errors are larger than 0.5.

For cases of 1-2 consecutive clipped samples, the Kriging method always performs
better than the similar waveform method. For three clipped samples, the similar
waveform method performs better with CCC value larger than 0.96. For four
consecutive clipped samples, the similar waveform method performs better with CCC
value larger than 0.91. For five consecutive clipped samples, the similar waveform
110
method performs better with CCC value larger than 0.88. This also holds for cases of
six or more consecutive clipped samples, although the amount of available data for such
tests is considerably smaller. Table 4.2 and Figure 4.12 summarize the comparison
results extracted from Figure 4.11. The green line in Figure 4.12 separates regions in
parameter-space, spanned by the number of consecutive slipped points and CCC value
of similar unclipped events, where the different methods produce better correction
results. The errors associated with corrections are denoted by colors, and the regions
that produce errors less than 0.3 (which would lead to reasonably small error in
magnitude derivation) are un-hatched.

Kriging Method
Cases
Median 97.5% confidence level
Similar Method
(with smaller median
value)
Single 0.0007 0.012 NA
2-Consecutive 0.009 0.2 NA
3-Consecutive 0.09 1.0 CCC > 0.96
4-Consecutive 0.21 1.6 CCC > 0.91
5-Consecutive 0.29 1.7 CCC > 0.88
6-Consecutive 0.7 2.2 CCC > 0.83
Table 4.2. A Statistical summary of correction errors for methods comparison

111

Figure 4.12. A comparison of the performance of the two correction methods. Better performance regions
are separated by the green line. The correction errors is colored with scale to the right. The correction
error in the un-hatched region is smaller than 0.3.

4.7 Discussion and Conclusions
We design and implement an automatic detection algorithm to identify clipped
waveforms in a near-fault data set generated by 26,080 earthquakes, where amplitudes
of clipped points set to zero instrumentally. Applications of the algorithm to other types
of instrumental clipping (e.g., Flat-Top) will require some adjustments, but the general
statistical features of clipped waveforms should remain the same. In our case, the
112
algorithm indicates that 3.6% of the waveforms are clipped, with the intensity of
clipping increasing with proximity to the fault and amplitude of local site effects.
To increase the range of available data it may be useful to correct some clipped
waveforms. Toward this end, we compare two possible correction methods associated
with the Kriging interpolation and using scaled-up version of similar unclipped
waveforms. The availability of the latter generally increases with increasing geometrical
simplicity of the fault (e.g., Ben-Zion and Sammis, 2003). As examples, along the
Parkfield section of the San Andreas fault more than 50% of the ongoing seismicity
belong to clusters of repeating events that have CCC values larger than 0.9 (e.g.,
Nadeau et al., 1994). In contrast, along the San Jacinto fault in southern California only
about 2% of the seismicity belong to clusters of repeating events with CCC > 0.9 (Aster
and Scott, 1993). For the seismicity along the Karadere-Duzce of the North Anatolian
fault used in this study, about 18% of the seismicity belong to clusters of repeating
events with CCC > 0.9 (Peng and Ben-Zion, 2005, 2006).
Typical median value of magnitude uncertainty in the Northern California Seismic
Network (NCSN) earthquake catalog of California is about 0.3 (Werner and Sornette,
2008). Since the correction error associated with Eq (4.1) is based on the log amplitude,
errors less than about 0.3 should produce errors in magnitude estimates that are within
the typical range. Our synthetic tests based on data with 100 Hz sampling frequency and
1-9 consecutive clipped samples can be summarized as follows. For clipped waveform
with 1-2 clipped samples, the Kriging method always performs better. For clipped
waveform with 3-5 consecutive clipped samples, the Kriging method produces
corrections with smaller errors unless unclipped similar waveforms with high CCC
113
values (0.88-0.96) are available. For clipped waveforms with 6 or more consecutive
clipped samples, the similar waveform method performs better but unclipped
waveforms with CCC > 0.85 are needed to produce errors smaller than 0.3. The CCC
values separating cases where the different methods perform better and calculated errors
will change somewhat for data with different frequency content and sampling rates.

114
Chapter 4 References
Aster, R.C. & Scott, J. 1993. Comprehensive characterization of waveform similarity
in microearthquake data sets, Bull. Seism. Soc. Am., 83, 1307-1314.

Ben-Zion, Y., Armbruster, J.G., Ozer, N., Micheal, A.J., Baris, S., Aktar, M., Peng, Z.,
Okaya, D. & Seeber, L., 2003. A shallow fault-zone structure illuminated by trapped
waves in the Karadere-Düzce branch of the North Anatolian Fault, western Turkey,
Geophys. J. Int., 152, 699-717.

Ben-Zion, Y. & Sammis, C.G., 2003. Characterization of Fault Zones, Pure Appl.
Geophys., 160, 677-715.

Karabulut H. and Bouchon, M., 2007. Spatial variability and non-linearity of strong
ground motion near a fault, Geophys. J. Int., 170, 262-274.

Kokaram, A.C., Morris, R.D., Fitzgerald, W.J. & Rayner, P.J.W., 1995. Interpolation of
missing data in image sequences, Image Processing, IEEE Transactions on, vol.4,
no.11, 1509-1519.

Nadeau, R.M., Antolik, M., Johnson P., Foxall, W. and McEvilly T.V., 1994.
Seismological studies at Parkfield III: microearthquake clusters in the study of fault-
zone dynamics, Bull. Seism. Soc. Am., 84, 247-263.

Peng, Z. and Y. Ben-Zion, 2005. Spatio-temporal variations of crustal anisotropy from
similar events in aftershocks of the 1999 M7.4 İzmit and M7.1 Düzce, Turkey,
earthquake sequences, Geophys. J. Int., 160, 1027-1043.

Peng, Z. and Ben-Zion Y., 2006. Temporal changes of shallow seismic velocity around
the Karadere-Düzce branch of the north Anatolian fault and strong ground motion,
Pure and Appl. Geophys., 163, 567-600.

Poupinet, G., Ellsworth W.L. & Frechet J., 1984. Monitoring velocity variations in the
crust using earthquake doublets: an application to the Calaveras fault, California, J.
Geophys. Res., 89, 5719-5731.

Sacks, J., Welch, W.J., Mitchell, T.J. & Wynn, H.P., 1989. Design and Analysis of
Computer Experiments, Statistical Science, 4, No. 4, 409-423.

Schaff D.P., Beroza G.C. & Shaw B.E., 1998. Postseismic response of repeating
aftershocks, Geophysical Research Letters, 25, 4549.

115
Seeber, L., Armbruster, J.G., Ozer, N., Aktar, M., Baris, S., Okaya, D., Ben-Zion, Y. &
Field, N., 2000, The 1999 Earthquake Sequence Along the North Anotolian
Transform at the Juncture Between the Two Main Ruptures: Published in Special
Volume The 1999 İzmit and Düzce Earthquakes: Preliminary Results, Eds. Barka,
A., Kozaci, O., Akyuz, S., and Altunel, S., ITU, Istanbul.

Shearer, P., Prieto, G. & Hauksson, E., 2006. Comprehensive analysis of earthquake
source spectra in southern California, J. Geophys. Res., 111, B06303,
doi:10.1029/2005JB003979.

Yang, W., Peng, Z. & Ben-Zion, Y., 2009. Variations of strain drops in aftershocks of
the 1999 Izmit and Düzce earthquakes along the Karadere-Düzce branch of the
North Anatolian fault, Geophys. J. Int., 177, 235-246.

Werner, M. J. & Sornette D., 2008. Magnitude uncertainties impact seismic rate
estimates, forecasts, and predictability experiments, J. Geophys. Res., 113, B08302,
doi:10.1029/2007JB005427.
116
Alphabetized Bibliography
Abercrombie, R., 1995. Earthquake source scaling relationships from −1 to 5 M
L
using
seismograms recorded at 2.5 km depth, J. Geophys. Res., 100, 24,015–24,036.

Aki, K., 1982. Scattering and Attenuation, Bull. Seism. Soc. Am., 72, 319-330.

Aki, K. & Richards, P.G., 2002. Quantitative Seismology, University Science Books.

Allmann, B.P. & Shearer, P.M., 2007. Spatial and temporal stress drop variations in
small earthquakes near Parkfield, California, J. Geophys. Res., 112, B04305,
doi:10.1029/2006JB004395.

Allmann, B.P. & Shearer, P.M., 2009. Global variations of stress drop for moderate to
large earthquakes, J. Geophys. Res., 114. B01310, doi:10.1029/2008JB005821.

Andrews, D.J., 1986. Objective determination of source parameters and similarity of
earthquakes of different size, in Earthquake Source Mechanics (S. Das, J.
Boatwright, and C.H. Scholz, editors), American Geophysical Union Monograph, 37,
259–267.

Aster, R.C. & Scott, J. 1993. Comprehensive characterization of waveform similarity
in microearthquake data sets, Bull. Seism. Soc. Am., 83, 1307-1314.

Bates, D.M. & Watts, D.G., 1988. Nonlinear Regression Analysis and Its Applications,
Wiley.

Ben-Zion, Y., 1996. Stress, slip, and earthquakes in models of complex single-fault
systems incorporating brittle and creep deformations, J. Geophys. Res., 101, 5677-
5706.

Ben-Zion, Y., 2003. Key Formulas in Earthquake Seismology, International handbook
of earthquake and engineering seismology, Part B, 1857-1875, Academic Press.

Ben-Zion, Y., 2008. Collective Behavior of Earthquakes and Faults: Continuum-
Discrete Transitions, Progressive Evolutionary Changes and Different Dynamic
Regimes, Rev. Geophysics, 46, RG4006, doi:10.1029/2008RG000260.

Ben-Zion, Y., Armbruster, J.G., Ozer, N., Micheal, A.J., Baris, S., Aktar, M., Peng, Z.,
Okaya, D. & Seeber, L., 2003. A shallow fault-zone structure illuminated by trapped
waves in the Karadere-Düzce branch of the North Anatolian Fault, western Turkey,
Geophys. J. Int., 152, 699–717.

117
Ben-Zion, Y. & Lyakhovsky, V., 2006. Analysis of aftershocks in a lithospheric model
with seismogenic zone governed by damage rheology, Geophys. J. Int., 165, 197-
210.

Ben-Zion, Y. & Sammis, C.G., 2003. Characterization of Fault Zones, Pure Appl.
Geophys., 160, 677-715.

Ben-Zion, Y. & Zhu, L., 2002. Potency-magnitude scaling relations for Southern
California earthquakes with 1.0 < M
L
< 7.0, Geophys. J. Int.,148, F1–F5.

Berckhemer, H., 1962. Die Ausdehnung der Bruchfläche im Erdbebenherd und ihr
Einfluss auf das seismische Wellenspektrum, Gerlands Beitr. Geophys., 71, 5–26.

Boatwright, J. & Fletcher J.B., 1984. The partition of radiated energy between P and S
waves, Bull. Seism. Soc. Am., 74, 361-376.

Boatwright, J., Fletcher, J.B., & Fumal, T.E., 1991. A general inversion scheme for
source, site, and propagation characteristics using multiply recorded sets of
moderate-sized earthquakes, Bull. Seismol. Soc. Am., 81, 1754–1782.

Bohnhoff, M., Grosser, H. & Dresen, G., 2006. Strain partitioning and stress rotation at
the North Anatolian Fault after the 1999 İzmit Mw=7.4 earthquake, Geophys. J. Int.,
166, 373-385.

Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from
earthquakes, J. Geophys. Res., 75, 4997–5009.

Bulut, F., Bohnhoff, M., Aktar, M. & Dresen, G., 2007. Characterization of aftershock-
fault plane orientations of the 1999 Izmit (Turkey) earthquake using high-resolution
aftershock locations, Geophys. Res. Lett., 34, L20306, doi:10.1029/2007GL031154.

Burdick, L.J., 1982. Comments on "the corner frequency shift, earthquake source
models, and Q," by T. C. Hanks, Bull. Seism. Soc. Am., 72, 1419-1426.

Castro, R.R., Anderson, J.G. & Brune, J.N., 1991. Origin of high P/S spectral ratios
from the Guerrero accelerogram array, Bull. Seism. Soc. Am., 81, 2268-2288.

Creager, K.C., 1997, Coral, Seismol. Res. Lett., 68, 269–271.

Dahlen, F.A., 1974. On the ratio of P-wave to S-wave corner frequencies for shallow
earthquake sources, Bull. Seism. Soc. Am., 64, 1159-1180.

Davis, S.D. & Frohlich, C., 1991. Single-link cluster analysis of earthquake aftershocks:
decay laws and regional variations, J. geophys. Res., 96, 6335–6350.

118
Enescu, B., Mori, J. & Miyazawa, M., 2007. Quantifying early aftershock activity of the
2004 mid-Niigata Prefecture earthquake (Mw6.6), J. Geophys. Res., 112, B04310,
doi:10.1029/2006JB004629.

Felzer, K.R., Abercrombie, R.E. & Ekstrom, G., 2004. A common origin for
aftershocks, foreshocks and multiplets, Bull. Seismol. Soc. Am., 94, 88-98.

Felzer, K.R. & Brodsky, E. E., 2006. Decay of aftershock density with distance
indicates triggering by dynamic stress, Nature, 441, 735-738.

Fletcher J.B. & McGarr, A., 2006. Distribution of stress drop, stiffness, and fracture
energy over earthquake rupture zones, J. Geophys. Res., 111, B03312,
doi:10.1029/2004JB003396.

Furuya, I., 1969. Predominant period and magnitude, J. Phys. Earth 17, 119-126.

Gross, S.J. & Kisslinger, C., 1994. Tests of models of aftershock rate decay. Bull.
Seismolo. Soc. Am., 84, 1571-1579.

Hanks, T.C., 1981. The corner frequency shift, earthquake source models, and Q, Bull.
Seism. Soc. Am., 71, 597-612.

Hanks, T.C., 1982. Reply to "comments on 'the corner frequency shift, earthquake
source models, and Q'", Bull. Seism. Soc. Am., 72, 1433-1444.

Hanks, T.C., and Wyss M., 1972. The use of body-wave spectra in the determination of
seismic source parameters, Bull. Seism. Soc. Am., 62, 561-589.

Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y. & Lockner, D., 2004. A Visco-
elastic damage model with applications to stable and unstable fracturing, Geophys.
J. Int., 159, 1155–1165, doi: 10.1111/j.1365-246X.2004.02452.x.

Helmstetter, A., Kagan, Y.Y. & Jackson, D.D., 2005. Importance of small earthquakes
for stress transfers and earthquake triggering, J. Geophys. Res., 110, B05S08,
doi:10.1029/2004JB003286.

Hill, D.P. & Prejean, S.G., 2007. Dynamic Triggering, in Earthquake Seismology -
Treatise on Geophysics ed. H. Kanamori and G. Schubert, V4, 257-291, Elsevier,
Amsterdam.

Hough, S.E., 1997. Empirical Green's function analysis: Taking the next step, J.
Geophys. Res., 102, 5369–5384.

119
Ide, S., & Beroza G.C., 2001. Does Apparent Stress Vary with Earthquake Size?,
Geophys. Res. Lett., 28(17), 3349–3352.

Imanishi, K. & Ellsworth, W.L., 2006. Source scaling relationships of microearthquakes
at Parkfield, CA, determined using the SAFOD pilot hole array. In Abercrombie,
R.E., McGarr, A., Di Toro, G., and Kanamori, H. (Eds.), Earthquakes: Radiated
Energy and the Physics of Faulting, Geophys. Monogr. 170, Washington D.C.
(American. Geophysical Union), 81–80.

Karabulut, H. & Bouchon, M., 2007. Spatial variability and non-linearity of strong
ground motion near a fault Geophys. J. Int., 170, doi:10.1111/j.1365-
246X.2007.03406.x

Kisslinger, C. & Jones, L.M., 1991. Properties of aftershock sequences in southern
California. J. Geophys. Res., 96, 11,947-11,958.

Kokaram, A.C., Morris, R.D., Fitzgerald, W.J. & Rayner, P.J.W., 1995. Interpolation of
missing data in image sequences, Image Processing, IEEE Transactions on, vol.4,
no.11, 1509-1519.

Langston, C.A., 1982. Comments on "the corner frequency shift, earthquake source
models, and Q," by T. C. Hanks, Bull. Seism. Soc. Am., 72, 1427-1432.

Lewis, M.A., Peng Z., Ben-Zion Y., & Vernon F.L., 2005. Shallow Seismic Trapping
Structure in the San Jacinto fault zone near Anza, California , Geophys. J. Int., 162,
867.881, doi:10.1111/j.1365-246X.2005.02684.x, 2005.

Lohman, R.B. & McGuire, J.J., 2007. Earthquake swarms driven by aseismic creep in
the Salton Trough, California. J. Geophys. Res., 112, B04405,
doi:10.1029/2006JB004596.

Lyakhovsky, V., Ben-Zion, Y. & Agnon, A., 1997. Distributed Damage, Faulting, and
Friction, J. Geophys. Res., 102, 27635-27649.

Lyakhovsky, V. & Ben-Zion Y., 2008. Scaling relations of earthquakes and aseismic
deformation in a damage rheology model, Geophys. J. Int., 172, 651-662, doi:
10.1111/j.1365-246X.2007.03652.x.

McGuire, J.J., Boettcher, M.S. & Jordan, T.H., 2005. Foreshock sequences and short-
term earthquake predictability on East Pacific Rise Transform Faults, Nature, 434,
457–461.

Mogi, K., 1967. Earthquakes and fractures. Tectonophysics, 5, 35-55.

120
Nadeau, R.M., Antolik, M., Johnson P., Foxall, W. & McEvilly T.V., 1994.
Seismological studies at Parkfield III: microearthquake clusters in the study of fault-
zone dynamics, Bull. Seism. Soc. Am., 84, 247-263.

Nauteau, C., Shebalin, P. & Holschneider, M., 2002. Temporal limits of the power law
aftershock decay rate, J. Geophys. Res., 107, B12, 2359,
doi:10.1029/2002JB001868.

Omori, F., 1894. On the aftershocks of earthquakes, J. Coll. Sci., Imp. Univ. Tokyo, 7,
111-200.

Otsuka, M., 1985. Studies on aftershock sequences-Part 1. Physical interpretation of
Omori formula, Sci. Rep. Shimabara Earthq. Volcano Obs., 12, 11-20 (in Japanese).

Park, J., Lindberg, C.R. & Vernon, F.L., 1987. Multitaper spectral analysis of high
frequency seismograms, J. Geophys. Res., 92, 12,675–12,648.

Peng, Z., Ben-Zion Y., Michael A.J. & Zhu L., 2003. Quantitative analysis of seismic
trapped waves in the rupture zone of the 1992 Landers, California earthquake:
Evidence for a shallow trapping structure, Geophys. J. Int., 155, 1021–1041.

Peng, Z., & Ben-Zion, Y., 2004. Systematic analysis of crustal anisotropy along the
Karadere-Düzce branch of the north Anatolian fault, Geophys. J. Int., 159, 253–274.

Peng, Z. & Ben-Zion Y., 2005. Spatio-temporal variations of crustal anisotropy from
similar events in aftershocks of the 1999 M7.4 İzmit and M7.1 Düzce, Turkey,
earthquake sequences, Geophys. J. Int., 160, 1027–1043.

Peng, Z. & Ben-Zion, Y., 2006. Temporal changes of shallow seismic velocity around
the Karadere-Düzce branch of the north Anatolian fault and strong ground motion,
Pure and Appl. Geophys. 163, 567–600.

Peng, Z., Vidale, J.E. & Houston, H., 2006. Anomalous early aftershock decay rates of
the 2004 M6 Parkfield earthquake, Geophys. Res. Lett., 33, L17307,
doi:10.1029/2006GL026744.

Poupinet, G., Ellsworth W.L. & Frechet J., 1984. Monitoring velocity variations in the
crust using earthquake doublets: an application to the Calaveras fault, California, J.
Geophys. Res., 89, 5719-5731.

Prieto, G., Shearer, P.M., Vernon, F.L. & Kilb, D., 2004. Earthquake source scaling and
self–similarity estimation from stacking P and S spectra, J. Geophys. Res., 109,
B08310, doi:10.1029/2004JB003084.

121
Prieto, G.A., Thomson D.J., Vernon F. L. & Shearer P. M., 2006. How much does P-
wave coda bias S-wave spectral estimates?, Seismol. Res. Lett., 77, 167.

Madariaga, R., 1976. Dynamics of an expanding circular fault, Bull. Seismol. Soc. Am.,
66, 639–666.

Mayeda, K., Hofstetter, A., O'Boyle, J.L. & Walter, W.R., 2003. Stable and
Transportable Regional Magnitudes Based on Coda-Derived Moment-Rate Spectra,
Bull. Seism. Soc. Am., 93, 224-239.

Molnar, P., Tucker, B.E. & Brune, J.N., 1973. Corner frequencies of P and S waves and
models of earthquake sources, Bull. Seism. Soc. Am., 63, 2091-2104.

Mueller, C.S., 1985. Source pulse enhancement by deconvolution of an empirical
Green's function, Geophys. Res. Lett., 12, 33–36.

Rolandone, F., Bürgmann R. & Nadeau R.M., 2004. The evolution of the seismic-
aseismic transition during the earthquake cycle: Constraints from the time-
dependent depth distribution of aftershocks, Geophys. Res. Lett., 31, L23610,
doi:10.1029/2004GL021379.

Sacks, J., Welch, W.J., Mitchell, T.J. & Wynn, H.P., 1989. Design and Analysis of
Computer Experiments, Statistical Science, 4, No. 4, 409-423.

Savage, J.C., 1972. Relation of Corner Frequency to Fault Dimensions, J. Geophys.
Res., 77, 3788–3795.

Sato, T. & Hirasawa, T., 1973. Body wave spectra from propagating shear cracks, J.
Phys. Earth, 21, 415-431.

Schaff D.P., Beroza G.C. & Shaw B.E., 1998. Postseismic response of repeating
aftershocks, Geophysical Research Letters, 25, 4549.

Seeber, L., Armbruster, J.G., Ozer, N., Aktar, M., Baris, S., Okaya, D., Ben-Zion, Y. &
Field, N., 2000. The 1999 Earthquake Sequence Along the North Anotolian
Transform at the Juncture Between the Two Main Ruptures: Published in Special
Volume The 1999 İzmit and Düzce Earthquakes: Preliminary Results, Eds. Barka,
A., Kozaci, O., Akyuz, S., and Altunel, S., ITU, Istanbul.

Sengor, A.M.C., Tuysuz, O., Imren, C., Sakınc¸ M., Eyidogan, H., Gorur, N., Le
Pichon, X. & Rangin, C., 2005. The North Anatolian Fault: A New Look, Annu.
Rev. Earth Planet. Sci., 33, 37–112, doi: 10.1146/annurev.earth.32.101802.120415.

122
Shearer, P., Hauksson E. & Lin, G., 2005. Southern California hypocenter relocation
with waveform cross-correlation, Part 2: Results using source-specific station terms
and cluster analysis, Bull. Seismol. Soc. Am., 95, 904-915, doi:
10.1785/0120040168.

Shearer, P., Prieto, G. & Hauksson, E., 2006. Comprehensive analysis of earthquake
source spectra in southern California, J. Geophys. Res., 111, B06303,
doi:10.1029/2005JB003979.

Thomson, D.J., 1982. Spectrum estimation and harmonic analysis, Proceeding of the
IEEE., 70, 1055–1096.

Trifunac, M.D., 1972. Stress estimates for the San Fernando, California earthquake of
February 9, 1971: Main event and thirteen aftershocks, Bull. Seism. Soc. Am., 62,
721-750.

Utsu, Y., Ogata, Y. & Matsu’uara, R.S., 1995. The centenary of the Omori formula for
a decay law of aftershock activity, J. Phys. Earth, 43, 1–33.

Warren, L.M. & Shearer, P.M., 2002. Mapping lateral variations in upper mantle
attenuation structure by stacking P and PP spectra, J. Geophys. Res., 107(B12),
2342, doi: 10.1029/2001JB001195.

Werner, M.J. & Sornette D., 2008. Magnitude uncertainties impact seismic rate
estimates, forecasts, and predictability experiments, J. Geophys. Res., 113, B08302,
doi:10.1029/2007JB005427

Wiemer, S. & Wyss, M., 2000. Minimum Magnitude of Completeness in Earthquake
Catalogs: Examples from Alaska, the western United States and Japan, Bull.
Seismol. Soc. Am., 90, 4, 859-869.

Wu, C., Peng Z. & Ben-Zion Y., 2009. Non-linearity and temporal changes of fault
zone site response associated with strong ground motion, Geophys. J. Int., 176,
265-278, doi: 10.1111/j.1365-246X.2008.04005.x.

Wyss, M. & Hanks, T.C., 1972. The source parameters of the San Fernando earthquake
inferred from teleseismic body waves, Bull. Seism. Soc. Am., 62, 591-602.

Yang, W. & Ben-Zion, Y., 2009. An algorithm for detecting clipped waveforms and
suggested correction procedures. Submit to Seismol. Res. Lett.,

Yang, W., Peng, Z. & Ben-Zion, Y., 2009. Variations of strain drops in aftershocks of
the 1999 Izmit and Düzce earthquakes along the Karadere-Düzce branch of the
North Anatolian fault, Geophys. J. Int., 177, 235-246.

123
Zöller, G., Holschneider, M. & Ben-Zion, Y., 2005. The role of heterogeneities as a
tuning parameter of earthquake dynamics, Pure Appl. Geophys., 162, 1027-1049.

Abstract (if available)

Abstract We study the statistical relationship between aftershock productivity and background geophysical features using an earthquake catalog with ~340,000 events that occurred between 1984 and 2004 in the Southern California region. We analyze properties of individual aftershock sequences and the average properties of stacked aftershock sequences in five regions. The results indicate that productivities of the aftershock sequences are inversely correlated with heat flow and existence of deep sedimentary covers, in agreement with the damage model predictions given by Ben-Zion and Lyakhovsky (2006). Using the observed ratios of aftershock productivities, along with simple expressions based on the damage model, we estimate the seismic coupling coefficients in the different regions. We estimate seismic strain-drops from a data set with ~26,000 aftershocks that occurred along the Karadere-Düzce branch of the North Anatolian Fault after the 1999 İzmit mainshock. We use a method that is associated with separation of source, travel-time and station spectral terms, and stacking results at several stages to enhance the signal-to-noise ratio. The strain-drops are obtained by fitting iteratively the separated source spectral of 201 nearest neighboring events in different amplitude bins to the omega square source spectra model. We analyze the variations of the obtained strain-drops along the fault zone, with the depth extension. Furthermore, we analyze seismic source properties using both P and S waves on seismograms. We develop a procedure to fit for P/S corner frequency ratios together with relative strain-drop systematically. We consider five different schemes to analyze the S wave, and the obtained median values of frequency ratios from all schemes are larger than 1.5. We describe an algorithm to detect clipped waveforms in the same ~26000 data set based on statistical characteristics of observed amplitude ranges, properties of neighboring points in seismograms.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Heterogeneity of earthquake stress drops, focal mechanisms and active fault zones

PDF

Volumetric interactions between major ruptures and fault zones illuminated by small earthquake properties

PDF

Dynamic rupture processes and seismic radiation in models of earthquake faults separating similar and dissimilar solids

PDF

High-resolution imaging and monitoring of fault zones and shallow structures: case studies in southern California and on Mars

PDF

Multi-scale imaging and monitoring of crustal and fault zone structures in southern California

PDF

Multi-scale imaging of the fault zone velocity structure: double-difference tomography, inversion of fault zone headwaves, and fault zone sensitivity kernels

PDF

Stress-strain characterization of complex seismic sources

PDF

Detailed properties of seismic waveforms and earthquake sources in Southern California

PDF

Integration and validation of deterministic earthquake simulations in probabilistic seismic hazard analysis

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Determining fault zone structure and examining earthquake early warning signals using large datasets of seismograms

PDF

Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure

PDF

Observations and modeling of dynamically triggered high frequency burst events

PDF

Multi-scale imaging of major fault zones in Southern California

PDF

Statistical analyses of ambient seismic noise spectra with applications to detecting imperfectly diffuse wave field and deriving attenuation and phase velocity information

PDF

Symmetry properties, pulverized rocks and damage architecture as signatures of earthquake ruptures

PDF

Spatiotemporal variations of stress field in the San Jacinto Fault Zone and South Central Transverse Ranges

PDF

Implications of new fault slip rates and paleoseismologic data for constancy of seismic strain release and seismic clustering in the eastern California shear zone

PDF

Microseismicity, fault structure, & the seismic cycle: insights from laboratory stick-slip experiments

PDF

Elements of seismic structures near major faults from the surface to the Moho

Asset Metadata

Creator Yang, Wenzheng (author)

Core Title Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Geological Sciences

Degree Conferral Date 2009-12

Publication Date 08/18/2009

Defense Date 07/17/2009

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

Tag aftershock,clipping,heat flow,OAI-PMH Harvest,seismic source,seismicity,waveform

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ben-Zion, Yehuda (committee chair), Sammis, Charles G. (committee member), Teng, Ta-Liang (committee member), Trifunac, Mihailo D. (committee member)

Creator Email wenzhengy@gmail.com,wenzheny@usc.edu

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

Unique identifier UC1154809

Identifier etd-Yang-3197 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-253278 (legacy record id),usctheses-m2570 (legacy record id)

Legacy Identifier etd-Yang-3197.pdf

Dmrecord 253278

Document Type Dissertation

Rights Yang, Wenzheng

Type texts

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

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu