Simulation of Near-Field Strong Ground Motions Using Hybrid Method

1. Introduction

Earthquake disaster investigations have shown that numerous strong earthquakes are caused by remobilization of active faults. Many casualties and severe damages to structures as well as huge economic losses have resulted from ground motions of strong earthquakes caused by active faults buried under urban areas. Recently, both potential hazard and defenses of active faults concealed under urban area has become a grand research subject paid highly attention to by the seismologists. Near-field strong ground motions, especially their high frequency content, are intensively affected by both slip heterogeneity on fault plane and rupture process of an earthquake fault. In the simulations of near-field strong ground motions, modeling effective finite fault source is very important. The gradually increasing number of recorded near source time histories has recently enabled strong motion seismologists to analyze more precisely the character of the near-fault ground motions and therefore contribute to the physical understanding of those features that control them (Malagnini et al. 2002, 2011; Akinci et al. 2010; D’Amico et al. 2010). Mavroeidis and Papageorgiou (2002) presented a comprehensive review and study of the factors that influence the near-source ground motions.

The stochastic method of synthesizing ground motion based on seismology interests engineers specifically in simulating higher-frequency ground motions (Akinci et al. 2001). The method is widely used to predict ground motions for regions, in which ground motion recordings from past earthquake are not available (Boore 2003). For far field, the point source model of Boore (1983) is very effective; however, for near-field, the method can not incorporate the factors which have significant effect on the near-field strong ground motions, and yields an overestimation of such ground motions. The stochastic ground motion modelling technique, also known as the band limited white-noise method, has been first described by Boore (1983). Ever since, many researchers have applied the method to simulate ground motions from point sources (e.g., Boore and Atkinson 1987; Atkinson and Boore 1995; Zafarani et al. 2005; D’Amico et al. 2012).

On the other hand, realistic acceleration time-histories should be employed in structural analysis to reduce the uncertainties in estimating the standard engineering parameters (Hutchings 1994), particularly for non-linear seismic behavior of structures. Thus, designers need to know the dynamic characteristics of predicted ground motion consistent with source rupture for a particular site to be able to adequately design an earthquake-resistant structure. Hall et al. (1995), Makris (1997), Chopra and Chintanapakdee (2001), Zhang and Iwan (2002) have experimentally as well as analytically studied the elastic and inelastic response of engineering structures subjected to actual near-fault records or simplified waveforms intending to represent the typical ground motion pulses observed in near-field regions.

During the past decades, much effort has been given in reliable simulation of strong ground motion from finite faults through methodologies that include theoretical or semi-empirical modeling of the parameters affecting shape, duration and frequency content of the strong motion records. Due to unavailability of strong recorded ground motion, simulation of ground motion has been carried out using the stochastic method proposed by Boore (2003). The ground motion spectrum has been generated by Atkinson and Boore model (1995). Even though the success of the point-source model has been pointed out repeatedly, it is also well known that it often breaks down, especially near the sources of large earthquakes. Recently, Beresnev and Atkinson (1997) have proposed a technique that overcomes the limitation posed by the hypothesis of a point source. Their technique is based on the original idea of Hartzell (1979) to model large events by the summation of smaller ones. In Beresnev and Atkinson (1997), the high-frequency seismic field near the epicentre of a large earthquake is modeled by subdividing the fault plane into a certain number of sub-elements and summing their contributions, with appropriate time delays, at the observation point. Each element is treated as a point source. A stochastic model is used to calculate the ground motion contribution from each sub-element, while the propagation effects are empirically modeled. Combining the stochastic method with the finite fault source model, Silva (1997), Beresnev and Atkinson (1998), Motazedian and Atkinson (2005) have proposed different methods, which could be effective for simulating or predicting near-field ground motions.

Two Californian earthquake events may be characterized as historical milestones related to near-source ground motions: the 1966 Parkfield and the 1971 San Fernando earthquakes. The 1966 Parkfield, California, event provided the now famous Station 2 (C02) record at a distance of only 80 m from the fault break (Housner and Trifunac 1967). Modern quantitative analysis of strong ground motion observations was started with this record. Aki (1968) and Haskell (1969) demonstrated that the observed transverse (i.e., fault-normal) displacement component of this ground motion record, which exhibited a simple impulsive form, was precisely what is expected for a right-lateral strike-slip rupture propagating from northwest to southeast. The 1971 San Fernando, California, earthquake provided the equally well-known Pacoima Dam (PCD) record. The strike-normal velocity component of this record also exhibited an impulsive character that several investigators attempted to model (e.g., Boore and Zoback 1974; Niazy 1975; Bouchon 1978). In addition, this record was the one that made earthquake engineers recognize the severe implications of the impulsive characteristics of near-source ground motions on flexible structures.

