List of the earthqaukes studied. Origin time (UT), latitude, longitude, and depth are taken from PDE catalog. GCMT parameters are from the GCMT catalog. W Phase results indicate frequency band, number of stations, and channels used for the final solutions. The solutions are the scenario of using PDE location and grid searching For td. The derived tdandMw are shown.
Ocean-wide tsunamis are almost always triggered by mega-earthquakes, rupturing substantial segments on the zones of interplate boundarybetween subducting slab and overriding plate. Whether, or specifically when, a mega-earthquake will occur on a subduction zone is crucial for preparedness of tsunami hazard mitigations. However, conditions controlling occurrences of mega-earthquakes are not well understood, as the 2004 Mw9.3 Sumatra earthquake contradicted the previous paradigm that the maximum-size earthquake in a given subduction zone can be predicted from its simple tectonic parameters such as slab age and convergence rates (Ruff & Kanamori, 1980). On top of that, the 100 years or so sampling periods of earthquake catalogue since instrumental era are not long enough to reveal the repeating periods, if any, of mega-earthquakes. Consequently, the fault dimensions necessary for a mega-earthquake to rupture in general become the most reliable constraint on its potential - for any subduction zone with length greater than 1,000 km, one cannot rule out the possibility that it is currently accumulating the strain energy of a mega-earthquake, which will release eventually to trigger ocean-wide tsunamis.
In the South China Sea (SCS) region, the Manila subduction zone, stretching over 1,000 km northward from Mindoro, Philippines to offshore SW Taiwan, poses a potential zone for the occurrence of a mega-earthquake and subsequently, thethreat of widespread SCS tsunami hazards. Indeed, it is a reasonable doubt that the historical tsunami hazards on coasts of SW Taiwan in 1661 and 1782(Soloviev & Go, 1984) might be caused by SCS tsunamis, amplified by SCS shelf when approaching Taiwan straits. Although it remains an open question as to when (or even whether) a mega-earthquake will occur in the Manila subduction zone, we endeavor to establish a warning system in Taiwan for SCS tsunamis.
Being able to predict the arrival times and amplitudes of the approaching tsunamis promptly after the occurrence of a tsunamigenic earthquake is key to the success of a tsunami warning system, which in turn hinges on modeling the evolution of earthquake-generated tsunami waves. There are three distinctive stages involved: generation, propagation, and run-up (Titov, 1997). For the generation stage, rapid determination of reliable earthquake source parameters (lon., lat., depth, moment, and focal mechanism) is crucial to infer the seafloor vertical displacements as initial conditions for tsunami wave propagations. Here, we adopt the source inversion using seismic W phase, which carries the long period information of earthquake source at high speeds and thus is suitable for tsunami warning purposes (Kanamori & Rivera, 2008). We propose using data of “extended BATS” for source inversion of SCS earthquakes. “Extended BATS” refers to the ongoing expansion of the Broadband Array in Taiwan for Seismology (BATS; Kao et al., 1998) network, leaded by the Institute of Earth Sciences Academia Sinica (IESAS) and incorporated with the Vietnamese Academy of Science and Technology (VAST) and the Philippines Institute of Volcanology and Seismology (PHIVOLCS) to construct a real-time network that encompasses SCS [see Figure 10 of Huang et al. (2009)]. In this study, we test the applicability of W phase inversion for past large SCS earthquakes, using data of stations that mimic distributions of “extended BATS”– the bulk of BATS network plus eight other stations surrounding SCS (Figure 1).
Given seismic source parameters, the vertical static seafloor displacements, as initial conditions of tsunami wave propagations, can be determined immediately by solving the surface deformations of an inclined fault in a half-space elastic medium, such as those described by Mansinha & Smylie (1971) and Okada (1985). However, the propagation stage that followed is the most time consuming stage and for the warning purposes of a regional scale, it is hardly practical for real-time calculations. Here, we employ the unit tsunami methods (Lee et al., 2005) to solve the dilemma. The methods divide all potential source regions into pixels and calculate the propagations of an initial unit amplitude uniformly assigned on the pixel. The resulting waves to a target station are stored in database and referred to unit tsunamis of that station-source pair. In the case of a real event, the predicted tsunamis of one station are synthesized by linear combinations of those unit tsunamis that are paired the station with all source pixels. The weightings of unit tsunamis in linear combination are determined by the vertical seafloor displacements averaged over the region of the corresponding pixel. We will elaborate the methods further in next section when building the database. Unit tsunamis are conceptually analogous to Green’s functions in Seismology and are applicable solely for linear systems. Therefore, in our warning system, we strategically focus on predicting the offshore arrival times and amplitudes of approaching tsunamis that is prior to the run-up stage and subsequently, the non-linear effects are minor and can be ignored.
2. Data and methods
2.1. Wphase inversion
We refer readers to Kanamori & Rivera (2008) regarding theory, modeling, and source inversion of W phase. We recapitulate here that W phase can be interpreted as superposition of the fundamental mode, first, second and third overtones of spheroidal modes at long period and can be synthesized by normal-node summation. The Green’s functions of six moment tensor elements are pre-computed for a distance range of 0°≦Δ≦90° with an interval of 0.1° and for a depth range of 0-760 km. The synthetic waveforms of an earthquake are derived by convolving the Green’s function with its moment rate function, which is a triangular function defined by two parameters, half duration, th, and the centroid delay, td(Figure 2). The broadband seismic data are deconvolved to displacement with instrument response removed and a band-pass filtered. A time period of 15Δs (Δ epicentral distance in degree) from the beginning of P wave is windowed to extract W phase. A time domain recursive method is used not only for real time operation but also for using available data to the point where it gets clipped at the large amplitude S or surface waves. The same procedures are applied on synthetic waveform and, together with data, a linear inversion is performed on concatenated time series using a given hypocenter location and origin time Kanamori & Rivera (2008).
We sorted out earthquakes from the GCMT catalog (Dziewonskiet al., 1981; Ekströmet al., 2005)that occurred between January 2000 and July 2009, were bounded by 10°N/26°N and 115°E/135°E, and had moments greater than 1025 dyn-cm (Figure1). Table 1 listed the PDE (Preliminary Determination of Epicenters) locations, td, th, and Mw values of the GCMT solutions for the 38 selected earthquakes. We divided the earthquakes into near and far groups relative to central Taiwan (Figure1). A virtual regional array, including the bulk of current BATS and a few stations from F-NET, Malaysian National Seismic Network, etc., was set up to mimic the distributions of extended BATS (Figure 1). The three component (ZNE) LH channel data (1 sample-per-second) were collected through Web site for W phase source inversions. Some stations were eliminated with their data continually exhibiting low quality (e.g., LYUB). We tested for six scenarios of two groups: group I using vertical component only and group II using all three ZNE components. In each group, three scenarios are tried with different level of knowledge on earthquake parameters: (1) using the GCMT centroid parameters (impractical in real events), (2) using the hypocenter parameters (lon., lat., depth, origin time) reported by the PDE catalogue, and a centroid time, td, determined by grid search; the source half duration, th, was set equal to td, (3) the same as (2) except the centroid location (lon., lat.) determined by a 2-dimensional grid search (depth fixed at PDE’s). We indicated the three scenarios as gCMT location, td location, and [td+xy] location, respectively. The cutting distance – within which data are removed – was determined to be 1.5° by examining results as a function of cutting distance. The inversion was iterated three times with an increasingly assigned threshold to discard stations of bad fitting data from previous iteration. An example of final fitting between synthetic and observed waveforms of vertical component for the first shock of the Mw7.0 Dec. 26, 2006, Pingtung earthquake is shown (Figure 3). The frequency band used to filter seismic waveforms basically followed Table 1 of Hayeset al. (2009) except that for a few Mw∼6.0 earthquakes, the band was fine-tuned to have better results.
2.2. Unit tsunami methods
The unit tsunami methods divide all potential source regions into pixels and assign an initial unit amplitude uniformly on each pixel, called unit source (Figure 4). The propagations of unit sources are pre-calculated with the resulting tsunamis at observing stations recorded as unit tsunamis to be stored in database. Each and every unit source has his unique unit tsunamis to a specific station stored and is pulled out to synthesize the composite tsunamis of that station in the case of a real event. The synthetics are done by linear combinationwith weighting factors determined by vertical seafloor displacements averaged over the region of corresponding unit source. In this way, the time consuming propagation stage is performed in advance, which makes prompt issuing of tsunami warnings more efficient. The principle of unit tsunami methods is conceptually analogous to Green’s function in Seismology and only works in linear systems. As compliance, we aim at predicting the offshore – prior to the run-up stage - arrival times and amplitudes of approaching tsunami waves. Under the circumstances, tsunami amplitudes are in general at least an order of magnitude less than the ocean depths and the nonlinear effects, such as convections and bottom frictions, are negligible. Furthermore, at open seas, the wavelengths (tens to hundreds of kilometers) of ocean-wide tsunamis are always much greater than ocean depths (a few kilometers). The vertical components of accelerations are thus neglected and horizontal motions of water mass are taken as uniform from top to bottom (Satake, 2002). Finally, we end up with the simplest form of hydrodynamic equations to govern the motions of incompressible fluids, i.e. the linear shallow water wave equations.
The potential source region of the Manila subduction zone was divided into 14×10 square pixels, each with 0.5° ×0.5° in size and an initial vertical seafloor displacement of 1 m was assigned to each pixel to constitute the group of unit sources (Figure 4). We employed Cornell Multigrid Coupled Tsunami Model (COMCOT; Liu et al., 1998) to simulate the propagations of each unit source. COMCOT is a finite difference scheme and in our case, we solved for linear shallow water wave equations in spherical coordinates with 1 minute (~1.8 km) and 1 sec in space and time, respectively. The boundary conditions were total reflection for ocean-land interfaces and radiation for map boundaries. We set up 32 virtual stations representing existing tidal stations of the Central Weather Bureau (CWB; Figure 5). The wave fields of one station from one unit source were referred to the unit tsunamiscorresponding to the station-unit source pair. For each unit source, we simulated for four-hour propagations with the resulting wave fields as a function of time at the 32 stations to be stored as 32 unit tsunamis in database. In the end, a total of 14×10×32 unit tsunamis were stored in database for synthetics of tsunami waves in the case of a real event. The stored unit tsunamis were readily available for arrival time predictions even without the occurrence of real events. We applied Short Time Average over Long Time Average (STA/LTA; Allen, 1982), a conventional scheme for picking seismic P and S phases, to automatically pick the arrival times of unit tsunamis with results stored in database for arrival time predictions.
In the case of a real tsunamigenic event, we calculated the vertical seafloor displacements caused by the earthquake (Okada, 1985), using source parameters inverted by fitting seismicW phase. The weighting factors of each unit source were derived by averaging the vertical displacements over area of the corresponding unit source. The tsunami waves at one station were then synthesized by summing the 14×10 unit tsunamis of the station with corresponding weighting factors.
3.1. W phase inversion
We judge the qualities of solutions based on the discrepancies between those of W phase inversion and those of GCMT solutions of the same event, in terms of both moment magnitudes (Mw) and focal mechanisms. The comparisons of Mw for six scenarios are presented in Figure 6 with group I (Z component only) in the first column and group II (ZNE components) in the second column. The numbers in each box (scenario) indicate the absolute means of magnitude differences and corresponding standard deviations. The absolute means of all six scenarios are less than 0.1 unit, validating the application of W phase inversion using data of regional network to determine source parameters of SCS earthquakes greater than Mw5.9. Among the three tried scenarios (gCMT location, td location, and [td+xy] location), only the last two are practicable in real-time W phase inversion, among which the group II td location is the best scenario (least discrepancies). The Kagan rotation angle refers to the solid angle rotating from one double couple to another (Kagan, 1991) and thus is a measurement of focal mechanism discrepancies. We present the Kagan angles - relative to the GCMT solutions of the same events - of the six scenarios in Figure 7 with mean and standard deviation indicated, using the same fashion as Figure 6. Again, the group II td location is the one with the minimum mean among all real-time practicable scenarios.
3.2. Unit tsunamis
Characteristics of tsunami wave propagations in SCS and around Taiwan depend on bathymetric features and can be learned from simulated propagations of unit sources. Figure 8 shows the propagating waves at different time frames of an exemplary unit source, from which we learned that tsunamis triggered by earthquake offshore north Luzon will hit south Taiwan at about 20 min., that the propagation speeds along east Taiwan coasts are much faster than those along west Taiwan coasts, meaning more warning time for the west coasts, and that almost all coasts around Taiwan will be attacked by SCS tsunamis within three hours after the generation of tsunamis. Figure 9a shows the resulting 32 unit tsunamis of the exemplary unit source, which demonstrate that stations 14 to 26 are the most affected ones with southern tip of Taiwan (stations 20, 21, and 23) being the most vulnerable (Figure 6). The STA/LTA scheme works well in determining the arrival times of unit tsunamis (Figure 9b) and we compile the data to produce one arrival-time map for each unit source (Figure 10), which will also be stored in database and are readily for prompt arrival-time predictions.
4.1. W phase inversion
Although the W phase source inversion was initially designed using data of global network and is now routinely operated in this fashion by the U.S. Geological Survey (USGS) National Earthquake Information Center (NEIC; Hayes et al., 2009), the W phase inversion is also applicable using exclusively regional data (e.g., the Full Range Seismograph Network of Japan) provided there is dense broadband network at regional distances (epicentral distances < 12°)[Kanamori H., personal comm.]. In this study, we have attested that, using data of regional network (the mimic network of extended BATS), the W phase is equally applicable for inverting source parameters of SCS earthquakes greater than Mw5.9. We expect the applicability can be even improved as a result of better azimuthal coverage once the construction of the extended BATS is finished.
The experiences learned by conducting the W inversion of past SCS earthquakes are valuable for future implementation of real-time warning and are addressed in the following. (1) The frequency band of W phase used for inversion as a function of earthquake size as suggested by Table 1 of Hayes et al. (2009) is robust and can be directly implemented. The only exception is that when the sizes of earthquake and, accordingly, the signal/noise ratios are relatively small (Mw∼6.0), under the circumstances, different frequency bands do matter and we need to manually try out the best scenario. Since earthquakes of such size rarely trigger significant tsunamis, we can just implement the empirical frequency band of Hayes et al. (2009) and focus on major events. (2) Having tried six scenarios in two groups, we find that grid search for epicenters of centroid location ([td+xy] location) does not necessarily decrease the discrepancies, compared to those of td location– probably due to the use of GCMT solutions as a reference. For the scenario of td location, the use of ZNE components (group II) only slightly reduces the discrepancies relative to those of Z component only, however, at the expense of processing time. Therefore, we can alternatively use group I td location as the scenario of real-time W phase inversion for a more efficient warning system. (3) The patterns of discrepancies exhibit similar trends for Mw and Kagan angles (Figures 7 and 8), presumably because both are derived from the same moment tensors – product of W phase inversion.
4.2. Unit tsunamis
We test with a scenario Mw9.0 earthquake in the Manila subduction zone to assess the applicability of predicting maximum amplitudes of tsunami waves using unit tsunami method. The vertical seafloor displacements caused by the scenario earthquake are calculated using Okada (1985)’s scheme [Figure 11(a)]. We first follow the conventional approach to simulate the propagation of the displacements, which takes time. Secondly, we derive the weighting factors of each unit source by averaging over the vertical displacements within each unit source [Figure 11(b)]. The tsunami waves are then synthesized by a linear combination of unit tsunamis – pulled from database – with corresponding weighting factors, which can be done instantaneously. Figure 12 shows resulting tsunami waves of the two approaches for comparison purpose. A maximum-amplitude map can also be created showing the magnitudes on the corresponding stations (Figure 13). By comparing (a) and (b) of Figures 12 and 13, we conclude that the unit tsunami method can produce waveforms and maximum amplitudes similar to those of conventional simulation method, however, only a much less time is needed.
Liu et al. (2009) have also proposed procedures to establish a tsunami early warning system in the SCS. The procedures also apply unit tsunami methods to pre-calculate tsunami waves of fault segments along the Manila subduction zone. The primary differences between the two systems are the ways in determining the weighting coefficients of unit tsunamis in synthetics - Liu et al. (2009)’s procedures invert the real-time observations of nearby ocean-bottom pressure sensors whereas we calculate from source parameters inverted by the W phase in this study. The two approaches actually complement each other in that the former is more accurate - direct observations - but rely on well-placed sensors and the latter has more observations – seismic stations – but indirect through source parameters.
As a final remark, we point out that the unit tsunami methods are not compromised on regions of land or land-ocean boundaries. For pixels on land exclusively, the calculated unit tsunamis are zero for all stations and thus have no contributions to the synthetic tsunamis despite their non-zero weighting factors. Likewise, for pixels on coast lines, the calculated unit tsunamis only constitute contributions from the ocean parts, thus automaticallyproportional to the ocean-land ratio within the pixel. As a result, they will contribute proportionally to the synthetic tsunamis in the processes of linear combinations. As for regions of complicated bathymetries such as trenches, the vertical seafloor displacements in the regions tend to deviate from those estimated by Okada (1985)’s schemes, which assume a half-space elastic medium. However, both conventional approach and unit tsunami methods will suffer from the same deviations. The unit tsunami methods are equally applicable to warnings of local, regional, and global scales. However, strategies in building database of unit tsunamis should be adjusted accordingly; in global scale with greater potential source regions, we may use greater pixels to reduce the size of database, at the expense of reducing spatial resolutions of seafloor displacements; in local scale, the pixels cannot be too small in order that the conditions of long wave validate.
In this study, we propose to combine W phase inversion and unit tsunami methods to build a tsunami warning system in Taiwan for earthquakes in the SCS region. The W phase inversion allows us to rapidly determine moment tensors of large earthquakes for the calculations of vertical seafloor displacements. The applicability of W phase inversion for SCS earthquakes using BATS stations and its extension has been attested and expects to be improved, pending future completion of extended BATS. We have built a database of unit tsunamis for the source region of the Manila subduction zone and the prediction of arrival times is readily available once the epicenter of tsunamigenic earthquake is known. Upon assessment of a Mw9.0 scenario earthquake, the prediction of maximum amplitudes on various coasts in Taiwan using unit tsunami method is comparable with that of conventional simulation, but with a much less time needed.
This work was supported by Central Weather Bureau with grant number MOTC-CWB-101-E-09 and National Science Council (NSC) in Taiwan with grant number 99-2116-M-008-040.