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 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.
Appendix 1: 4D-Model
In the frame of the 4D-model, we assume that the FDDs propagate inside the crust with no dissipation and parameters of a particular FDD (amplitude, velocity, width, displacement distribution, etc.) are the same along the whole extent of the FDD and depth. Bearing in mind these restrictions, the formulation of the model for the process of the Earth’s crust deformation may be outlined as follows.
When crystalline structure of a solid body contains a defect, there is a corpus-cular dislocation around the defect – so-called “dislocation Frenkel-Kantorova”. It has been shown (Frenkel, 1959) that the dislocation, dimension of which greatly exceeds the distance between atoms, could easily move inside the crystal. The function describing the propagating dislocation (1-D case) may be written as
ϕ(t, x)= π − 4arctg[e−(t−x/v)/lv] (A1)
ϕ(t, x)= −(t − x/v)/lv
e−(t−x/v)/lv+ e(t−x/v)/lv, (A2)
where v is the velocity of propagation, lv is half-width of the dislocation that depends on the velocity v, ϕ(t, x) is the angular measure of the dislocation dis-placement with respect to equilibrium lv= l0
1− v2/v20, v0is the sound velocity in the crystal.
Various wave equations including so-called KdV equation (Korteweg and de Vries, 1895) were studied by many authors, starting from Russell (1845), for a liquid medium. The solution of the KdV equation allows to find so-called Russell or KdV soliton.
The Russell’s soliton, or solitary wave, is described by exact solution of the KdV equation (Kruskal et al., 1970)
y(t, x)= y0/ch2
x− vt l
. (A4)
Let us transform the equation to the form that is usual for application in geophysics When comparing the equations for the Russell and Frenkel–Kantorova solitons, it could be seen that they are essentially the same. Thus, the idealisation of the
Earth’ crust as a liquid or elastic medium will not result in significant errors when describing the processes of deformation by a propagating dislocation.
When the deformation of the Earth’ crust is considered as a dislocation inside a crystal, it is assumed that the parts of the crust could move relatively each other due to tectonic fracturing. The faults, fractures and cracks, therefore, could be considered as the “atoms” of dislocation. Relative displacements up to 10–100 meters along the distance of 10–20 km, and more than 100 meters for the larger distances, may appear through the systems of faults and cracks. The values of rel-ative deformation are about 10−5–10−6that are comparable with the data observed in reality. Let us consider a simplified form of equations (A1)–(A2)
ϕ(t, x) = −(t − x/v)/lv
This asymptotic equation for a single object coincides with the equations for dynamic deformational displacements, which are applied in this work. The equa-tions were obtained by Ovcharenko (1997, 1998) on the basis of heuristic reasons with no consideration of corresponding differential equation of motion for a dis-location. The important difference consists in the use of multidimensional space (x) = (x, y, z). The typical equation for Russell–Frenkel–Kantorova’s soliton has been simplified by decomposition of the exponential member in denominator into the power series, and by excluding the series members of more than fourth exponent. The simplification however leads to a slower attenuation of the deform-ational displacements with distance that is not significant in our case. As it has been mentioned above, we do not consider all properties of the classic soliton and we use only asymptotic approximation. Therefore, we apply the term “front of deformation” (FDD).
Thus, considering the aforementioned statements, the potential function of elastic strain displacement for dynamic model of a set of plane FDDs, propagating within the infinite media, could be described as a simple summation
E(x, y, z, t)=
Description of the parameters is given in Table A1.
In many cases the observation networks are characterised by sparse location of stations and duration of observation does not exceed 10 years. Thus we have to utilise various information including numerous “indirect” geophysical and geodetic data, which are characterised by longest period of observation, in some cases more than 100 years. Due to technical reasons, different quantities (vertical and hori-zontal displacement, velocity of displacement, deformation, variations of distance,
Table AI. Basic formulae for 4D-model of dynamic deformation Data Element Time T0 Equation
(yes= T0,
inclination) had being measured in different periods of time. The basic mathemat-ical expressions should cover all kinds of the data we use in modelling. Table AI summarised the interconnected formulae, which are based on description of the process of deformation by a system of plane fronts of dynamic deformation (FDD).
Besides these quantities (Table AI), which can be directly observed or calcu-lated from the observations, during the modelling and analysis there is a necessity to calculate various characteristics of dynamic tensor of deformation Tdef.
Tdef=
It is supposed that the all variables are functions of location and time. The in-variants of the tensor of deformations (E1, E2, E3)which also depend on location and time, are important characteristics of the process
E1 = x+ y+ z
The invariants may be calculated by numerical differentiation of displacement functions. For example, for the first invariant we have the following
E1(x, y, z, t) ∼=
As far as the displacement functions are known in analytical mode after determ-ination of kinematic parameters (number of the waves, velocities and azimuthal parameters of propagation) of the model and after evaluation of amplitudes of the FDDs, there are no problems with numerical differentiation. Finally, for dilatation or first invariant of the tensor of deformation we have the following
E1(x, y, z, t) ∼=
where coef is the coefficient for transition between distance (e.g., kilometres) and displacement (e.g., millimetres).
For calculation of variation of distance between two stations (Table AI, VLBI) it is possible to use the following approximate equation
δL = L(x, y, z, t) − L0=
As it follows from the formal mathematical description, it is necessary to eval-uate a set of parameters{a4i, b4i, c4i, d4i, m4i, i = 1, N}. At the same time we do not know the number of parameters N in advance. For correct evaluation of the parameters, which describe the model, the procedure should satisfy to several con-ditions. The general condition is the following: the number of unknown parameters should be less (or equal) than the number of observations.
In the frame of considered model and the used data, the basic optimising functional may be written in the following way
M[z,α,βi] = β1· SGPSlon+ β2· SGPSlat+ β3· SGPSVE + β4· SGPSVN
β5· SLevelling+ β6· SVLBI+ β7· Sσ x/σ y+ α. (A13) where M is the Tikhonov’ regularising function (Tikhonov and Arsenin, 1979); z is the set of unknown parameters or the solution to be found; α is the parameter of Tikhonov’ regularisation; βi is the weight that considers influence of various types of observed data on the solution; Siis the mean square residual function for various types of observed data Si = fi − Fi(z)2; fi is the given type of the used data (modelled); Fi is the type observed data; is the Tikhonov’ stabilising function, which is accepted in the simplest form = z 2.
Optimal consideration of influence of various types of observed data is caused by
βi ≈ 1/SiNI, β1/S1≈ β2/S2≈ · · · ≈ βi/Si ≈ . (A14) The iteration procedure, which is used in the modelling, could provide the condition easily.
Let us consider the structure of non-linear inverse problem of optimisation. If we know the set of parameters{a4i, b4i, c4i, d4i, i= 1, N}, then the problem (A13) is transformed to a linear inverse problem relatively the group {m4i, i = 1, N}.
The transformation provides a basis for decomposition of the general task into two stages, namely: evaluation of kinematic parameters (characteristics of propagation and orientation) and evaluation of amplitude characteristic (elastic displacement) of the model. The inverse kinematic problem, as a completely non-linear one, requires a special procedure. When the kinematic parameters are found, the amplitude char-acteristics may be evaluated using standard procedure for linear multi-parametric regularisation (Tikhonov and Arsenin, 1979).
For solution of the first stage of general problem we use specific transformation of dynamic objects that is described by equations in Table AI. The formulae take specific values (zero or local extreme) when a single front of dynamic deformation (FDD) passes through the given point of observation. If the peculiar (or indicator) local points are found during preliminary processing of the observed data, the evaluation of kinematic parameters may be sufficiently simplified. Description of technique for detection of the indicator points for various kinds of observed data, as well as the iterative procedure for evaluation of kinematic parameters, may be found in (Ovcharenko 1995, 1997a, b, 1998, 1999; Nusipov and Ovcharenko 1997;
Ovcharenko et al. 2000). Here, when speaking about the indicator points, we mean a set of points in space and time domains that hold fixed parameters (location and time moment) of certain local peculiarity of the geodynamic phenomenon (see description of the term in Section 4.1 “Basic principles”).
The procedure of evaluation of kinematic parameters using the set of indicator points (IP) may be shortly described as follows. The translational movement of solitary plane wave of deformation may be described by parametric equation for the axis plane of maximum tension or compression.
y = ax + bz + ct + d, (A15)
where (a, b, c, d) are the arbitrary constant values; (x, y, z) are the Cartesian co-ordinates; t is the time parameter. Parameter c is the velocity Vyof the axis plane.
where (a, b, c, d) are the arbitrary constant values; (x, y, z) are the Cartesian co-ordinates; t is the time parameter. Parameter c is the velocity Vyof the axis plane.