At high frequencies (f > 1 Hz), ground motions become increasingly stochastic in nature. The stochastic methods are generally capable of matching the spectral amplitudes of high frequency ground motions, but are generally not capable of matching the recorded waveforms (Somerville 1998). Firstly, ground motions are estimated by identifying the major regional faults and propagating seismic waves generated at these potential sources to the site of interest. The two commonly used techniques, finite-fault and point source methods of Boore and Atkinson (1987) and Beresnev and Atkinson (1997, 1998) are used for simulation of earthquakes. Both techniques have an omega-squire spectrum.

The main objective of the chapter is to simulate the near-fault strong motion records. Simulation of ground motions is carried out using the hybrid method proposed by Mavroeidis and Papageorgiou (2003) and the stochastic model of Boore (2003). Due to unavailability of strong recorded ground motion, the stochastic method proposed by Boore (2003) is applied to simulate the acceleration time histories. The ground motion spectrum is generated by Atkinson and Boore model (1995). Firstly, macro-source parameters characterizing the whole source area, i.e. global source parameters such as fault length, fault width, rupture area, and average slip on the fault plane are estimated; secondly, slip distributions characterizing heterogeneity or roughness on the fault plane, i.e. local source parameters are reproduced by the hybrid slip model; finally, the finite fault source model, which is developed based on the global and local source parameters is combined with the stochastic method. A simple, yet effective, analytical model proposed by Mavroeidis and Papageorgiou (2003) is also used to adequately describe the impulsive character of near-fault ground motions both qualitatively and quantitatively. The calculated response spectra are compared with those, mentioned in International Building Code (IBC 2000) and Iranian Code of Practice for Seismic Resistance Design Building (Standard No. 2800) to validate the availability and practicability of the proposed method for near-field Tombak site at south-eastern part of Iran. This site includes massive LNG storage plants near to fault. Then, the response of mentioned site under simulated ground motion has been studied by conducting one dimensional ground response analysis. According to the above study, the bed rock and ground surface accelerations of the site are provided. The results can be used in hazard analysis of specific sites in the considered region, particularly for the performance analysis of structures.

2. Simulation method

A simple and powerful method for simulating ground motions is to combine parametric or functional descriptions of the ground motion’s amplitude spectrum with a random phase spectrum modified such that the motion is distributed over a duration related to the earthquake magnitude and to the distance from the source. This method of simulating ground motions often goes by the name ‘‘Stochastic method’’. It is particularly useful for simulating the higher-frequency ground motions of most interest to engineers (generally, f > 1 Hz), and it is widely used to predict ground motions for regions of the world in which recordings of motion from potentially damaging earthquakes are not available. One of the essential characteristics of the method is that it distills what is known about the various factors affecting ground motions (source, path, and site) into simple functional forms.

2.1. Stochastic finite-fault simulation method

In this study, the Stochastic Method is used for simulating the strong ground motion. The method assumes that the far-field accelerations on an elastic half space are band-limited, finite-duration, white Gaussian noise, and that the source spectra are described by a single corner-frequency model whose corner frequency depends on earthquake size (Mayeda and Malagnini 2009). The ground spectrum Y (M₀, R, f) is conveniently broken into several simple functions – the Earthquake Source (E); the Path (P); the Site (G) and the instrument or type of motion (I):

$$Y\left(M_0,R,f\right)=E\left(M_0,f\right)P\left(R,f\right)G\left(f\right)I\left(f\right)$$E1

where, M₀ is the seismic moment, R is the shortest distance from the fault to the site and f is the frequency. Atkinson and Boore model (1995) is used to obtain the ground motion spectrum. The process of strong ground motion simulation is depicted as a flowchart in Figure 1.

Figure 1.

Flow chart showing the structure of the FOTRAN program for Atkinson and Boore model (1995)

The source spectrum, E, is obtained by the following equations specifying both the shape and the amplitude as a function of the earthquake size:

$$E\left(M_0,f\right)=C{M}_{0}S\left(M_0,f\right)$$E2

$$S\left(M_0,f\right)=S_a\left(M_0,f\right)\times S_b\left(M_0,f\right)$$E3

By adopting the source spectrum model AB95 (Atkinson and Boore model 1995), the above equation for source spectrum is rewritten considering the seismic moment dependence of the above factors S_a in terms of corner frequencies f_a and f_b:

$$E\left(M_0,f\right)=C{M}_0\left\{\frac{1-\epsilon }{1+{\left[\raisebox{1ex}{$f$}\!\left/ \!\raisebox{-1ex}{$f_a$}\right.\right]}^2}+\frac{\epsilon }{1+{\left[\raisebox{1ex}{$f$}\!\left/ \!\raisebox{-1ex}{$f_b$}\right.\right]}^2}\right\}$$E4

where, C is a constant given by

$$C=\frac{<R_{\Theta \Phi }>VF}{4\pi {\rho }_s{\beta }_s{R}_o}$$E5

