Synopses of Master Papers

Bulletin of IISEE, 47, 19-24, 2013

SOURCE PROCESS OF THE 2003 PUERTO PLATA
EARTHQUAKE USING TELESEISMIC DATA AND
STRONG GROUND MOTION SIMULATION
Fabricio Moquete Everth*
MEE11604

Supervisor: Yuji YAGI **Toshiaki YOKOI ***

Abstract
Using teleseismic data, the spatial and temporal slip distribution of the 2003 Puerto Plata earthquake
in the Dominican Republic, with the recent method, was determined. We obtained the slip distribution
over the fault plane, with the following source parameters: (strike, dip and rake) = (110.0º, 8.0º and76.6 º); seismic moment of M 0 = 1.446×1019 Nm (M w 6.7), source duration = 14 sec. These results
obtained from the inversion analysis were used as parameters for simulating strong ground motion in
the Puerto Plata area and other cities by empirical attenuation formulas for Peak Ground Acceleration(PGA) and Peak Ground Velocity (PGV) with the shortest distance to the coseismic rupture and its
moment magnitude. The fault area was estimated to be 30 km x 27 km = 810km2 from the aftershock
distribution obtained by a local seismic network. These PGA and PGV values are converted to themodified Mercalli intensity by using the empirical formula. The obtained values are quite consistent
with an isoseismal map with the registered intensity that was felt during the Puerto Plata earthquake.
Keywords: Inversion, Puerto Plata Earthquake, Source Process, Strong Ground Motion.

1. INTRODUCTION
On September 22nd, 2003, at 00:45 (4:45 UTC) the northern part of the Dominican Republic on theHispaniola Island was severely shaken by an M w 6.5 (NEIC-IRIS) earthquake (Figure 1). A very wide
damage was caused on buildings in the major cities of Puerto Plata and Santiago along with landslides
in outer areas followed by a large aftershock of M w 5.6 (USGS) about an hour after the main shock.
The seismic activity in the aftershock zone increased that year; only in the month of September 208
earthquakes were registered. Also aftershocks ranging up to nearly M w 5.0 continued for over a month,alarming local inhabitants.. The focal mechanism of the main shock and the aftershock zone, alongwith geological evidence for an active fold-thrust belt offshore of Hispaniola Island, suggests that theearthquake occurred beneath Puerto Plata, on the south dipping thrust fault located offshore (Dolanand Bowman, 2004).

Atlantic Ocean
Puerto Plata
Santiagode los C.
Santo Domingo

Figure 1. Focal mechanism of the 2003
Puerto Plata earthquake and its largest
aftershock, showing the location of thecity of Puerto Plata and the two main
cities of the country, the Capital city ofSanto Domingo and the second mostimportant Santiago de los Caballeros(Global CMT Catalog Search, 2003).

__________________________________________________________________________________* Seismological Institute of Universidad Autonoma de Santo Domingo.
** Assoc. Professor, Graduate School of Life and Environmental Sciences, University of Tsukuba, Japan.
*** Chief Research Scientist, International Institute of Seismology and Earthquake Engineering, BRI, Japan.

2. METHOD
Yagi and Fukahata (2011) proposed an improved method in which the uncertainty of the Green’s
Function is included in the inversion analysis. This new method considers the modeling error of theGreen’s function as shown below.
2.1. Observation equation
In general the observed seismic waveform for the far-field term at a station j due to shear dislocation
source on a fault plane S is given as follows:
200(𝑡, 𝜉) ∗ 𝐷̇𝑞𝑞(𝑡, 𝜉)d𝜉 + 𝑒𝑏𝑏
𝑢𝑗 (𝑡) = � � 𝐺𝑞𝑞
𝑞=1 𝑆

(1)

For Green’s function calculation, P, pP and sP phases are considered together with thereceiver response function and inelastic attenuation adding the possibility for the mechanism of thosesub-events to change (Figure2).

Figure 2. Left: Spatial basis function 𝑋𝑘 (𝜉). Right: Finite fault model with multiple point
sources of sub-events distributed on the fault plane.
The amplitude and timing of reflection phases in heterogeneous media are the mostdifficult part in the estimation of the Green’s function. The modeling error 𝛿𝛿𝑞𝑞𝑞 of the Green’sfunction is introduced as:0(𝑡) = 𝑔𝑞𝑞𝑞 (𝑡) + 𝛿𝛿𝑞𝑞𝑞 (𝑡),𝑔𝑞𝑞𝑞(2)This error adds a new term to the observation error 𝑒𝑏𝑏 as follows:
2𝐾𝐞𝑗 (𝐚) = � � 𝐏𝑞𝑞𝑞 (𝐚)δ𝐠 𝑞𝑞𝑞 + 𝐁𝑗 𝐞𝑏𝑏
𝑞=1 𝑘=1(3)2.1.1. Prior Constraints
we use smoothness of slip rate in space and time as prior constraints.
∇2 𝐷̇(𝑡, 𝜉) + 𝑒𝑠 = 0;

