Seismic Hazard Assessment for the Taiwan Region
on the Basis of Recent Strong-Motion Data and
Prognostic Zonation of Future Earthquakes
VLADIMIR SOKOLOV1, ARKADY OVCHARENKO2, CHIN-HSIUNG LOH3and KUO-LIANG WEN4
1Geophysical Institute of Karlsruhe University, Hertzstr. 16, 76187, Karlsruhe, Germany;2Institute
of Geophysics, Urals Division, Russian Academy of Science, Amundsen str. 100, 620016,
Ekaterinburg, Russia;3National Science and Technology Center for Hazard Reduction, 200, Sec. 3, Hsinhai Rd., Taipei, Taiwan, R.O.C.;4Institute of Applied Geology, National Central University, Chung-li, Taiwan, R.O.C.
(Received: 13 January 2003; accepted: 11 December 2003)
Abstract. The paper describes an integrated approach to seismic hazard assessment, which was applied to the Taiwan region. First, empirical models for ground motion estimation in the region were obtained on the basis of records from recent (1993-1999) earthquakes. The database includes strong-motion data collected during the recent Chi-Chi earthquake (M= 7.6, 21 September 1999) and large (M = 6.8) aftershocks. The ground-motion database was also used for evaluation of generalised site amplification functions for typical soil classes (B, C and D). Second, the theoretical seismic catalogue (2001–2050) for the Taiwan region had been calculated using the 4D-model (loc-ation, depth, time) for dynamic deformation of the Earth’ crust and 5D-model (loc(loc-ation, depth, time, magnitude) for seismic process. The models were developed on the basis of available geophysical and geodynamic data that include regional seismic catalogue. Third, the region and site and time-dependent seismic analysis, which is based on schemes of probable earthquake zones evaluated from the theoretical catalogue, regional ground motion models, and local site response characteristics, has been performed. The seismic hazard maps are compiled in terms of Peak Ground Acceleration (PGA) and Response Spectra (RS) amplitudes. The maps show distribution of amplitudes that will not be exceeded with certain probability in condition of typical soil classes during all possible earthquakes that may occur in the region during time period of 2003–2025. The approach allows introduction of a new parameter that describes the dependency of seismic hazard on time, the so-called “period of maximum hazard”. The parameter shows the period, during which every considered site will be subjected by the maximum value of ground motion characteristic (PGA or RS).
Key words: future earthquake zonation, ground-motion models, integrated approach, time-dependent seismic hazard
1. Introduction
During the last one-two decades large earthquakes caused massive loss of human life and extensive physical destruction throughout the world (Iran, 1990; Japan,
1995, Turkey, 1999; India, 2001, etc.). At the present state of knowledge, the proper level earthquake-resistant design seems to be one of the most effective ways to prevent grave consequences of the earthquake. The goal of design is to produce a structure or facility that can withstand a certain level of shaking without excessive damage. That level of shaking is described by a design ground motion (design earthquake), which can be characterised by engineering ground motion parameters (peak amplitudes, Fourier and response spectra, seismic intensity, ac-celeration time histories, etc.). The specification of engineering (or design) ground motion parameters is the goal of Seismic Hazard Analysis (SHA). It involves the quantitative estimation of ground-shaking hazard at a particular site taking into account characteristics of potentially dangerous earthquakes around the site.
There are two types of uncertainty in seismic hazard analysis, namely: (I) the uncertainty connected with effect of a future earthquake, or uncertainty in ground motion parameter evaluation for a particular earthquake, region and site condition; (II) the uncertainty due to the unpredictable nature of a future earthquake (when, where, how large). The first uncertainty is caused by incomplete knowledge and may be reduced by using additional data. For example, the Chi-Chi earthquake (ML = 7.3, September 21, 1999, Taiwan) and its aftershocks produced a rich set
of more than 1000 ground motion recordings, which are currently under extensive study (Shin and Teng, 2001).
At present there is a common opinion (e.g., Rundle et al., 2000) that the only way to reduce the second type of uncertainty is the development of numerical models, which approach the reality. In most cases the definition of the earthquake source zone parameters (configuration, maximum possible magnitudes, and earth-quake recurrence) is carried out on the basis of the historical data, i.e., using events which already occurred. In other words, statistical properties of the known seismi-city are assumed to be stationary in time and space, and therefore should apply in the foreseeable future. There are many examples, including the case of the above-mentioned Chi-Chi earthquake, when the erroneous estimation of these parameters led to grave consequences.
In this paper we describe an integrated approach to seismic hazard assess-ment, which was recently (1999–2002) applied to the Taiwan region. The region is characterised by high seismicity, (Figure 1) and earthquakes occur up to about 200–300 km. Taiwan has repeatedly been hit by large earthquakes and at least 19 damaging earthquakes occurred there during the twentieth century (Shin and Teng, 2001). The worst of the earthquakes, the Chi-Chi earthquake, caused heavy damage in the central and western part of Taiwan. The official casualty figure included 2470 dead, more than 11,000 injured, and more than 100,000 structures destroyed (Shin and Teng, 2001). The most recent, so-called 331 Hualien earthquake (ML = 6.8, 31
March 2002) occurred approximately 100 km to southeast from the city of Taipei. The earthquake caused collapsing of residential building in the city and falling of construction facilities from the top (56 stories) of Taipei Financial Centre that is currently under construction. Five construction workers were killed.
Figure 1. Seismicity of Taiwan. (a) A 3D plot showing hypocentres of earthquakes occurred
during 1973–2000 (M > 3). (b) A scheme of epicentres for events of M > 3 (1990 to 1999). The 19 damaging earthquakes in Taiwan during the 20th century are shown as open stars. The figure is taken from Shin and Teng (2002), and two recent events were added, namely: the Chi-Chi earthquake (21 September, 1999) and the Hualien earthquake (31 March, 2002).
When performing assessment of seismic hazard for the region, we tried to re-duce both the uncertainties mentioned above. The parameters of future earthquakes are evaluated on the basis of dynamic 4D (x, y, z, t – location, depth, time of occurrence) model for the process of Earth crust deformation and 5D (x, y, z, t, M) – location, depth, time of occurrence, magnitude) statistical model of seismicity. The final purpose of the model application is to compile a “theoretical seismic catalogue” of future earthquakes, which could be used as a basis for evaluation of the future seismic zone configuration and maximum possible magnitude. The information provides necessary information for seismic hazard assessment in terms of design input ground motion parameters, which is based on Fourier Amplitude spectra (FAS) of ground acceleration. The approach allows consideration of re-gional (source scaling, attenuation relation) and local (site condition) peculiarities of seismic ground motion. It is possible to obtain the results in terms of seismic in-tensity, peak acceleration, response spectra and characteristic accelerograms using a common ground motion model (Sokolov, 2000a,b). The features of input ground motion models, both the regional and local ones, were analysed and quantitatively described for the Taiwan region on the basis of recent earthquakes including the large MW = 7.6, 1999 Chi-Chi earthquake.
The schemes of seismic source zones (future earthquakes), regional spectral models, as well as generalised site response characteristics, are used, as input parameters, for evaluation of ground motion characteristics for soil classes B, C and D throughout Taiwan Island (Taiwan zonation maps). The maps show Peak Ground Acceleration and Response Spectra that will not be exceeded with a certain probability (50%, 84% and 97%) during time period of 2003–2025.
2. Ground Motion Models
Empirical FAS models for estimating design input ground motion parameters in the Taiwan region have been obtained recently (Sokolov et al., 1999, 2000a) on the basis of the recordings of small and moderate earthquakes in 1993–1998 (5.0 < ML <6.5, 1380 accelerograms). Two models were developed. The first model is
represented by the empirical spectra of ground acceleration estimated for reference distance R = 1 km (Figure 2a). The records were obtained on various types of soil: from rock sites to soft soils of different thickness. Therefore, averaging the variety of site conditions, the spectral model corresponds to the “average soil”. Answering the question what does the “average soil” really mean, we assume that the “generic soil” class used by Boore and Joyner (1997) seems to be the most appropriate analogue by means of amplification amplitude. On the other hand, the “average soil” class, reflecting response of different typical soil classes, should be characterised by a broadband amplification.
To extend the “average soil” spectral model to the larger magnitudes, the log-arithmic increments (γM)of spectral amplitudes (A) per unit of magnitude were
calculated for frequencies 0.2–12 Hz in the following form
γM(f )= log A(f )/M, (1)
where M is the increment of magnitude between the above mentioned mag-nitude ranges, and A(f ) is the reference distance spectral amplitude (mean values, Figure 2a). The approach was suggested by Khalturin et al. (1975) and further developed by Chernov (1984, 1989). Figure 2b shows the values γM, which were
determined using the empirical “average soil” spectra for the Taiwan region. The distribution of γM values versus frequency shows that the values of γM increase
with decreasing frequency and the coefficients γM are nearly constant in the
low-frequency range. In other words, as the source energy increases, the intensity of low-frequency vibrations increases more rapidly than that of high-frequency vibra-tions. The fact agrees with the results of the previous studies of the γMcoefficients
(Chernov, 1989; Sokolov, 1998a) on the basis of empirical data obtained in various seismic regions (California, northern Italy, Caucasus, and Central Asia).
The spectrum at a given distance R is calculated by multiplication of the “reference distance” spectra from a proper magnitude range by
Figure 2. “Average soil” spectral model for various intervals of magnitude M. (a) Empirical
spectra (solid lines) calculated for reference distance R = 1 km. Numbers near the curves denote standard deviation and number of the used records (in parentheses). Dashed lines show the spectra estimated using the coefficients γM (see text). (b) Frequency dependence of the coefficients γM (logarithmic increment of spectral amplitudes) calculated using different
magnitude ranges (1) 5.5 < M < 6.0 versus 5.0 < M < 5.5; (2) 6.0 < M < 6.5 versus 5.5 < M < 6.0; (3) averaged values).
where b= −1.0 for distances R < 50 km and b = 0.0 for 50 < R < 150–170 km; VS is the shear wave velocity in the Earth’s crust (3.6 km/s); Q(f ) = 225f1.1
for deep (depth more 35 km) earthquakes, Q(f ) = 125f0.8 for shallow (depth between 10–15 and 30-35 km) and Q(f ) = 80f0.8 for very shallow (depth less 10–15) earthquakes. It was shown that the stochastic calculation technique based on the spectral model allows for a satisfactory prediction of the peak ground ac-celeration and response spectra for distances up to about 200 km in the region. Comparison of the model with ground motion data obtained during the recent Chi-Chi earthquake and large aftershocks (Sokolov et al., 2002a, b) allowed us to conclude that the local magnitude-based “average soil” spectral model, which was recently developed for the Taiwan region, may be applied for prediction of spectral amplitudes (within frequency range 0.3–12 Hz) and peak ground acceleration for larger reverse-fault earthquakes in the region.
The second model, so-called “Very Hard Rock” or VHR spectra, was obtained as follows. The general model for the Fourier acceleration spectrum A at frequency f is given by
A(f )= (2πf )2CS(f )D(R, f )I (f ), (3)
where C is the scaling factor; S(f ) is the source spectrum; D(R, f ) is the attenu-ation function, and I (f ) represents frequency-dependent site response. The scaling factor is calculated as
C = ([Rθ φ]F V )/(4πρβ3Rb), (4)
where Rθ φis the radiation coefficient, F is the free surface amplification, V
repres-ents the partitions of the vector into horizontal componrepres-ents, ρ and β are the density and shear velocity in the source region, R is the source-site distance, and b = 1.0 for distances R < 50 km and 0.0 for 50 < R < 150–170 km. A commonly used source function S(f ) in the Brune (1973) single-corner-frequency model is
S(f ) = M0/[1 + (f/f0)2]. (5)
For the model, the source acceleration spectrum at low frequencies increases as f2 and approaches a value determined by f0 (corner frequency) and M0 at frequencies f > f0. The value of f0 can be found from the relation f0 = 4.9× 106β(σ/M
0)1/3. Here σ is the stress parameter in bars, M0is the seismic moment in dyne-cm and β in km/s.
The function D(R, f ) accounts for frequency-dependent attenuation that mod-ifies the spectral shape. It depends on the hypocentral distance R, regional crustal material properties, and the frequency-dependent regional quality factor Q that represents anelastic attenuation. These effects are represented by
where P (f, fmax)is the high-cut filter.
In the Taiwan region for hypothetical very hard rock (VHR) sites (ρ = 2.8 g/cm3, β = 3.8 km/s, I (f ) = 1), it has been proposed (Sokolov et al., 2000a) to use the single-corner frequency Brune source model (Equations (3)–(6)). High-frequency amplitudes are reduced, through the kappa operator (Anderson and Hough, 1984), by multiplying the spectrum by the factor (f ) = exp(−πκf ) (κ = 0.03). The Q-model is described above. The seismic moment M0 and stress parameter σ are determined using recently proposed regional relationships between M0and magnitude (ML)(Li and Chiu, 1989)
log M0= (19.043 ± 0.533) + (0.914 ± 0.035)ML (7)
and between σ and M0(Tsai, 1997)
log σ = −3.3976 + 0.2292 log M0(±0.6177). (8)
Figure 3 shows the VHR spectra for various magnitudes. When analysing strong motion data obtained during the Chi-Chi earthquake and the largest aftershocks (ML = 6.8), it has been found (Sokolov et al., 2003a, b) that the high- and
in-termediate frequency (more than 0.3–0.4 Hz) radiation during the earthquakes are generated mainly from the fracture of small-scale fault heterogeneities. Therefore, the attenuation models for ground motion acceleration, which are based on long-period seismic moment values (moment magnitude) and were developed for other type of faulting, fail to predict the peak ground motions for the mainshock (EERI, 1999; Tsai and Huang, 2000; Boore, 2001). It seems to be more preferable, at least for the case of shallow thrust earthquakes, to use so-called “short-period seismic moment” proposed by Koyama and Zheng (1985), which reflects the properties of heterogeneities. On the other hand, the local magnitude ML determined using
Wood–Anderson seismograph with dominant period 0.8 sec may be considered as a magnitude related to the “short-period seismic moment”. Therefore the VHR model, which was developed using records from small and intermediate earth-quakes (ML < 6.5) and which is based on the regional relationship between M0
and local magnitude ML, may be applied for magnitudes up to ML= 7.3–7.5.
3. Site Amplification Characteristics for Typical Soil Classes in the Taiwan Region
The most common procedure (Borcherdt, 1970) for estimation of the site response from earthquake data is to determine the sediment-to-bedrock ratio by dividing the Fourier spectrum of a site by that of a nearby reference (rock) site using earthquake records. However, in many cases it is not possible to find a station, which is suitable as a “hard-rock reference” station. It was proposed recently (Sokolov, 1998b) to use so-called “very hard rock” (VHR) spectral model for estimation of site response characteristics in terms of frequency dependent amplification (spectral ratios). The
Figure 3. Very Hard Rock (VHR) theoretical ω-square spectra of ground acceleration for the
Taiwan region, distance R= 10 km.
approach consisted in calculating spectral ratios between spectra of actual earth-quake records (horizontal components) and those modelled for a hypothetical VHR site. On the one hand, besides local site response, the spectral ratios include effects of source rupture peculiarities and inhomogeneous propagation path. On the other hand, when using a large enough number of events varied by magnitude, source depth and azimuth, the effects of focal mechanism and directivity are expected to be averaged out. The approach was applied for evaluation of characteristics of spectral amplification for particular sites in Taipei basin (Sokolov et al., 2000b, 2001).
Empirical amplification functions for site classes B, C and D in the Taiwan Island were evaluated (Sokolov et al., 2003a, b) as the ratio between Fourier amp-litude spectra of recorded accelerograms and spectra modelled for hypothetical very hard rock (VHR). Site classes B, C, and D can be distinguished by geo-logical age and rock type. Site class B may include igneous and metamorphic rocks, limestone and hard volcanic deposits. Sandstones, shales, conglomerates of Miocene age or older are all classified as class B. Most class B sites are located at hilly and mountainous highlands, except for some on limestone, igneous and metamorphic rock sites at lower elevations. Pliocene and Pleistocene sandstones, shale/mudstones, and conglomerates, volcanic breccias and pyroclastic rock of a similar age are classified as class C. Most class C sites are to be found on lateric terraces and highlands. Late Quaternary deposits, other than lateric terraces, such as loose sand, silt, clay, and gravel deposits, are classified as class D (stiff soil) or E (soft soil). Soils in class D are fluvial terraces, stiff clays, and sandy gravel deposits. To differentiate between stiff soil and soft soil sites, SPT-N values and geomorphologic classifications are used (Lee et al., 2001).
For evaluation of the generalised site response functions we used recordings obtained during earthquakes in 1993–1998 (5.2 < M < 6.5, 21 events) and three large events of the Chi-Chi earthquake sequence, namely: mainshock (M = 7.6) and two aftershocks (M = 6.8). The distribution of recordings regarding the soil
classes is as follows: class B – 146 records; class C – 180 records; class D – 770 records. Comparison between mean-amplitude amplification functions for gener-alised site classes is shown in Figure 4a. Note that the curves for the site classes are similar for frequencies greater than 2–3 Hz. When comparing distribution of standard deviation versus frequency for various soil classes (Figure 4b), we can note the following. On the one hand, soil class B stations are characterised by the highest variability of site response at frequencies more than 2 Hz. On the other hand, several stations located on class C and D sites exhibit high variability of site response in low-frequency range (f < 0.8−1.0 Hz). Most likely, the phenomenon, which is especially expressed for the case of class D sites, may be explained by influence of surface waves (shallow earthquakes, alluvium-filled valleys).
Let us compare (Figure 4c) generalised mean-amplitude amplifications result-ing from the Taiwan data with that proposed for western US (Boore and Joyner, 1997) and Greece (Klimis et al., 1999). On the one hand, class C and class D amplification functions proposed for western US, as well as the Hellenic class C amplification function, exhibit lower values than correspondent Taiwan functions for the whole frequency range considered. The Hellenic class D amplification is characterised by the higher amplitudes than the Taiwan class D amplification at frequencies more 3-5 Hz. Extremely high weathering rate in Taiwan due to the wet climate, which leads to high attenuation of high-frequency radiation, may be considered as the reason of the difference.
The amplification functions demonstrated their reliability when compared with independent data obtained by Lee et al. (2001). It is also necessary to note the prominent influence of geologic and geomorphologic factors on site amplification function for the considered cases of rock (class B), soft rock or very dense soil (class C), and stiff soil (class D) sites. The influence is reflected by large variations of amplitudes and predominant frequencies between particular stations with the same site class that was assigned on the basis of the rock age and geological clas-sification. Therefore, as it has been noted by Lee et al. (2001), further studies on site classification should be carried out using more specific subsurface geotechnical data including thickness of soil, shear wave velocity, and density.
4. Future Earthquake Zonation 4.1. BASIC PRINCIPLES
The slow deformation of the Earth crust is the main source of the crust fracture and fast energy release in the form of earthquakes. Geodetic and strain-metric ob-servations often reveal non-periodic slow undulations in the Earth crust. Numerical examples of instrumental observation of the dynamic objects in the form of charac-teristic non-linear trend of displacement may be found in GPS (Global Positioning System) data or in PSWL (Persistent Sea Water Level) monitoring (Figure 5a– c). Direct instrumental observations along very long (VLBI) and short (2-color EDM) baselines also track the wave-like processes of various periodicity. Precise
Figure 4. Generalised site amplification functions B, C and D. (a) Mean amplification
func-tions for site classes C, B and D, Taiwan region. (b) Distribution of standard deviation of the amplification values versus frequency. (c) Comparison of amplification functions evaluated in this study (curves #1) with generic site amplification functions proposed for western America by Boore and Joyner (1997) (curves #2, soil classes C and D) and with amplification functions proposed for Greece by Klimis et al. (1999) (curves #3, soil classes C and D).
measurements performed by dense networks allow concluding that (a) deformation of the crust may be, in first approximation, considered as having elastic nature and (b) deformation of the crust is complex and dynamic process, which includes wave-like components (Figure 5d).
When constructing the models of the complicated geodynamic processes inside the Earth crust (deformation and seismicity), it is necessary to use analytical ap-proximation. Two such approximations (4D and 5D-models) and the technique for the models adjustment to the observed data have been recently proposed
(Ovchar-Figure 5. Dynamic nature of the Earth crust deformation. (a) typical wave-like deformation
that could be revealed from the GPS data (station USNA, 38.98N, 76.48W). (1) Observed latit-udinal displacements; (2) curve, which obtained after smoothing (14-points running average), reveals five fast-propagating waves of (probably) seasonal nature; (3) third-power polyno-mial shows a slow-propagating wave of deformation. (b) Examples of synchronous changes (longitudinal displacement) for two closely located stations USUD (36.133N, 138.362E) and TSKB (36.105N, 140.087E). (c) Example of slow-propagating wave of deformation that was revealed after filtration of the PSWL (sea water level) data (zero lines were changed). (d) Modelled example of the arbitrary separate waves of deformation (curves (1)–(4) and the summarised (observed) deformation (curve 5). Period of GPS observation is bounded within the lines; arrows show duration of various types of observation.
enko 1997a,b, 1988, 1999; Ovcharenko et al., 2000, 2001). The data base that was used for modelling contains all basic data of geophysical (geodynamic) observa-tion: seismic catalogues, data obtained by global and local GPS networks, very long base line interferometery (VLBI) data, sea water level observations, high precision levelling, etc.
Historically, the first technique of instrumental observation, which reflects the dynamic process of the Earth crust deformation, is sea water level measurements conducted by coastal network (e.g., Tolkatchev 1996). In the course of time, the large scientific projects, which were devoted to study of recent geodynamics of continents, were performed. Among then we should note repeated levelling of large areas; tracking of local deformation using special equipment (deformographs); tracking of surface inclination using special network; tracking of variations of distance between two points located at the surface (VLBI, 2-color geodimeter, etc.); permanent and repeated tracking of location of reference points (GPS-stations, global and local network); evaluation of direction and amplitude of main stresses for earthquake source using seismic recordings. These kinds of observa-tions may be called direct tracking of various characteristics of the process of crust deformation.
Tracking of sea water level, observation of variations of gravitational, mag-netic and electromagmag-netic fields, measurements of elastic properties of geological medium and water level in the boreholes, changes of the coastlines, volcanic activ-ity, etc, may be considered as the indirect data. The indirect data are frequently characterised by the cheapest, longest and high-quality time series (Figure 5d).
The general approach for development of the 4D (x, y, z, t – geographical location, depth and time) models of geological and geophysical processes is as follows. A function F[p, x, y, z, t] is introduced on the basis of available data. The vector p describes parameters of the process. Most of these parameters, as well as the type of the function F, are not initially known. It is supposed that the function F reflects the considered geological phenomenon adequately. The problem consists in a search of the vector p on the basis of uncompleted observational data {Fi, i =
1, n}. When determining the function F, we assume, as a basic principle, that the Earth’ crust contains the dynamic objects, namely: the wave-like slow-propagating variations of tectonic deformation. Many authors, starting with Russell (1845) have studied the possibility of appearance of solitary waves in various media (liquid, solid and vapour). Many authors also suggested the presence of these waves of deformation (stress) inside the Earth crust. For example, it has been supposed that waves determine the response of the lithosphere to changes of horizontal stresses (Elsasser, 1969), or that waves appear due to the variations of the Earth’s rotation (Guberman, 1987).
In the 4D-model, for a formalised description of complicated non-linear phe-nomena of the crust deformation, we introduce a special dynamic object, which is further called “Front of Dynamic Deformation” (FDD) (Figure 5d). We interpret the term “FDD” as a lengthy line-stretched area of extreme or anomalous values
of deformation, which propagate inside the crust with velocities ranging from 0.01 km/year up to 1000 km/year. Direction of propagation of the FDD is parallel to direction of displacement and distribution of deformation within the FDD is characterised by alternating signs. In the axial area the absolute values of deform-ational displacement increase from the area toward the edges of the FDD. When reaching the local maximum, the displacement decreases in an asymptotic way. The distance between points on the asymptotic branches, where the amplitudes of displacement are equal to a half of the maximum values, is called “width of FDD”. In essence, the front of dynamic deformation is described in a similar way, as a single planar seismic wave. The major difference consists in much lower velocity of the FDD propagation, and much greater amplitudes of displacement and width.
The data, which were collected during various types of geodynamic and geo-physical observations, are used for compilation of generalized databases for the modelling, namely (1) catalogue of geodynamic indicators or “indicator points”, and (2) file for amplitude calibration. The “indicator point” is used for fixing of parameters (location and time moment) of certain local peculiarity of geodynamic phenomenon. The essence of peculiarity depends on the type of used data. In the case of horizontal displacement, the “indicator point” fixes the moment of change of the sign of measured value (crossing of zero). For the case of seismic data it is the seismic event; for vertical displacement and sea water level data it is the moment of local extreme.
The dynamic 5D-model for seismic process (x, y, z, t, M – geographical location, depth, time, magnitude) had been developed on the basis of seismic cata-logue only (Ovcharenko et al., 2001). Seismic process is described by a system of dynamic objects, so-called planar dynamic “fronts of seismicity”, which were eval-uated during specific iterative procedure on the basis of initial observed catalogue. Every seismic event (earthquake) is described by, at least, five different “fronts of seismicity”.
4.2. SHORT DESCRIPTION OF THE MODELS 4.2.1. 4D-model for the crust deformation
It has been shown empirically (Nusipov and Ovcharenko, 1997; Ovcharenko, 1997a,b, 1998, 1999; Ovcharenko et al., 2000, 2001) that elastic deformation of crust may be approximated, with high accuracy, by a collection of multi-parametric functions (4D-model), which represent a set of dynamic objects (fronts of dynamic deformation or FDD). The idea is supported by results of theoretical studies of deformation of elastic media. The basic formulation of the model is given in Appendix 1. The 4D-model describes process of dynamic volumetric deforma-tion of Earth crust, which formally contains sources of earthquakes in peculiar discrete points. In other words, it is accepted, as a main principle, that seismic events are considered a priori as singular points of continuous dynamic process of deformation.
Summarising mathematical description of the 4D-model (Appendix 1), the function of elastic deformational displacement for a set of plane FDDs in infinite medium may be determined as
E(x, y, z, t)= n i=1 mi 1 qi (9) qi = p2i + 2 pi = (y − aix− biz− cit− di) (10)
where (x, y, z, t) are the Cartesian coordinates and time; (a, b, c, d) are the kinematic parameters of the model, which describes velocity and azimuthal para-meters of propagation; 2 is the effective width of the FDD; mi are the amplitude
(energy) characteristics of the model; n is the number of FDDs.
The dynamic model for a set of the FDDs may be constructed using various kinds of observed data by solution of the inverse problem. The goal of solution is determination of parameters{ai, bi, ci, di, mi, i = 1, n} which satisfy the variety
of observed data. It is proposed to decompose the initial task into two parts (see also Appendix 1). The first part, or inverse kinematic problem, consists in determination of kinematic parameters of the model{ai, b1, c1, d1, i = 1, n}. All the data that describe the dynamic nature of the process of deformation may be used for solution of the problem. The second part, or inverse dynamic problem, consists in determ-ination of amplitude parameters of the model {mi, i − 1, n} using the data that
were obtained during the direct observation of the Earth’s crust deformation (GPS, geodetic measurements, etc). The corresponding problem can be solved using the Tikhonov’s smoothing function
Mα[E, K, m] = k j=1 Ej − n i=1 miKij 2 + α n i=1 (mi)2, (11)
where Ej(x, y, z, t) are the components of deformation observed at stations
{xj, yj, zj, tj}, j = 1, k; mi (i = 1, n) are the unknown parameters of the model;
αis the parameter of regularisation; K is the kernel function that depends on the observed parameter of deformation. The solution of the similar tasks in geophys-ics has been thoroughly investigated and widely applied (Tikhonov and Arsenin, 1979).
The resulting 4D-model allows us to extend the process in space and time and to evaluate the components of tensor of deformation. In other words, if the parameters of the model have been found on the basis of experimental (limited in space and time) data, it is possible to evaluate the tensor of dynamic deformation outside the period of observation in future and in past. Obviously, the accuracy of the evaluation depends on duration of observation, density of network and dimensions of area covered by monitoring.
The possible seismic events, on the one hand, could be revealed by analysis of distribution of deformation inside the Earth’s crust. Figure 6 shows the results of application of the preliminary 4D-model, which was based on the data that were available before the large Chi-Chi earthquake (central part of Taiwan Island, MW = 7.6, September 21, 1999; thrust faulting, see also Shin and Teng, 2001).
The data include GPS (e.g., Yu and Chen, 1994, Yu et al., 1997) and sea water level observations, and seismic catalogue (Central Weather Bureau, Taiwan). The development of a narrow zone of intensive compression (depth about 10–20 km) in the central part of Taiwan Island is clearly seen (Figure 6a). The location, shape and absolute values of deformation correspond to the parameters of the Chi-Chi earthquake, and period of the maximum values of deformation corresponds to 1997–2000 years. On the other hand, it is proposed, as working hypothesis, to consider seismic events as the peculiar points of the field of dynamic deformation – the interaction of four or more fronts of dynamic deformation. The location of the Chi-Chi earthquake source is clearly fixed by the correspondent zone of interaction of several revealed FDDs (Figure 6b).
4.2.2. 5D-Model for Seismic Process
The 4D-model, as a model of deformation that causes the rupture processes inside the crust, uses the phenomena of rupture (earthquakes) as a geodynamic indicator of deformation. Earthquake magnitude is a characteristic of energy radiated during the rupture. Besides the amplitude of deformation (slip), magnitude also relates with other parameters (fault area, physical properties of rock, etc.) that are not dir-ectly considered in the 4D-model. The 4D-model determines location of idealised point that corresponds to interaction of several FDDs. However, it is necessary to consider a certain volume of the crust that is affected by the FDDs. The volume is characterised by a complex distribution of deformation and stress and development of earthquake rupture depends on mechanical and chemical processes in the crust. Therefore, it is not possible to establish a simple correlation between earthquake energy (magnitude) and parameters of deformation determined by the 4D-model. Bearing in mind the dynamic nature of the crust deformation, it is obvious that the rupture may occur before or after, near or far off the idealised point of the FDDs interaction. These factors lead to the uncertainties in magnitude, time and location of potential seismic events predicted by the considered deformational model.
It is natural to assume that “background” seismicity reflects the above men-tioned parameters of the Earth’ crust in the given region. Therefore, a statistical formalised 5D (x, y, z, t, M – geographical location, depth, time, magnitude) model of seismic process, which is based only on seismic data, has been developed. The seismic catalogue may be considered as a 5D table. Every event from a general set should belong to at least five subsets with the linear relations between location
Figure 6. Example of application of the preliminary 4D-model for the northern part of Taiwan
Island. (a) Absolute vector of deformation; dark areas correspond to the zones of tension, light areas denote the compression zones. The epicentre of the Chi-Chi earthquake (1999, reverse faulting) is shown by black star. (b) The fronts of dynamic deformation (FDD) which were modelled for September 1999, depth 13.8 km. Big star – the Chi-Chi earthquake hypocentre; small stars denote earthquakes occurred before and after the Chi-Chi earthquake mainshock.
and time of occurrence. For a single statistical subset, the following generalised equation may be written
M = aX + bY + cZ + dT + e, (12)
where (a, b, c, d, e) are the unknown parameters, and (x, y, z, t, M) are the para-meters of seismic event in catalogue. Obviously, it is not possible to describe the whole seismic catalogue by a single equation. The algorithm of determination of the parameters{ak, b,ck, dk, ek, k= 1, n} represents the essence of approximation
The basic Equation (12) contains five unknown parameters{a, b, c, d, e} and it is quite enough to use five seismic events that satisfy the condition
M1= aX1+ bY1+ cZ1+ dT1+ e M2= aX2+ bY2+ cZ2+ dT2+ e M3= aX3+ bY3+ cZ3+ dT3+ e M4= aX4+ bY4+ cZ4+ dT4+ e M5= aX5+ bY5+ cZ5+ dT5+ e . (13)
If the parameters of linear relationship{a, b, c, d} are found, it possible to select all events that satisfy the condition
|M − aX + bY + cZ + dT + e| ≤ , (14)
here is the a priori accuracy of approximation of the catalogue. The iterative algorithm, which belongs to non-linear simplex optimisation, has been constructed for consecutive selection of the particular subset (12) in seismic catalogue. Every seismic event should belong to at least five independent subsets. The basic events are revealed by direct sorting of all possible combinations (five events from the initial catalogue). For the case of a very large number of events in the catalogue, the method of stochastic search for the best solution is used. The solution is found from a subset, which has been selected, in stochastic manner, from the entire large data set. It is also possible to apply a criterion of weighted minimal Euclidean distance between seismic event and 5D plane {x, y, z, t, magnitude}. The weights assigned to particular events are determined on the basis of analysis of magnitude distribu-tion in the studied catalogue. Small weights correspond to small magnitudes and vice versa.
The model is considered statistical, but not in the traditional sense, bearing in mind general statistics of seismic catalogue. All statistics for linear relations between parameters are calculated after special separation of the catalogue into subsets. The number of these subsets is a minimal quantity, which means that it is not possible to separate the initial catalogue into a lower number of subsets with the same accuracy. The subsets are interconnected with each other in such way that every event in the catalogue does belong to five or more subsets. This provides a basis for optimal statistical continuation of catalogue. Every event may be tested whether it does belong to the already constructed subsets with certain statistical relations. The event will be assigned to, at least, five subsets simultaneously, and parameters of the subsets may be slightly re-evaluated. On the other hand, the particular event jointly with other events may produce a new statistical subset. However, in condition of large number of events and statistical relations, the second case is believed to be an extreme rarity. Practically all the potential seismic events should satisfy to the already revealed statistical relations in the initial catalogue. Thus, if an event, which does belong to five or more subsets, has been found in space-time domain and this point does not belong to the initial catalogue, it
means that this new event is characterised by a very high statistical probability of occurrence in correspondent space-time domain.
When using the dynamic model of seismicity (5D-model), the process of predic-tion consists in solupredic-tion of a sequence of local inverse problem based on equapredic-tion (12). Let us assume that we know parameters of five subsets{a, b, c, d, e}, which determine the seismic event to be predicted. Then we can construct the following system of linear equations
M= a1X+ b1Y + c1Z+ d1T + e1 M= a2X+ b2Y + c2Z+ d2T + e2 M= a3X+ b3Y + c3Z+ d3T + e3 M= a4X+ b4Y + c4Z+ d4T + e4 M= a5X+ b5Y + c5Z+ d5T + e5 . (15)
Solution of the system produces parameters{x, y, z, t, M} of the required event. All seismic events in the initial catalogue are described by at least five relationships (12), therefore it is possible to calculate the “theoretical” parameters for every event in the initial catalogue and to estimate the accuracy of approximation. The points of intersection of the linear subsets, which do not belong to the initial catalogue, are considered as the potential seismic events or prognostic earthquakes. The complete set of the prognostic earthquakes is evaluated using the algorithm that generates all combination of linear functions (12).
When the 4D-model for process of dynamic deformation and the 5D-model for seismic process are used jointly for prediction purposes (Figure 7a), it is necessary to determine the peculiar points that coincide each other in the both models (see also Ovcharenko, 1998; Ovcharenko et al., 2000, 2001). First, the 5D-model is applied to find a potential event (x, y, z, t, M). Second, a time selection is performed, and the event (x, y, z, t, M) is tested whether the corresponding point exists in the 4D-model. The testing is performed by analysis of all fronts of dynamic deformation (FDD) and by calculation of the criterion
|Y − αX + βZ + γ T + λ| ≤ , (16)
where (α, β, γ , λ) are the kinematic parameters of the FDDs; is the model accur-acy. The point does belong to the 4D-model if the criterion is fulfilled for at least four FDDs. Thus, if the point does belong to the both models (5D and 4D) it is considered as a prognostic event. The magnitude of the event is evaluated by the 5D-model, by definition.
4.2.3. Application of the Models for Description of Seismic Process
In this section we are going to discuss the question “how well the complex pro-cesses of the Earth’s crust deformation and seismicity may be described by the developed models?” using several typical examples of application of the models.
Figure 7. Modelling of seismic process for the Taiwan region. (a) Scheme of the forecasting
using the 4D- and 5D-models. (b) Comparison of location of epicentres of the real earthquakes (upper picture) and the modelled events (lower picture) for the Taiwan region (1900–1998). (c) Evaluation of accuracy of approximation (magnitude and location – latitude) of the real seismic catalogue (Taiwan region, 1900-1999) by the modelling. (d) Comparison of total number of large observed earthquakes and theoretical events. (e) Comparison of paramet-ers (location, time, and magnitude) of strong earthquakes (stars) occurred in Taiwan region during September–December 1999 (aftershocks of the Chi-Chi earthquake are not shown) and possible seismic events predicted by the 4D/5D-models (crosses).
The 4D-model was constructed using the analytical basis described in Appendix 1 and using various types of observed data. When analysing correspondence between the initial data and the result of modelling, it is possible to compare time-dependent series (GPS, or PSWL data) or distribution of various parameters on the studied area. On the other hand, it is possible to compare location of the empirically de-rived “indicator points” and results of the modelling. In this paper we consider the second case. For the Taiwan region, the constructed 4D-model contains 359 FDDs and 1795 unknown parameters (359 amplitude and 1436 kinematic parameters). The accuracy of evaluation of the kinematic parameters of the model is estimated using the catalogue of geodynamic indicators or “indicator points”. The catalogue is compiled by analysis of GPS data, sea water level observations, generalised seismic catalogue, etc. Table I lists statistical characteristics of approximation of the “indicator” catalogue by the modelling.
From the seismological point of view, the main goal of modelling of the process of dynamic deformation and seismicity is to compile a complete (M > 2–3) prognostic or theoretical seismic catalogue for a sufficiently long time period. The catalogue may be further used for purposes of seismic hazard zonation. The schemes of seismic source zones are not strongly influenced by statistical errors of prognostic time of occurrence, location and magnitude of the events. Obviously, the models intended to be used for forecasting of future events should adequately describe the entire set of occurred events.
Figure 7. Continued.
Table I. Statistical characteristics of the difference between
parameters of observed and modelled “indicator points”, Taiwan region
Parameter Standard deviation Maximum absolute (32 319 points) errors (97%)∗
Longitude, deg 0.0025 0.22
Latitude, deg 0.001 0.2
Depth, km 2.35 20
Time, year 0.63 15
∗Difference between parameters of observed and modelled events for the subset containing 97% of all points (31 350 points).
Figure 7. Continued.
The generalised catalogue for the region contains more than 29000 events (M > 3). The regional seismic catalogue collected by Central Weather Bureau contains information only for limited area (119–123E, 21–26N). Therefore we supplemen-ted the data by NEIS catalogue (http://neic.usgs.gov/neis/bulletin/bulletin.html). The depth for many shallow earthquakes in the NEIS data is equal to 33 km that can be seen in the Figure 7b. At the same time, we should note that the general-ised catalogue is incomplete for interval 1900–1972 when there was no telemetric seismographic network in the region. We assume that the notion does not concern intermediate and large earthquakes (M > 5).
The results of the modelling (Figures 7b, c, d; Table II) of the observed cata-logue (1900–1999) by the developed models (Sokolov et al., 2003b) show that the
Table II. Statistical characteristics of the difference between parameters
of observed and modelled (5D-model) earthquakes, Taiwan region Parameter Standard deviation Maximum absolute
(29 986 events ML= 3–8) errors (97%)∗ Longitude, deg 0.0042 −0.4; +0.42 Latitude, deg 0.0067 −0.50; +0.55 Depth, km 2.63 −10, +30 Time, year 0.51 −3.8; +4.5 Magnitude, unit 0.03 −0.3; +0.3
Difference between parameters of observed and modelled events for the subset containing 97% of all events (29 100 events).
overwhelming majority (94–95%) of the events are characterised by the following errors of approximation: location (longitude and latitude) – less than 1–2 km; depth – 2–3 km; date of occurrence – 0.5 years; magnitude – 0.03 units. Maximum ab-solute errors do not exceed 0.5–0.4 degrees for location; 11 km for depth; 4 years for date of occurrence; 0.3 unit for magnitude. However, there are large errors for small number (5–6%) of theoretical events, which are caused by peculiarities of the applied algorithm. Moreover, the observed catalogue contains foreshocks and aftershocks, which could not be entirely considered in the modelling.
Figure 7c,d shows examples of comparison between the observed and modelled seismic catalogues. It is seen (Figure 7c) that the modelled catalogue provides a good agreement with general characteristics (total number of earthquakes of various magnitudes) of the observed catalogue. Let us also compare parameters (location, time of occurrence, and magnitude) of the largest earthquakes occurred in the Taiwan region during September–December 1999 and the theoretical seismic events predicted by the models (Figure 7d). Location of extended source of the Chi-Chi earthquake with non-uniform distribution of dislocation coincides with two large theoretical consecutive events. Magnitude, location and time of the other strong earthquake, which occurred near the SE coast, almost completely fit the parameters of the predicted event. However, the earthquake of M = 6.2 that oc-curred a month later to the South from the Chi-Chi earthquake area does not exist in the theoretical catalogue. It is possible to assume that large earthquake could produce a local fast-propagating spherical front of deformation (post-earthquake front). The front, in turn, causes new events due to interaction with the exist-ing FDDs. Therefore, the aftershocks may be divided, at least, into two groups, namely: “triggered” events that are caused by the additional post-earthquake front and “independent” earthquakes resulted from interaction of the existing FDDs.
4.3. THEORETICAL CATALOGUE AND PARAMETERS OF FUTURE LARGE EARTHQUAKES
The theoretical catalogue for future 40–50 years for Taiwan region (119–123 E, 21–26 N) was calculated using the described models. When developing the mod-els, we used GPS data (global and local networks), sea water level observation, and generalised seismic catalogue provided by National Weather Bureau, which contains more than 29000 events (M > 3) and covers time period from 1900 to 1999.
The theoretical catalogue has been compiled on the basis of empirical data obtained before June 2001 (2001.5). When the paper was prepared, it was pos-sible, for testing purposes, to compare parameters of observed earthquakes and theoretical events for time period between January 2001 and June 2003 (the cata-logue of observed events for January–December 2002 contains earthquakes with magnitudes more 4.0). It has been noted above that real seismic event may occur before or after the corresponding theoretical event. Thus, when comparing real and theoretical catalogues, it is necessary to use the theoretical catalogue calcu-lated for a period that covered a sufficiently long interval before and after the real observation. The condition is not realised at present: we have a theoretical catalogue that covers only the interval after the observation. Nevertheless we made direct comparison between observed and theoretical catalogues (Figure 8) for the following time periods: observed earthquakes – July 2000–June 2003 (more than 480 events with M > 4.5); theoretical events – January 2001–June 2003 (more than 550 events with M > 4.0). There is, in general, a good agreement between groups (clusters) of observed earthquakes and theoretical events in the central part of Taiwan and eastern coastal areas. The discrepancies are observed for the remote boundary zones where the necessary observational (input) data are insufficient and even absent.
Actually, the comparison of large observed and theoretical seismic catalogues is not a simple procedure. It necessary to resolve the problem: how it is possible to identify two correspondent events in 5D-space (geographic location, depth, time and magnitude)? The following approach was used. Let us denote the parameters of event as follows: longitude as X1, latitude as X2, depth as X3, time as X4, and magnitude as X5. During the first step, all parameters of the catalogues are normalised and centralised as
x = (x − x0)/σx, σx= (x− x0)2/n, x 0= n i=1 xi /n. (17)
The “distance”, or a quantity related to correspondence of the events, between two events is evaluated as follows
Figure 8. Comparison of location of observed earthquakes (epicentres, ML > 4.5, July 2000–June 2003) and events in the theoretical catalogue (M > 4.5, January 2001–June 2003. Upper picture – location of epicentres; lower picture – generalised vertical projection (South-North).
δE(xF, xP ) = 5 i=1 (xiF − xiP)2/5 δM(xF, xP ) = 5 i=1 |(xiF − xiP)| /5. (18)
What events (xiF, i = 1, 5), (xiP, i = 1, 5) should be used in these equations?
Obviously, the first one is an arbitrary earthquake in the observed catalogue. The second one – is the nearest, in terms of (17), event from theoretical catalogue. Thus, when evaluating the error of prognosis, we should determine “the best prognostic” events as
Mink=1,mδE(xF, xP ), Mink−1,mδM(xF, xP ). (19)
Here m is the number of events in theoretical catalogue. Note that all parameters of earthquakes (location, depth, time, and magnitude) in Equation (18) are char-acterised by the equal weights (influence on the choice). Obviously, the different weights may be chosen to increase or decrease influence of a specific parameter.
Let us increase the influence of the time parameter on the choice of “best prognostic” event by increasing the weight (WT = 10). Table III summarises
statistical characteristics of difference between parameters of observed earthquakes (66 events of magnitudes > 4.99) and corresponding (“best prognostic”) events in theoretical catalogue. In this case, some Real-Theoretical (R-T) pairs are charac-terised by relatively large difference (more than 40–50 km) between location of earthquakes. The averaged relative differences between location and magnitude of the real and theoretical earthquakes are close to zero. The averaged absolute differences do not exceed 10–20 km for location, 14–16 km for depth, 0.1 unit for magnitude and 1.0 years for date of occurrence.
Figure 9 compares location of epicentres of large (M ≥ 5.9) earthquakes oc-curred in the region during the considered period and the correspondent (WT=10) events in theoretical catalogue. The difference between location almost for all R-T pairs does not exceed 0.3-0.4 degree (30-40 km). Three R-T pairs are characterised by large difference between location (more than 0.5 degree) and the real earth-quakes are shown as open stars. Two real earthearth-quakes (M 6.7 and 5.9) of these pairs occurred in the boundary zones of the studied area. The third real earthquake (M = 6.1) occurred in the central part of the island in July 2000 and, most prob-ably, it should be also considered as aftershock of the Chi-Chi (September 1999) earthquake.
It is necessary to note that the hypocentre of a real earthquake is determined as a point of rupture initiation. Location of the point on the whole ruptured area (earthquake source) is unpredictable. The models determine location and depth of
Table III. Statistical characteristics of the difference between parameters of
observed and modelled earthquakes, M > 5.0, Taiwan, July 2000 – June 2003,
WT = 10 (see text)
Parameter Relative Relative Absolute Absolute difference, difference, difference, difference, earthquakes earthquakes earthquakes earthquakes with M > 5.0 with M > 5.5 with M > 5.0 with M > 5.5
Longitude, deg 0.00 0.03 0.11 0.19
Latitude, deg 0.01 0.03 0.12 0.15
Depth, km −2.75 −6.5 13.0 14.0
Time, year −0.4 −0.65 0.84 1.0
Magnitude, unit 0.05 0.08 0.11 0.14
a possible seismic event as idealised point that corresponds to interaction of four or more fronts of dynamic deformation (FDD). Thus, the prediction should be con-sidered as successful, if the difference between hypocentre of observed earthquake and location of theoretical seismic event does not exceed largest dimension of the earthquake source.
Let us consider the following example. One of the largest earthquake during the considered period (ML= 6.8, 31 March 2002) occurred approximately 100 km to
south-east from the city of Taipei (see Figure 9). The earthquake caused collapsing of residential building in the city and falling of construction facilities from the top (56 stories) of Taipei Financial Center that is currently under construction. Five construction workers were killed. The corresponding theoretical event (ML= 6.5,
August-September 2002) is located at 30–40 km to the south-eastern direction, and the source dimensions of earthquake of such magnitude (ML=6.5-6.8) are larger than the difference of location. It is necessary also to note that the earthquake occurred in the boundary zone of the studied area where the differences between parameters of observed and theoretical events are larger than that for the central and eastern part of the island.
Strictly speaking, instead of idealised point, it is necessary to consider a certain volume of the crust that is affected by the fronts of dynamic deformation (FDDs). Depending on the parameters of the FDDs, the volume is characterised by complex distribution of deformation and stress. Bearing in mind the unpredictable charac-ter of rupture initiation and propagation that depend on mechanical and chemical processes in fractured, porous, fluid-saturated, chemically evolving, thermoplastic crust, it is possible to consider the following phenomenon: a single real earthquake may correspond to several theoretical events. By the other words, the rupturing zone of the real earthquake may extend through several points that are considered as theoretical events and that may occur within a certain volume during a certain
Figure 9. Comparison of location of large (M > 5.9) observed earthquakes (stars, July
2000–June 2003) and events in the theoretical catalogue (triangles).
time period (see Figure 7d). In this case, when comparing the real and theoretical catalogues, it is necessary to check whether it is possible to find the correspondence between a single real earthquake and a group of smaller and nearby theoretical events that should occur within a certain time period around the real earthquake.
Let us consider distribution of summarised seismic moments (SSM) of the-oretical events. Average error of evaluation of the time parameter for large and intermediate earthquake was estimated as 1 year and average error location is about 20 km (Table III). Therefore, the values of seismic moments were calculated, within volume of crust with dimension 0.4 degrees × 0.4 degrees in plane and 40 km in depth, for all theoretical events that may occur during time period of 2 years. The well-known relation between moment magnitude MW and seismic
moment M0 was used, namely: MW = 2/3 log10M0 − 10.7. Figure 10a shows distribution of SSM values (base 10 logarithm) along the territory for shallow events (depth 0–40 km) and location of large (M > 5.9) real earthquakes. All real large earthquakes (up to June 2003) did occur within the areas with increased SSM values: log10SSM > 25 corresponds to earthquakes with MW > 6.0. The
distribution of SSM values for the nearest future is shown in Figure 10b. Using the information, one can make suggestion about development of seismic process in the region. On the other hand, the schemes calculated for various consecutive periods of time may be used as a basis for time-dependent hazard assessment. In this case, the SSM value may be considered as seismic moment of “maximum magnitude” or “characteristic” earthquake, which can occur within the given elementary volume.
5. Time-Dependent Hazard for the Taiwan Region
The region and site and time-dependent seismic analysis is based on schemes of probable earthquake zones, regional source scaling and attenuation models, and information on local site response on seismic motion. The following scheme was applied for the analysis (Figure 11).
The characteristic earthquakes were used for calculation of ground motion parameters along the territory of Taiwan Island. We should take into account un-certainty in location, geometry and dimensions of future earthquake sources. The uncertainty was considered as follows: five points inside the elementary volume were accepted, with equal probability, as location of the characteristic events. For shallow (depth less 40 km) earthquakes, two values of characteristic depth (5 km and 15 km) were used. For deep earthquakes we used a single value of character-istic depth – 40 km. The source-to-site distances, ten (shallow events) or five (deep events) values for every site and elementary volume, were evaluated using these parameters of the source location.
Ground motion parameters (Peak Ground Acceleration and Response Spectra at selected frequencies), which may be expected during future earthquakes, were evaluated for grid points (10 km× 10 km) covered the territory of Taiwan Island. For every grid point the probability that ground motion parameter X will not exceed a given value x may be estimated as follows:
P(M=m;R=r)[X ≤ x] = 1 σx √ 2π x xmin exp{(x − a)2/2σx2}dx, (20) where a is the mean value of log10X for an earthquake of given M and R; σx
is standard deviation describing the scatter of ground motion parameter for the earthquake, and xmin is of sufficiently small value (xmin ≈ a − 5σx). As far as
we have ten or five values of characteristic distance, the resulting function P was calculated as an average of the particular functions Pi.
Peak Ground Acceleration (PGA) values and Response Spectra amplitudes (RSA) at selected frequencies (parameter a in Equation (20) were calculated on
Figure 10. Distribution of summarised seismic moment values (base 10 logarithm) for
shallow theoretical events (depth 0–40 km). Considered time period includes observations. Epicentres of large observed earthquakes, magnitude and date of occurrence are shown for comparison. The dates of earthquakes are shown in decimal format: for example, June 2002 corresponds to 2002.5.
Figure 11. Scheme of time-dependent seismic hazard assessment.
the basis of regional spectral models (FAS) and generalised site response functions (mean amplitude values, soil classes B, C and D) using stochastic approach. A set of 40 synthetic acceleration time functions was generated for every (M, R) pair using effective duration (regional estimations). The resulting parameters (PGA and RSA) were estimated as the average values calculated from the set.
The cumulative probability functions P for the grid point (observation site) were calculated from all (M, R) pairs (elementary volumes of crust). The values of ground motion parameters were determined for three values of probability, namely: 0.5 (mean value), or 50% of being exceeded during given (M, R) event; 0.84 (mean + 1 standard deviation), or 16% of being exceeded during given (M, R) event; 0.97
(mean + 2 standard deviation), or 3% of being exceeded during given (M, R) event. The maximum value of ground motion parameter resulting from all considered (M, R)pairs was finally assigned to the site.
It can be seen that the approach may be considered as a modification of so-called “scenario earthquake” analysis, or as a variant of “logic tree” approach. In this case we did not consider the distribution of earthquakes in time in a probabilistic manner. However, it is possible to introduce “the time dependency” by calculation of the ground motion parameters distribution for the certain periods of time. In this study we determined maximum value of ground motion parameter for the points of the grid using all considered periods (2003–2025). It is possible to apply various weights for various time periods that should reflect the increasing uncertainty of earthquake prediction in the course of time. However, in this study we apply similar weights for all considered time periods.
Figures 12 and 13 present, as an example, the seismic hazard zonation maps in terms of Peak Ground Acceleration and Response Spectra (the contours are given in cm/s2). The schemes show distribution of amplitudes that will not be exceeded with 84% probability in condition of generalised soil classes during all possible earthquakes that may occur in the region during time period of 2003– 2025. We should note that all earthquakes in this study were considered as point sources, even with five possible locations, therefore the contours are characterised by a concentricity. Of course, the 3D of representation the source (length, width, strike and deep characteristics) is more reliable. However, for this case a study of predominant source parameters in the Taiwan region should precede the ground motion calculations. This is one topic for future research.
The approach used in the study allows us to introduce a new parameter that describes dependency of seismic hazard on time, so-called “period of maximum hazard” or PMH (Figure 14), which is of particular interest of seismic risk man-agement and insurance business. As far as we have several schemes of seismic zonation (source zones) for various time periods, it is possible to evaluate the period, during which every considered site will be subjected by the maximum value of ground motion parameter. When using jointly, the three types of zona-tion maps (future seismic source zones (Figure 10), distribuzona-tion of ground mozona-tion parameters (Figures 12, 13), and period of maximum hazard (Figure 14)) allow optimisation of engineering decisions, and may be considered as a basis of seismic code development.
6. Conclusion
The research described in this paper may be considered as the parts of uni-fied approach to region and site and time-dependent seismic hazard assessments, namely: from evaluation of characteristics of future earthquakes to ground motion parameters, which are used in seismic design and risk management.
Figure 12. Distribution of Peak Ground Acceleration values (cm/s2), which will not be ex-ceeded with probability P (84%) in conditions of generalised soil classes (B and D) during time period of 2003–2025.
The theoretical seismic catalogue is calculated using the dynamic modelling of the Earth’ crust deformation (4D-model) and seismicity (5D-model). Various available geophysical and geodetic data (GPS, sea water level observation, seis-mic catalogues, etc.) are used for creation of the models. The developed models adequately describe the observed seismicity in the studied region. The theoretical catalogue is used as a basis for future earthquake zonation and determination of maximum-magnitude earthquakes for various time periods and various locations.
When evaluating ground shaking parameters for purposes of seismic hazard and seismic risk assessment, we suggest using Fourier amplitude spectrum of ground acceleration (FAS) as a universal input parameter in seismic hazard analysis, both deterministic and probabilistic ones. For the case of Taiwan region, we developed empirical models of ground motion spectra on the basis of several thousands re-cords from recent (1993–1999) earthquakes. The database includes strong-motion data collected during the Chi-Chi earthquake (M = 7.6, 21 September 1999) and large (M= 6.8) aftershocks.
Empirical amplification functions for site classes B, C and D in Taiwan were evaluated as the ratio between Fourier amplitude spectra of recorded accelerograms and spectra modelled for hypothetical very hard rock (VHR). The amplification functions demonstrated their reliability when comparing with independent data. It
Figure 13. Distribution of Response Spectra amplitudes (cm/s2) for various frequencies (peri-ods), which will not be exceeded with probability P (50%) in conditions of generalised soil classes (B and D) during time period of 2003–2025.
Figure 14. Distribution of “period of maximum hazard” (PMH) – the parameter that describes
dependency of seismic hazard on time. PMH denote the period, during which every considered site will be subjected by the maximum value of considered ground motion parameter (PGA in this case).
is also necessary to note the prominent influence of geologic and geomorphologic factors on site amplification function for the considered cases of rock (class B), soft rock or very dense soil (class C), and stiff soil (class D) sites. The influence is reflected by large variations of amplitudes and predominant frequencies between particular stations with the same site class that was assigned on the basis of the rock age and geological classification. Therefore, as it has been noted by Lee et al. (2001), further studies on site classification should be carried out using more spe-cific subsurface geotechnical data including thickness of soil, shear wave velocity, and density.
The schemes of probable (future) earthquake source zones, regional source scaling and attenuation models, and information of local site response are used for region and site and time-dependent seismic hazard analysis. The approach al-lows us to introduce a new parameter that describes dependency of seismic hazard on time, so-called “period of maximum hazard” or PMH, which is of particular interest of seismic risk management and insurance business.
We should note that the results obtained in this study (Seismic Hazard Zonation) should be considered as preliminary variant. Bearing in mind several shortcomings and unresolved problems, the described approach should be developed further, both in general and particular aspects. The future tasks may be outlined as follows.
1. Evaluation of possible sources of the wave-like objects (FDDs), which may consist in the following. Nonlinear differential equation of motion is analysed for the case of heterogeneities of density that are located inside the crust and that are influenced by varying gravity field. The numerical solution of the equa-tion allows checking possibility of appearance of slow undulate deformaequa-tional processes inside inhomogeneous crust.
2. Bearing in mind dynamic nature of the process of the Earth crust deform-ation, the 4D- and 5D-models should be continually supplemented by new data. The prognosis should be revised after every large earthquake and when accumulating new experimental data (2–3 years).
3. Study of typical parameters of earthquake sources (dimensions, orientation, etc.) for different areas in the Taiwan region to be used in strong ground motion modelling.
4. Analysis of peculiarities of seismic waves attenuation and propagation from shallow large earthquakes located in different areas, for example, near the north-eastern coast.
5. Studies on site classification using more specific subsurface geotechnical data including thickness of soil, shear wave velocity, and density. Evaluation of generalised site response characteristics for classes E and area-dependent characteristics.
6. Detailed analysis of response of alluvium filled valleys (Taipei basin, Ilan area, etc.) on earthquakes occurred in different location and direction from the basin. 7. Analysis of peculiarities of soil response during strong seismic excitation and quantitative description of possible non-linear soil response in terms of site amplification function, consideration of non-linear effects during large earthquakes.
Acknowledgements
The constructive comments and suggestions from anonymous reviewers, which helped improve the article, are gratefully acknowledged. The authors would like to express their gratitude to National Science Council of the Republic of China for opportunity to perform the study, as well as to scientific workers and employees of National Center for Research on Earthquake Engineering and National Taiwan University. We are very grateful to Mr. Hsu Chien-Chin, Dr. Jean Wen-Yu and Mr. Huang Shin-Kai for providing valuable and friendly assistance. Ground-motion records of the Chi-Chi earthquake sequence were provided by National Weather Bureau of the Republic of China, and local GPS data – by Prof. Yu Shui-Beih, Institute of Earth Science, Academia Sinica.