Here, $$<R_{\Theta \Phi }>$$ accounts for the radiation pattern (≈ 0.55); V represents the partition of total shear wave energy into horizontal components (= 0.707); F accounts the effect of free surface (≈ 2); $${\rho }_s$$ and $${\beta }_s$$ are the density and shear wave velocity of the bedrock; R_o is a reference distance and usually taken as 1 km. The corner frequencies f_a and f_b are obtained from the seismic moment using the following relations

$$\log f_a=2.41-0.533M_0$$E6

$$\log f_b=1.431-0.188M_0$$E7

The Source duration is evaluated as 0.5/f_a.

The path effects are represented by simple functions that describe the geometric spreading function, attenuation (intrinsic and scattering attenuation), and the general increase of duration with distance due to wave propagation and scattering. The simplified path effect, P, is given by the multiplication of the geometrical spreading and Q functions:

$$P=Z\left(R\right)\exp \left\{\frac{-\pi Rf}{Q\left(f\right)C_Q}\right\}$$E8

The relation between distance and geometrical spreading function, Z(R), is given by the following function

$$Z\left(R\right)=\frac{1}{70}\sqrt{\frac{130}{R}}$$E9

and Q(f) is the frequency dependent quality factor which is given by the following equation

$$Q\left(f\right)=100f^{0.8}$$E10

The path duration function of 0.05R is calculated from Atkinson and Boore (1995).

The attenuation or diminution operator D(f) accounts for the path-independent loss of high frequency in the ground motions. A simple multiplicative filter can account for the diminution of the high frequency motions. Here, f_max is 10 Hz. The diminution factor is calculated based on the following equation

$$D\left(f\right)={\left\{1+{\left(\raisebox{1ex}{$f$}\!\left/ \!\raisebox{-1ex}{$10$}\right.\right)}^8\right\}}^{-0.5}$$E11

The particular type of ground motion resulting from the simulation is controlled by the filter I(f). If ground motion is desired, then

$$I=-{\left(2\pi f\right)}^n$$E12

where, I = (-1)^0.5. n = 0, 1, 2 for ground displacement, velocity and acceleration, respectively.

A time domain simulation is carried out to get the actual Fourier amplitude spectrum. A White Gaussian Noise (WGN) is produced and windowed off using a windowing function given below

$$W\left(t;\epsilon ,\eta ,t_{\eta }\right)=a{\left(t/t_{\eta }\right)}^b\exp \left\{-c\left(t/t_{\eta }\right)\right\}$$E13

where,

$$\begin{array}{l}a={\left\{\exp \left(i\right)/e\right\}}^b\\ b=-\left(\epsilon \ln \eta \right)/\left[1+\epsilon \left(\ln \epsilon -1\right)\right]\\ c=b/\epsilon \\ t_{\eta }=f_{Tgm}\times T_{gm}\end{array}$$E14

Boore suggested the values of ε and η to be 0.2 and 0.05, respectively and f_Tgm = 2 based on Saragoni and Hart (1974). The windowed WGN is converted to frequency domain and normalized by its root mean square amplitude. The entire process of obtaining WGN is shown in Figure 2. Then, the ground motion spectrum, shown in Figure 2, is multiplied with the normalized windowed noise to get the Fourier amplitude spectrum as shown in Figure 2.

Figure 2.

Basis of procedure for simulating ground motions using the stochastic method

2.2. Analytical model proposed by Mavroeidis and Papageorgiou (2003)

For near-field strong ground motions, most of the elastic energy arrives coherently in a single, intense, relatively long period pulse at the beginning of record, representing the cumulative effect of almost all the seismic radiation from the fault. The phenomenon is even more pronounced when the direction of slip on the fault plane points toward the site as well. The Mavroeidis and Papageorgiou model (2003) adequately describes the impulsive character of near-faults ground motions both qualitatively and quantitatively. In addition, it can be used to analytically reproduce empirical observations that are based on available near-source records. The input parameters of the model have an unambiguous physical meaning. The proposed analytical model has been calibrated using a large number of actual near-field ground motion records. It successfully simulates the entire set of available near-fault displacement, velocity, and (in many cases) acceleration time histories, as well as the corresponding deformation, velocity, and acceleration response spectra. An “objective” definition of the pulse duration is given based on model input parameters. In addition, Mavroeidis and Papageorgiou (2003) investigate the scaling characteristics of the model parameters with earthquake magnitude. Also, Mavroeidis et al. (2004) derive the Fourier transform of the analytical model and identify the parameters that have the most significant effect on the spectral characteristics of the model. Finally, a simplified (adequate for engineering purposes) method is proposed for the synthesis of near-fault ground motions. The pulse duration (or period), the pulse amplitude, as well as the number and phase of half cycles are the key parameters that define the waveform characteristics of near-fault velocity pulses. Therefore, an analytical model with four parameters in principle should suffice to describe the entire set of velocity pulses generated due to forward directivity or permanent translation effects. Seismologists have used “wavelets” (also referred to as “signals,” “signatures,” or “pulses”), particularly in fields such as seismic filtering, wavelet processing, wave-propagation modelling, and trace inversion (Hubral and Tygel 1989). Although, various wavelets have been proposed in the literature, only a limited number of them are popular and frequently used in practice. In this study, the analytical wavelet signal proposed by Mavroeidis and Papageorgiou (2003) is chosen and expressed by