𝜕2
𝐷̇ (𝑡, 𝜉) + 𝑒𝑡 = 0,
𝜕𝑡 𝑡

(4)

where the errors 𝑒𝑠 and 𝑒𝑡 are Gaussian with zero and covariance 𝜌12 𝐈 and 𝜌22 𝐈, respectively.

2.1.2. Akaike’s Bayesian Information Criterion
The Akaike’s Bayesian Information Criterion (ABIC, Akaike, 1980) expression for a given initialmodel parameter vector 𝐚𝑖 is as follows:(𝛼 2 , 𝛾 2 ) = 𝑁𝑁𝑁𝑁 𝑠(𝐚∗ ) − 𝑙𝑙𝑙𝛼 2 |𝐆1 + 𝑥 2 𝐆2 | + 𝑙𝑙𝑙�𝐇 𝑡 𝐂𝑑−1 �𝛾 2 , 𝐚𝑖 �𝐇 + 𝛼 2 (𝐆1 + 𝑥 2 𝐆2 )� +log�Gd �γ2 , ai �� + C ′(5)Here we used:

𝜎𝑔2 = 𝑠(𝐚∗ )/𝑁,(6)which is derived from the necessary condition that the partial derivative of the ABIC respecting σ2gshould be zero.
To reconstruct the spatio-temporal slip-rate distribution 𝐷̇𝑞 is the interpolation using afinite number of basis functions (Figure 2 Left):𝐾𝐿𝐷̇𝑞0 (𝑡, 𝜉) =� � � 𝑎𝑞𝑞𝑞 𝑋𝑘 (𝜉)𝑇𝑙 (𝑡 − 𝑡𝑘 )𝑘=1 𝑙=1(7)

3. DATA
Teleseismic body-wave data was retrieved from the Incorporated Research Institutions for Seismology(IRIS) (http://www.iris.edu/dms/wilber.htm); the following seismic networks were selected:
 Global Telemetered Seismograph Network (USAF/USGS).
 Global Seismographic Network (GSN) stations and arrays (IRIS/USGS - IRIS/IDA).
The BHZ channels of more than 16 stations were downloaded and their quality, (e. g. thesignal to noise ratio) was checked. The number of stations chosen was 11 (Table 1 and Figure 3) andonly P-wave was picked for the analysis and only the stations with epicentral distance between 30 and90 degrees from the event’s location were selected.
Table 1. List of the 11 chosen stations.

Figure 3. Location ofthe chosen stations
on the planet and theearthquake’sepicenter location.

4. INVERSION ANALYSIS
We adopted the epicenter (19.85N 70.67W) and depth (10 km) from the USGS catalog, with a fault
mechanism parameters: strike, dip, rake = (110.0º, 10.0º, 85º) from the Global CMT catalog, withmagnitude M w 6.5 from NEIC - IRIS as an input parameter for the inversion process.
The minimum ABIC in the tested data was obtained with the depth value of 10 km. Figure 4b shows
the graphic with the ABIC values obtained when the dip angle was changed from 8 to 11 degrees at afixed depth of 10 km. With the dip of 8 degrees the ABIC becomes minimum. For this reason it was
decided that the depth should remain 10 km and the dip 8º (Figure 4, a) and b)).

Figure 4. a) ABIC variation with different depth values. b) ABIC variation withdifferent dip angle values; the red dot represents the minimum ABIC obtained.

5. RESULTS AND DISCUSSION
In the inversion method, the fault plane can be identified along with the calculation of the moment
magnitude, the dimensions of the fault, the source time function, and a detailed temporal and spatial
distribution of slip on the fault.
The time window duration for the Green’s function was 50 seconds. There is a series of
15 triangle functions for the slip rate function of every sub-fault with a rise time window of 5 seconds
each.
5.1.

Inversion Process Outcome

The results are shown in Figure 4. The focal mechanism shows the fault geometry with a rupture
duration of 14 seconds and total seismic moment of M 0 = 1.446×1019 Nm (M w 6.7) (Figure 5 Left).
This seismic moment is higher than the one calculated by the NEIC-USGS (M 0 = 7.1×1018 Nm)(M w 6.5) and the Global CMT (M 0 = 5.25×1018 Nm) (M w 6.4).
In Figure 5 (left and center) we can see that the coseismic slip extension along the strikezone coincides with the aftershock distribution of approximately 30 km along the strike zone by which thesize of the fault can be roughly estimated.
Figure 5 (Right) shows the waveform obtained from the inversion process. The black
waveform represents the observed data while the red waveform represents the synthetic waveformfrom the inversion output. As it is shown in the image, the synthetic and observed waveform matches.

Figure 5. Left: Focal mechanism, Source time function and Coseismic slip vectors distributionof the 2003 Puerto Plata earthquake. Center: Coseismic slip distribution shown on the map of
the Puerto Plata area. Right: Comparison between the observed (black) and the synthetic (red)waveform data obtained with the new formulation.