$$f\left(t\right)=A\frac{1}{2}\left[1+\cos \left(\frac{2\pi f_p}{\gamma }t\right)\right]\cos \left(2\pi f_{p}t+\nu \right)$$E15

This problem is easily resolved by limiting the time interval of the signal as follows

$$-\frac{\gamma }{2f_p}\le t\le \frac{\gamma }{2f_p}$$E16

The period of the harmonic oscillation should be smaller than the period of the envelope represented by the elevated cosine function in order to produce physically acceptable signals; that is,

$$\frac{1}{f_p}<\frac{\gamma }{f_p}\Rightarrow \gamma >1$$E17

The combination of equations (15) to (17) yields the formulation of the proposed analytical model for the near-fault ground velocity pulses:

$$\nu \left(t\right)=\left\{\begin{array}{l}A\frac{1}{2}\left[1+\cos \left(\frac{2\pi f_p}{\gamma }\left(t-t_0\right)\right)\right]\cos \left(2\pi f_p\left(t-t_0\right)+\nu \right),\\ 0,otherwise\end{array}\right\}$$E18

where, parameter A controls the amplitude of the signal; f_p is the frequency of the amplitude-modulated harmonic (or the prevailing frequency of the signal); $\nu$ is the phase of the amplitude-modulated harmonic (i.e., $\nu =0$ and $\nu =\pm \pi /2$ define symmetric and antisymmetric signals, respectively); $\gamma$ is a parameter that defines the oscillatory character (i.e., zero crossings) of the signal (i.e., for small $\gamma$ the signal approaches a deltalike pulse; as $\gamma$ increases, the number of zero crossings increases); and t₀ specifies the epoch of the envelope’s peak. The analytical expressions for the ground acceleration and displacement time histories, compatible with the ground velocity given by equation (18), are

$$a\left(t\right)=\left\{\begin{array}{l}-\frac{A\pi f_p}{\gamma }\left[\begin{array}{l}\sin \left(\frac{2\pi f_p}{\gamma }\left(t-t_0\right)\right)\cos \left[2\pi f_p\left(t-t_0\right)+\nu \right]+\\ \gamma \sin \left[2\pi f_p\left(t-t_0\right)+\nu \right]\left[1+\cos \left(\frac{2\pi f_p}{\gamma }\left(t-t_0\right)\right)\right]\end{array}\right],t_0-\frac{\gamma }{2f_p}\le t\le t_0+\frac{\gamma }{2f_p}with\gamma >1\\ 0,otherwisw\end{array}\right\}$$E19

$$d\left(t\right)=\left\{\begin{array}{l}\frac{A}{4\pi f_p}\left[\begin{array}{l}\sin \left(2\pi f_p\left(t-t_0\right)+\nu \right)+\frac{1}{2}\frac{\gamma }{1-\gamma }\sin \left[\frac{2\pi f_p\left(1-\gamma \right)}{\gamma }\left(t-t_0\right)+\nu \right]+\\ \frac{1}{2}\frac{\gamma }{1+\gamma }\gamma \sin \left[\frac{2\pi f_p\left(1+\gamma \right)}{\gamma }\left(t-t_0\right)+\nu \right]\end{array}\right]+C,t_0-\frac{\gamma }{2f_p}\le t\le t_0+\frac{\gamma }{2f_p}with\gamma >1\\ \frac{A}{4\pi f_p}\frac{1}{\left(1-{\gamma }^2\right)}\sin \left(\nu -\pi \gamma \right)+C,t<t_0-\frac{\gamma }{2f_p}\\ \frac{A}{4\pi f_p}\frac{1}{\left(1-{\gamma }^2\right)}\sin \left(\nu +\pi \gamma \right)+C,t>t_0+\frac{\gamma }{2f_p}\end{array}\right\}$$E20

A parametric study in terms of $\nu$ and $\gamma$ of the normalized (with respect to f_p and A) acceleration, velocity, and displacement pulses can be performed based on the equations presented previously. We define the normalized time variable as

$$\overline{t}=2\pi f_p\left(t-t_0\right)$$E21

Then, the normalized acceleration and displacement time histories can be expressed by rewriting equations (19) and (20), as

$$\overline{a}\left(\overline{t}\right)=\frac{a\left(t\right)}{A{f}_p}=\left\{\begin{array}{l}-\frac{\pi }{\gamma }\left[\sin \left(\frac{\overline{t}}{\gamma }\right)\cos \left(\overline{t}+\nu \right)+\gamma \sin \left(\overline{t}+\nu \right)\left(1+\cos \left(\frac{\overline{t}}{\gamma }\right)\right)\right],-\pi \gamma \le \overline{t}\le \pi \gamma with\gamma >1\\ 0,otherwisw\end{array}\right\}$$E22

$$\overline{d}\left(\overline{t}\right)=\frac{d\left(t\right)}{\left(A/f_p\right)}=\left\{\begin{array}{l}\frac{1}{4\pi }\left[\sin \left(\overline{t}+\nu \right)+\frac{1}{2}\frac{\gamma }{\gamma -1}\sin \left(\frac{\gamma -1}{\gamma }\overline{t}+\nu \right)+\frac{1}{2}\frac{\gamma }{\gamma +1}\sin \left(\frac{\gamma +1}{\gamma }\overline{t}+\nu \right)\right],-\pi \gamma \le \overline{t}\le \pi \gamma with\gamma >1\\ \frac{1}{4\pi }\frac{1}{\left(1-{\gamma }^2\right)}\sin \left(\nu -\pi \gamma \right),\overline{t}<-\pi \gamma \\ \frac{1}{4\pi }\frac{1}{\left(1-{\gamma }^2\right)}\sin \left(\nu +\pi \gamma \right),\overline{t}>\pi \gamma \end{array}\right\}$$E23

Assuming that the duration of the pulse is independent of the source–station distance for stations located within ~10 km from the causative fault, the pulse period and the moment magnitude are related through the following empirical relationship obtained by least-squares fit analysis:

$$\log T_p=-2.2+0.4M_W$$E24

In this section, we propose a very simplified methodology for generating realistic synthetic ground motions that are adequate for engineering analysis and design. We exploit the simple analytical model introduced in the present work to describe the coherent (long-period) component of motion and the stochastic (or engineering) approach to synthesize the incoherent (high-frequency) seismic radiation (for a review of the stochastic approach of ground motion synthesis; see Boore (1983) and Shinozuka (1988)). For the latter component of motion, due to the proximity of the point of observation to the source, it is necessary to use a source model that provides guidance as how to distribute the available seismic moment of the simulated event on the fault plane. Such a source model is the specific barrier model of Papageorgiou and Aki (1983). According to this model, an earthquake is visualized as a sequence of equal-size sub-events uniformly distributed on a rectangular fault plane. At the present time, the proposed mathematical model along with its scaling laws can take into account (with confidence) only for the forward directivity effect. Even though the analytical expression can replicate near-fault ground motion records that manifest the permanent-translation effect as well, the limited number of recordings with permanent translation does not permit the derivation of appropriate scaling laws. Therefore, the proposed analytical model should be utilized with caution for the generation of synthetic long-period ground motions that intend to incorporate the permanent translation effect. In these cases, the permanent offsets of the synthetic displacement time histories should be compatible with the tectonic environment and earthquake magnitude of the simulated event. The proposed methodology is written in MATLAB with the following steps:

1. Select the moment magnitude, M_W, of the potential earthquake and calculate the prevailing frequency, f_P, by f_P=1/T_P. For selected values of the parameters A, $\gamma$ and $\nu$ (or for a suite of values of these three parameters), generate the coherent component of acceleration time history (or a suite of time histories) using equation (19).

2. For the selected fault–station geometry, generate the synthetic acceleration time histories for the moment magnitude, M_W, specified previously, using the specific barrier model.

3. Calculate the Fourier transform of the synthetic acceleration time histories generated in steps 1 and 2.

4. Subtract the Fourier amplitude spectrum of the synthetic time history generated in step 1 from the Fourier amplitude spectrum of the synthetic time history produced in step 2.

5. Construct a synthetic acceleration time history so that (a) its Fourier amplitude spectrum is the difference of the Fourier amplitude spectra calculated in step 4; and (b) its phase coincides with the phase of the Fourier transform of the synthetic time history generated in step 2.

6. Superimpose the time histories generated in steps 1 and 5. The near-source pulse is shifted in time so that the peak of its envelope coincides with the time that the rupture front passes in front of the station.

3. The seismotectonic and seismicity of Tombak region

The Zagros region is one of the most seismically active regions in Iran. The Tombak LNG terminal is located along the Persian Gulf northern coast, south of the Zagros Mountains, which mark the deforming zone separating Arabia (Arabian plate) and Central Iran (Eurasian plate) (Figure 3(a)). Location of Tombak area is presented in Figure 3(b). The massive LNG storage tanks exist in this terminal. These tanks have high importance from engineering and economical point of view so seismic loads should be considered in their analysis and design. The relevant codes of LNG storage containers emphasize that a comprehensive seismic hazard investigation should be conducted for regional seismicity and earthquake events of known near-fault.

Figure 3.

Location map: (a) Zagros folded zone, and (b) Tombak area