From the results of the inversion it can be seen that the slip amount is higher in the areathat surrounds the hypocenter (seismogenic zone) with a maximum slip of 0.30 m (red area) and 0.21m (yellow area). The red and yellow areas shown in Figure 5 (left) are considered as the asperities; therupture starting point (hypocenter) is almost at the center of this asperity associated with maximum
slip. It can be seen that the slip amount in the shallower part ranges from 0.6 m to 0.18 m.
5.2.

Strong ground motion simulation

Two different methods were applied in order to simulate the strong ground motion of the Puerto Plata
earthquake: Empirical attenuation-distance curves for PGA (peak ground acceleration) and PGV (peakground velocity): Here the required source information is solely limited to the seismic magnitudeand fault area, as types of outer faults parameters.
 Stochastic Green’s Function Method (SGFM): This method simulates the strong ground motionof a large seismic event using synthetic ground motion records of a relatively small event thatshares the same source location with that large one in place of the real records in the EmpiricalGreen’s Function Method (Irikura and Miyake 2011-2012).
5.2.1 Empirical attenuation relations
The obtained PGA and PGV values, with both methods, when converted to Modified Mercalli
Intensity (I MM ), using the relationship between PGA, PGV and I MM (Wald et al, 1999) are very close
to the values shown in the isoseismal map (Table 3). The calculation was done for latitude andlongitude range: (18.0N 20.0N and 68.5W 71.5W). Also, we conducted the same strong motion
simulation method to the Latitude and Longitude of the cities of Puerto Plata, Santiago de losCaballeros and the capital city of Santo Domingo using the same parameters mentioned above, andthe obtained PGA and PGV values when converted to I MM are also very close to the values registeredin the isoseismal map.
Table 3. PGA, PGV and I MM values obtained with the empirical-attenuation curve.

Note: 1 Fukushima and Tanaka (1990), 2 Si and Midorikawa (1999), 3 Average of theModified Mercalli Intensity calculated from an average of both PGA values using Wald et al.
(1999) and the observed Modified Mercalli Intensity.
5.2.2 Stochastic Green’s Function Method
Stochastic Green’s Function is calculated using the way of Ohnishi and Horike (2004), and used inplace of the real records, afterwards, in the Empirical Green’s Function Method.
The synthetic accelerogram and velocity seismograms calculated for the area of Puerto
Plata. The PGA and PGV averaged over two horizontal components are 295 cm/s2 and 15 cm/s. The
Fourier spectra of the ground acceleration and pseudo-velocity response are shown in Figure 6; thisfigure shows an intensity of 5.2 in the JMA scale, which converted to the Mercalli scale is 7.6, thisvalue is quite close to the intensity felt on the Puerto Plata area in the isoseismal map.
Since specific geological data of the ground in the Puerto Plata area was not available, avelocity structure for firm ground (AVS30=300m/s) was assumed, so that the results of the simulation
would be closer to real

Figure 6. Up: Crustal structure (Type J5 of CRUST2.0, top-left). Ground structureused for simulation (top-center). Down: Fourier amplitude spectrum of groundacceleration (bottom-left) and Pseudo velocity response (bottom-center) of the Puerto
Plata earthquake: Synthetic strong motion waveform data for acceleration (top-right)
and velocity.

6. CONCLUSSION
A rupture duration of 14 seconds and total seismic moment of M 0 = 1.446×1019 Nm (M w 6.7). This
seismic moment is higher than the one calculated by the NEIC-USGS (M 0 = 7.1×1018 Nm) (M w 6.5)
and the Global CMT (M 0 = 5.25×1018 Nm) (M w 6.4), taking around 10 iterations in order to obtain an
acceptably small value of ABIC. The slip amount is higher in the area that surrounds the hypocenter(seismogenic zone) with a maximum slip of 0.30 m. This slip distribution along strike matches thefault size estimated on the aftershocks distribution.
With the empirical attenuation curve we can see that the area of the epicenter where thePGA values are within 350 and 400 for both methods and the PGV between 20 and 25; these values
converted to I MM are very close to the I MM value near the epicenter in the isoseismal map which is in
between 8 and 9. Despite the Stochastic Green’s Function Method preciseness, this simulation
outcome is not very reliable due to the lack of important data.

7. REFERENCES
Dolan, J. F., and Bowman, D. D., 2004, Seismological Research Letters, 75, 587-591.
Fukushima, Y., and Tanaka, T., 1990, Bull. Seism. Soc. Am. 80, No. 4, 757-783.
Irikura, K., and Miyake, H., 2011-2012, Lecture notes on Strong Motion Seismology, IISEE/BRI.
Kuge, K., 2010-2011, Earthquake Source Process, IISEE/BRI Lecture note.
Si, H., and Midorikawa, S., 1999, J. Struct. Constr. Eng. AIJ, 523, 63–70 (in Japanese withEnglish Abstract).
Website: Global CMT catalog search http://www.globalcmt.org/CMTsearch.html
Website: U. S. Geological Service (USGS). (http://neic.usgs.gov/neis/eq_depot/2003/eq_030922/).
Yagi, Y., and Fukahata, Y., 2011, Geophysics J. Int. 186, 711-720.