West of the Makran coast, where oceanic crust is subducting beneath Eurasia, the collision of the Arabian shield with Iran has uplifted the Zagros Mountains. The Zagros Mountains belt represents the early stage of a continental collision between the Arabian plate and the central Iran continental blocks. The Zagros Mountains are a seismically active region. Seismicity is restricted to the region between the Main Zagros Thrust and the Persian Gulf. Strong earthquakes are thought to occur on blind active thrust faults, which do not reach the surface. Fault plane solutions of these earthquakes indicate displacement mainly on low to high-angle reverse faults at depth of 6-12 km in the uppermost part of the basement. Most of the earthquakes for the region have generally M = 5.0 to 6.5, and have originated on sources beneath the decollement (Berberian 1995). Subduction on the main Zagros thrust has now ceased and it is seismically inactive (Ni and Barazangi 1986) except for the northern Zagros, where the surface trace of the thrust has been reactivated as right-slip main recent fault. The Zagros active fold-thrust belt lies on the north-eastern margin of the Arabian plate, on Precambrian (Pan-African) basement. It is composed of Cambrian to Neogene's folded series and is the result of five major tectonic events (Berberian and King 1981; Berberian 1983). The Zagros fold-thrust belt is composed of five units. The folds are parallel to the thrust faults. The axial part of the folds, striking NW SE, appears as broad asymmetrical folds with axial planes dipping to the NE and North. Their north-eastern limbs gently dip (20°) to the NW whereas their south-western limbs are steeper (40°) to the SE reaching 60 to 80°down slope and in some cases are nearly vertical, overturned or thrusted. The Main Zagros Thrust Fault (MZTF) indicates a fundamental change in sedimentary and structural evolution and seismicity. It marks the geosuture between the two colliding plates of the Eurasia and the Arabia. The global zone taken into account lies between 32°N and 26°N in latitude and 50°E and 58°E in longitude.

4. Estimation of the model parameters

The source and earthquake parameters have been obtained from the pervious seismic hazard study which has been conducted for Tombak area. The most of the models based on the stochastic method are fundamentally point-source models. Although it is true that near and intermediate-field terms are lacking, in most applications the frequencies are high enough that the far-field terms dominate, even if the site is near the fault. Furthermore, the effects of a finite-fault averaged over a number of sites distributed around the fault (to average over radiation pattern and directivity effects) can be captured in several ways: 1) using the closest distance to faulting as the source-to-site distance; 2) using a two-corner source spectrum; 3) allowing the geometrical spreading to be magnitude dependent. The material properties described by density ρ, and shear wave velocity β, are estimated to be 2.8 gr/cm³ and 3.5 km/sec, respectively. All parameters used for simulation are summarized in Table (1).

Parameters	Values
${\rho }_s,{\beta }_s,V,<R_{\Theta \Phi }>,F,R_0$	2.8, 3.5, 0.707, 0.55, 2.0, 1.0
Geometrical spreading (including factors to insure continuity of function)	r<40 km: 1/r; r$\ge$40 km: (1/40)(40/r)^0.5
Q, cQ	180f^0.45, 3.5 km/s
Source duration	0.5/fa
Path duration	0.05 R
Site amplification	Boore and Joyner (1997) generic rock
Site diminution parameters (fmax ,$\kappa$)	100.0, 0.03

Table 1.

Model parameters

In the methodology of Beresnev and Atkinson (1997,1998), modelling of finite source requires information of the orientation and dimensions of fault plane, as well as information of the dimensions of sub-faults and the location of hypocenter. The trends of epicenteral and hypocenteral distribution are in accordance with the strike and dip angle of the focal mechanism (strike, dip, slip) = (175, 85, 153) of the mainshock (Yamanaka 2003). The source dimension is therefore roughly estimated to be 20 km x 16 km (Yamanaka 2003).

Source parameters can be classified into two types (Irikura 2000): global source parameters and local source parameters. They represent different features of the fault source and are determined by different methods. The global source parameters characterize the macro feature of the entire source area and include spatial orientation of fault (location, attitude, buried depth), fault size (length, width, area), and both average slip and average rupture velocity on the fault plane. In the global source parameters both the slip type and spatial orientation are determined by seismogeology investigation and geophysical exploration; while the moment magnitude of the scenario earthquake caused by an active fault is estimated from its seismic hazard assessment. Fault size and average slip on the fault plane are also estimated by seismic scaling laws. In this study, the information for generating near-field strong ground motion such as magnitude related to each return period and the epicentral and hypocentral distances for stochastic method have been extracted from the seismogeology investigation that have been presented in Table (2). Table (3) lists the basic parameters used in the strong ground motion predictions.

Zone	Seismic source	Tombak
Magnitude	Epicentral distance (km)	Depth (km)	Hypocentral distance (km)
ZFF-C	OBE (475 years)	6.5	5	11	12
SSE (5000 years)	7.0	5	14	15

Table 2.

Parameter values of finite fault source model

Parameters	Values
Fault orientation	Strike 122°, Dip 40°
Fault dimensions along strike and dip (km)	28 × 16
Burial depth of upper limit of the fault (km)	5.0
Moment magnitude (MW )	6.7
Subfault dimensions along strike and dip (km)	1 × 1
Stress drop (bar)	50
Q(f)	150f ^0.5
Geometrical spreading	1/R
Windowing function	Saragoni-Hart
Kappa	0.05
Crustal shear wave velocity (km/s)	3.7
Rupture velocity (km/s)	0.8 × shear wave velocity
Crustal density (g/cm^3 )	2.8

Table 3.

Basic parameters used in the strong ground motion predictions

Using the magnitude for each return period, the duration of pulse is defined for both levels of earthquakes. The other parameters for the pulse are extracted from Mavroeidis and Papageorgiou (2003). The data has been obtained by the calibration procedure of actual near-fault strong ground motion records. The chosen parameter values have been summarized in Table (4).

ω	Return period (year)	Amplitude of pulseA (cm/s)	Oscillatory character of the signal $\gamma$	Phase of Pulse$\nu$ (degree)	Time shiftt0 (second)
6.5	475	50	1.5	134°	5
7.0	5000	90	2	100°	6.4

Table 4.

Input parameters for defining long period pulse

5. Results and discussion

A simple and powerful method for simulating ground motions is based on the assumption that the amplitude of ground motion at a site can be specified in a deterministic way, with a random phase spectrum modified such that the motion is distributed over a duration related to the earthquake magnitude and to distance from the source. This stochastic method is particularly useful for simulating the higher-frequency ground motions, and it is used to predict ground motions for regions of the world in which recordings of motion from damaging earthquakes are not available. This simple method has been successful in matching a variety of ground motion measures for earthquakes with seismic moments spanning more than 12 orders of magnitude. SMSIM (StochasticModel SIMulation or Strong Motion SIMulation) is a set of programs for simulating ground motions based on the stochastic method. Programs are included both for time-domain and for random vibration simulations. In addition, programs are included to produce Fourier amplitude spectra for the models used in the simulations and to convert shear velocity versus depth into frequency-dependent amplification. The necessary parameters in these models are distinguished for defining the theoretical relationships. In this study, the near-field strong motion time histories, obtained from the stochastic method, are presented for both levels of earthquakes (475y and 5000y) as shown in Figure 4. The long period pulse has been calculated and presented in Figure 5 for both return periods.

In this stage, the pulse acceleration superimposes to the synthetic acceleration time history. The near-source pulse is shifted in time so that the peak of its envelope coincides with the time of rupture front of station. The final acceleration time histories for both levels of earthquakes are shown in Figure 6. The response spectra, obtained from simulated strong ground motion analysis, for 5% damping ratio are shown in Figure 7 for return periods 475 and 5000 years. As it is seen in Figures 7(a) and 7(b), the response spectrum has a sudden increasing for period ranges of 1 to 4 and 2 to 6 for 475y and 5000y, respectively. Therefore, structures which their periods settle in these ranges are influenced from near-field due to the long period pulse ground motion. In Figure 8, the smoothed response spectrum, obtained for return period of 475 years, is compared with IBC 2000 and Standard No. 2800. This figure shows that the response spectra are close to each other in the period range of 0 to 1 second. When the period exceeds than 1 second, the effect of long period pulse becomes apparent in the response spectra.

Figure 4.

Synthetic acceleration time histories for return periods: (a) 475 years, and (b) 5000 years

Figure 5.

Long period pulses for return periods: (a) 475 years, and (b) 5000 years

Figure 6.

Final acceleration time histories for return periods: (a) 475 years, and (b) 5000 years

Figure 7.

Response spectra of final acceleration time histories for return periods: (a) 475 years, and (b) 5000 years

Figure 8.

Comparison between response spectra (Blue Curve: 475y, Red Curve: Standard No. 2800, Gray Curve: IBC 2000)

6. Site response analysis

Soil investigation and borings are carried out for the detail design of storage tank foundations and spill basin structures of NIOC LNG Project. The site is located at the Persian Gulf coast in Tombak region, approximately 60 km away from Assaluyeh city, in Bushehr province. Standard penetration test has been running since initial stages of drilling operation. SPT tests have been performed approximately in each 1.5 m advance of drilling. Thus, there is a full set of data covering whole area with SPT results. Seismic tests of down-hole are done for full depth (80 m). Shear and compression waves velocities (Vs & Vp) are determined, accordingly. The results are presented in Table (5). The longitudinal wave velocity increases due to water table. The water table in borehole is approximately 9.0 m below the ground. Generally, the shear wave velocity increases versus depth due to an increase of soil density. According to Standard No. 2800, the shear wave velocity more than 760 m/s is assigned as rock; therefore the seismic bed rock is located at 10 m below the ground. The mentioned standard is used to classify the soil type which is type II and class C in this project. One-dimensional ground response analysis of the site is carried out by the equivalent linear approach using SHAKE 91 program (Idriss et al. 1992).

Depth	Density	Vp	Vs	E	G	K	$\upsilon$
m	gr/cm^3	m/s	m/s	MPa	MPa	MPa
0-1.3	2.1	570	330	571	229	377	0.25
1.3-4	2.1	850	500	1297	525	817	0.24
4-5.5	2.1	1040	600	1891	756	1263	0.25
5.5-10	2.1	1500	700	2801	1029	3353	0.36
10-20	2.25	1520	760	3466	1300	3466	0.33
20-30	2.3	1600	800	3925	1472	3925	0.33
30-42	2.15	1600	835	3936	1499	3505	0.31
42-52	1.96	1650	870	3879	1484	3358	0.31
52-62	2.2	1600	820	3911	1479	3660	0.32
62-78	2.17	1650	900	4529	1758	3564	0.29

Table 5.

Soil properties

According to the geotechnical site investigation, the soil type and the thickness of each layer are defined. From depth 0.0 m to 8.0 m, there are very diverse layers including boulder, gravel with some silty sand or sandy silt. SPT tests for these layers are refused because of coarse size of grains. From 8.0 m to 41.6 m, a very dense light brown sandy gravel with some cobbles is observed. SPT values are very high for this layer and some trace silt is seen from depth 32.0 m to 32.5 m. After this layer, a very dense light brown silty sand with some gravel is exist with 2.0 m thick. From 43.5 m to 50.0 m, we can see a very dense light brown sandy gravel with some cobbles and trace silt and clay. This layer changes to very dense gray silty sand with trace boulder (at depth: 50.0 m to 52.0 m). From 52.0 m to 52.6 m, borehole drilling machine interfaces to a piece of rock presenting by "boulder" phrase in borehole log. From 52.6 m to 54.0 m, borehole log shows a very dense to medium dense of gray silty sand. This layer is become very dense and its color changes to light brown from depth 54.0 m to 61.0 m. From depth 61.0 m to 64.0 m, some gravel and trace silt is added to last mentioned layer, so we can see a silty sand with some gravel and trace silt. From 64.0 m to 73.6 m, borehole log shows a very dense light brown sandy gravel with trace cobbles and clay. At the end of borehole, to depth 80.0 m, borehole log shows a very dense light brown silty sand with trace gravel and clay. Considering the results, particularly the seismic down-hole and SPT data, it is found that the average shear wave velocity equals approximately to 650 m/s and the SPT number is more than 50. For site specific response analysis, it is recommended to use shear wave velocity measured in the field. The average shear wave velocity of soil within 30 m was found to be around 650 m/s. The simulated ground motion is used as the input motion for ground response analysis. The results of ground response analysis are presented in Figure 9. The variation of maximum acceleration with depth as shown in Figure 9 indicates that the increase of PGA at the surface is about 1.06 and 1.11 times higher than those on the bed rock for return periods of 475 years and 5000 years, respectively. The response spectra obtained from ground response analysis for 5% damping ratio are shown in Figure 10.

Figure 9.

Acceleration time histories: (a) simulated on bedrock for return period = 475 years, (b) obtained on ground surface for return period = 475 years, and (c) simulated on bedrock for return period = 5000 years, (d) obtained on ground surface for return period = 5000 years

Figure 10.

Spectral accelerations obtained by response analyses for return periods: (a) 475 years, and (b) 5000 years

7. Conclusions

The present study has been focused on simulating near-field strong ground motions using hybrid method for Tombak area in south-eastern part of Iran. The simulation of ground motion has been carried out using the stochastic method proposed by Boore (2003). Afterwards, the analytical model proposed by Mavroeidis and Papageorgiou (2003) is applied to consider the impulsive character of near-fault ground motion. Then, the response of mentioned site under simulated ground motion has been studied by conducting one dimensional ground response analysis. The results can be used for estimating the probable ground motion acceleration time-histories to be used in the hazard analysis of specific sites in the region under study, particularly for performance analysis of exiting structures. The ability of this hybrid method in simulating strong motions is also shown in this study. The simulation parameters, obtained in this study, can be used to asses the strong-motion level at a much larger number of sites, where no record is available, to investigate how different characteristics of motion affect the damage distribution in the Tombak region. Based on the above study the following conclusions are derived:

• The maximum PGAs on bedrock at the Tombak site are established as 0.32 g and 0.35 g for return periods of 475 and 5000 years, respectively.

• Ground response analysis indicates small increases of PGAs at the surface which are about 1.06 and 1.11 times higher than those on the bedrock for return periods of 475 and 5000 years, respectively.

• The response spectra obtained from the analyses indicate that the effect of long period pulse appears in the period ranges of 1 to 4 and 2 to 6 for 475y and 5000y, respectively, which are the typical ranges of period for LNG storage tank structures constructed in the Tomabk region.

• The obtained results demonstrate that the employed models adequately describe the nature of the impulsive near-fault ground motions both qualitatively and quantitatively.

• The proposed procedure mentioned in the present work will facilitate the study of the elastic and inelastic response of structures subjected to near-source seismic excitations within the considered region.