Site-specific strong ground motion prediction using 2.5-D modelling

Summary

An algorithm was developed using the 2.5-D elastodynamic wave equation, based on the displacement–stress relation. One of the most significant advantages of the 2.5-D simulation is that the 3-D radiation pattern can be generated using double-couple point shear-dislocation sources in the 2-D numerical grid. A parsimonious staggered grid scheme was adopted instead of the standard staggered grid scheme, since this is the only scheme suitable for computing the dislocation. This new 2.5-D numerical modelling avoids the extensive computational cost of 3-D modelling. The significance of this exercise is that it makes it possible to simulate the strong ground motion (SGM), taking into account the energy released, 3-D radiation pattern, path effects and local site conditions at any location around the epicentre. The slowness vector (p_y) was used in the supersonic region for each layer, so that all the components of the inertia coefficient are positive.

The double-couple point shear-dislocation source was implemented in the numerical grid using the moment tensor components as the body-force couples. The moment per unit volume was used in both the 3-D and 2.5-D modelling. A good agreement in the 3-D and 2.5-D responses for different grid sizes was obtained when the moment per unit volume was further reduced by a factor equal to the finite-difference grid size in the case of the 2.5-D modelling. The components of the radiation pattern were computed in the xz-plane using 3-D and 2.5-D algorithms for various focal mechanisms, and the results were in good agreement. A comparative study of the amplitude behaviour of the 3-D and 2.5-D wavefronts in a layered medium reveals the spatial and temporal damped nature of the 2.5-D elastodynamic wave equation.

3-D and 2.5-D simulated responses at a site using a different strike direction reveal that strong ground motion (SGM) can be predicted just by rotating the strike of the fault counter-clockwise by the same amount as the azimuth of the site with respect to the epicentre. This adjustment is necessary since the response is computed keeping the epicentre, focus and the desired site in the same xz-plane, with the x-axis pointing in the north direction.

2.5-D modelling, 3-D radiation pattern, moment tensor source formulation, strong ground motion

Introduction

The scale of socio-economic damage caused by an earthquake depends to a great extent on the SGM characteristics. The ground motion characteristics depend very much on the radiation pattern and the local site condition. Therefore the site-specific ground motion prediction taking into account the energy released, 3-D radiation pattern, path effect and local site conditions is of prime importance in earthquake engineering. Furthermore, the development of 2.5-D modelling is justified since the available 3-D algorithms are limited to very low frequencies in modelling the anisotropic structures, owing to the large amount of memory and computational time required. The configuration for 2.5-D wave simulation is achieved when the medium properties vary only in two dimensions and are assumed constant in the third direction, and the source–receivers are confined to the same plane.

Bleistein (1986) described a mathematical formulation for a 2.5-D ray-based acoustic wave equation for constant- and variable-density media. A different approach was followed by Liner (1991) and Stockwell (1995), who developed the 2.5-D acoustic wave equations using the 3-D Green's function. Williamson & Pratt (1995) proposed a filter operator to determine the 2.5-D response by simply applying the filter to the 2-D response. Pedersen et al.(1994) developed a 2.5–D indirect boundary element method based on moving Green's functions to study the 3-D scattering effects using 2-D parameters. Narayan (1998, 1999) developed 2.5-D algorithms for acoustic wave propagation using constant- and variable-density media and reported that the last term of Liner's (1991) and Williamson & Pratt's (1995) equations are superfluous and that the 2.5-D computation time is only 12–15 per cent more than that for the 2-D computation. Furumura & Takenaka (1996) developed a 2.5-D formulation for the pseudo-spectral point source excitation. Takenaka & Kennett (1996a) developed a 2.5-D time-domain elastodynamic wave equation for plane wave incidence, as well as an extension of this equation for general anisotropic media using the radon transform (Takenaka & Kennett 1996b).

In the past, most SGM simulations were done using explosive sources (Virieux 1986; Kim et al. 1995), and numerical implementations of double-couple point sources were never completely presented, except by Vidale & Helmberger (1987), Coutant et al. (1995), Graves (1996) and Pitarka (1999). Vidale & Helmberger (1987) used the Alterman & Karal (1968) procedure to generate a double-couple source. The source implementation based on the moment tensor formulation is simpler and contains the complete source description (Aki & Richards 1980; Pitarka 1999).

3-D and 2.5-D elastodynamic wave equations were approximated using the parsimonious staggered grid approximation (Ohminato & Chouet 1997), and corresponding algorithms were developed. The double-couple point shear-dislocation source with a particular focal mechanism based on the moment tensor formulation (Aki & Richards 1980; Pitarka 1999) was implemented into the numerical grid. The radial, transverse and vertical components of the radiation pattern for various focal mechanisms were computed using 3-D and 2.5-D algorithms. The components of the ground displacement were computed using the moment per unit volume by varying the effective volume of the finite-difference cell. This paper gives the site-specific SGM prediction using a 2.5-D simulation taking into account the 3-D radiation pattern and the energy release. A vacuum formulation was used at the free surface for stability. A sponge boundary condition was implemented on the model edges to avoid edge reflections (Israeli & Orszag 1981).

The 2.5-D Elastodynamic Wave Equation

The 3-D elastodynamic wave equation for a linear, isotropic elastic medium without a source is as follows:  

(1)

and the stress–strain relation is  

(2)

where (u, v, w) are the displacement components; (τ_xx, τ_yy, τ_zz, τ_xz, τ_xy, τ_zy) are the stress components; ρ is the density; λ and µ are Lame's coefficients; and the symbols ∂_x, ∂_y, ∂_z and ∂_tt are shorthand representations of the differential operators ∂/∂_x, ∂/∂_y, ∂/∂_z and ∂²/∂t², respectively.

Takenaka & Kennett (1996a) derived a 2.5-D elastodynamic wave equation for a 3-D wavefield in a medium varying in two dimensions only, assuming the medium to be invariant in the y-direction. They took advantage of the invariant nature of the horizontal slowness vector along the invariant axis of the medium. They replaced the space derivatives (with respect to y) of the displacement and stress components of eqs (1) and (2) by the time derivatives of the corresponding displacement and stress components multiplied by the slowness vector in the y-direction (see Takenaka & Kennett 1996a for details). The following 2.5-D time-domain wave equation was developed using the above replacements:  

(3)

and the stress–strain relation is  

(4)

where p_y is the slowness vector in the y-direction. It can be noted that all the variables in eqs (3) and (4) are real-valued, so the propagation characteristics can be efficiently calculated using 2-D time-domain numerical techniques.

Numerical Approximation

The 2.5-D elastodynamic wave equation developed by Takenaka & Kennett (1996a) was approximated using the parsimonious staggered grid scheme (Ohminato & Chouet 1997). This scheme is stable for the full range of Poisson's ratio and is free from the spatial derivative of the elastic parameters. The presence of time derivatives of stress in the 2.5-D wave equation requires the stress components at two times, so there is no parsimony in storage requirement in the present study. However, the parsimonious scheme based on the displacement–stress relation is more suitable for the computation of dislocation than that based on the velocity–stress relation (Coutant et al. 1995).

The time and spatial derivatives of eqs (3) and (4) were replaced by the second-order staggered grid differencing operator:  

(5)

 

(6)

The slowness vector (p_y) has been retained in the supersonic region so that all the components of the inertia coefficient are positive (see eq. 21 of Takenaka & Kennett 1996a). The stresses at time ‘n’ used in eq. (5) were computed using the differenced stress–strain relation given in eq. (6), where Δt is the time step, and Δx and Δz are the grid intervals in the x- and z-directions, respectively. ‘l’ and ‘m’ are the grid indices in the x- and z-directions, and ‘n’ is the time index.

Free Surface and Absorbing Boundary

Free-surface boundary conditions often require careful consideration in numerical schemes from the point of view of stability and accuracy. In the present study, a vacuum formulation was adopted in which α, β and ρ → 0 in the region above the free surface (Graves 1996). Graves (1996) reported that this formulation is stable for second-order schemes and unstable for fourth-order schemes. The vacuum formulation is very useful in modelling the effects of surface topography, since it can be implemented with the same difference equation as used in the interior of the model. A sponge boundary condition, developed by Israeli & Orszag (1981), was implemented on the model edges to avoid edge reflections.

Moment Tensor Source Implementation

The source implementation is based on the moment tensor source formulation that uses stress tensor components at the source location (Pitarka 1999). The equivalent body forces are added to each component of the stress tensor at every time step, in contrast to Graves (1996), who added body forces to each velocity component of the elastodynamic wave equation. 3-D stress tensor components for a source location at the grid point (l, m, k) are

 

where M_xx(t), M_yy(t), M_zz(t), M_xy(t), M_xz(t) and M_zy(t) are the moment tensor components. V is the volume of a finite-difference cell, rather than the volume of a material cell, since the effective finite-difference grid size is half of the material grid size (Ohminato & Chouet 1997). Moment per unit volume was used in the 2.5-D modelling instead of moment per unit area.

A simple homogeneous model was simulated using 3-D and 2.5-D algorithms to prove the accuracy and the validity of the 2.5-D elastodynamic wave equation, using various values of moment. Moment was varied in the simulation by varying the finite-difference grid size (cell volume) for a fixed value of moment per unit volume. The slowness vector (p_y) was kept in the supersonic region (p_y = 0.016/v_yi; v_yi is the velocity of P and S waves in the ith layer) so that all the components of the inertia coefficient are positive. The second derivative of convolution of a polynomial window [{1 − (τ −1)²}³] and the Gaussian function [exp{–α(T– T₀)²}] with a dominant period T₀ was used as a source function. τ is equal to T/T_0. The model parameters, namely Lamé's parameters (λ and µ), density and dominant frequency, were taken as 10 GPa, 10 Gpa, 2.5 g cc⁻¹ and 1.5 Hz, respectively. Radial, transverse and vertical components of the ground displacement were computed at an offset of 6.3 km using a double-couple point source with a focal mechanism with a strike of 90°, a dip of 45° and a rake of 90° and at a depth of 6.5 km.

Table 1.

Maximum ground displacement caused by a source with a focal mechanism of strike=90°, dip=45° and rake=90° over a half-space.

3-D and 2.5-D responses were computed using various grid sizes (90m, 100m and 110 m) and a fixed amount of moment per unit volume (Table 1). Therefore the total moment used in each case was different, depending on the effective volume of the finite-difference cell. Fig. 1 shows the 2.5-D and 3-D radial, transverse and vertical components of the amplitude responses for the different grid sizes. This figure shows that the variation in ground displacement is in agreement with the input moment used. The maximum ground displacements in the 3-D and 2.5-D radial, transverse and vertical responses are noted in Table 1, for the different grid sizes. The 2.5-D amplitude responses corroborate the 3-D responses. It was concluded on the basis of the results of the numerical experiments that the moment per unit volume should be used in 2.5-D modelling, and it should be scaled down by the finite-difference grid size in order to obtain an amplitude response comparable with the 3-D amplitude. However, this conclusion requires further detailed study. It proves the validity and accuracy of the 2.5-D wave equation.

Figure 1.

3-D and 2.5-D radial, transverse and vertical responses of a half-space model for various grid sizes.

3-D and 2.5-D Double-Couple Radiation Pattern

The radiation patterns for various source mechanisms were generated using the relation between the fault parameter strike, dip, rake and the moment tensor components (Aki & Richards 1980) in a Cartesian co-ordinate system. 3-D and 2.5-D snapshots were computed for individual components of the radiation patterns using various source mechanisms for a homogeneous model with the same Lamé parameters of 10 GPa and a density of 2.5 g cc⁻¹. The model size was 15 × 15 × 3 km, and the source was kept in the centre at a depth of 1.0 km. The snapshots were computed after 2.0 s. The grid size, time step and dominant frequency for the modelling parameters were 100 m, 0.01 s and 2.0 Hz, respectively. Figs 2(a) and (b) show the radial and vertical components of the 3-D radiation pattern on the surface for a double-couple point shear-dislocation source with focal mechanism with a strike of 90°, a dip of 90° and a rake of 90°. There are four shear lobes in the radial component and only two lobes in the vertical component; and there are only two compressional lobes in both the radial and vertical components of the 3-D radiation patterns at the surface. It should be noted that the maximum positive amplitude in the snapshots is shown in blue, and negative amplitude is shown in yellow.

3-D and 2.5-D radiation patterns were computed using the same source (strike = 90°, dip = 90° and rake = 90°) in the xz-plane, with the x-axis pointing in the north direction. The model size was 10 × 7 × 10 km, and the source was kept at 5.0 km depth. The snapshots were taken at 1.5 s in the xz-plane, lying in the centre of the model. The other computational parameters were the same as in the previous case. Figs 2(c) and (d) show the 3-D radial and the vertical components of the radiation pattern. The radial and vertical components of 2.5-D radiation pattern are shown in Figs 2(e) and (f), respectively. A comparison of the 3-D and 2.5-D snapshots reveals that the four-lobed 2.5-D radiation pattern is similar to the 3-D radiation pattern.

Similarly, Fig. 3 shows the 3-D and 2.5-D radial, transverse and vertical components of the radiation pattern in the xz-plane for a point source with a focal mechanism with a strike of 90°, a dip of 45° and a rake of 90°. Again, the compressional and shear lobes in the individual components of the 2.5-D radiation pattern match with the corresponding 3-D radiation pattern. Analysis of Figs 2 and 3 reveals that the compressional and shear lobes in the various components of the 2.5-D and 3-D radiation patterns are identical. The similarity of the 2.5-D and 3-D radiation patterns demonstrates that the 3-D radiation pattern can be generated using the 2.5-D wave equation in a 2-D environment.

Site-Specific SGM Prediction

Recent advances in the performance of computers have brought 3-D modelling of seismic wave propagation just within reach. 3-D computations are still very expensive to perform, however, due to the large memory requirements. In the following paragraphs, it is shown how 2.5-D modelling can be used to predict the site-specific ground motion, taking into account the 3-D radiation pattern. For this purpose, two simple cases are considered. In the first case, the ground response is computed at a site having 45.0° azimuth with respect to the epicentre (Fig. 4b) and is caused by a source with a focal mechanism with a dip of 90.0°, a rake of 0.0° and a strike of 45.0°. First, the 3-D vertical radiation patterns are computed for sources with focal mechanisms with a dip of 90.0°, a rake of 0.0° and a strike of 45.0°, and a dip of 90.0°, a rake of 0.0° and a strike of 0.0°, using a homogeneous model of 15 × 15 × 3 km. Figs 5(a) and (b) show the radiation patterns for the above two cases, respectively. These radiation patterns are identical except for the fact that the second one is rotated 45° counter-clockwise.

Similarly, in the second case, the response is computed at a site having 60° azimuth with respect to the epicentre using a source with a focal mechanism with a dip of 90.0°, a rake of 0.0° and a strike of 30.0°. The radiation patterns shown in Figs 5(c) and (d) are computed using sources with focal mechanisms with a dip of 90.0°, a rake of 0.0° and a strike of 30.0°, and a dip of 90.0°, a rake of 0.0° and a strike of 330.0°, respectively. These are identical except for the fact that the second one is rotated 60.0° counter-clockwise. The analysis of these two cases reveals that ground motion can be predicted at any site around the epicentre with the help of 2.5-D modelling, just by decreasing the strike of the fault by the same amount as the azimuth of the site with respect to the epicentre.

Furthermore, to support the above conclusion, 2.5-D and 3-D responses of a model (Fig. 4a) were computed corresponding to the first case. The epicentre, focus and the receivers were kept in the xz-plane in both the 3-D and 2.5-D cases, with the x-axis pointing in the north direction. The model dimensions, material properties and source–receiver positions are shown in Fig. 4(a). The grid size, time step and dominant frequency were taken as 100.0 m, 0.01 s and 1.5 Hz, respectively. The responses were computed at nine equidistant receivers, 1.0 km apart. Figs 6(a) and (b) show the 2.5-D and 3-D responses for a point source with parameters dip = 90.0°, rake = 0.0° and strike = 45.0° at 7.0 km depth. Similarly, Figs 6(c) and (d) show the 2.5-D and 3-D responses for a point source having parameters dip = 90.0°, rake = 0.0° and strike = 0.0°. The computed results show that the 2.5-D responses are very similar to the corresponding 3-D responses.

It is concluded that site-specific strong ground motion can be predicted, taking into account the 3-D radiation pattern, just by rotating the strike of the fault counter-clockwise by the same amount as the azimuth of the site with respect to the epicentre. This conclusion was drawn on the basis of the detailed analysis of 3-D radiation patterns and the 3-D and 2.5D responses. This adjustment is necessary since the 2.5-D response can only be computed in the xz-plane, with the x-axis pointing in the north direction, and we have to keep the epicentre, focus and desired site in the same xz-plane.

Figure 2.

3-D radiation pattern (radial and vertical components) on the surface; and 3-D and 2.5-D radial and vertical components of the radiation pattern in the xz-plane for a point source with a focal mechanism of dip=90.°, rake=90.° and strike=90.°. Figs 2, 3 and 9 may be viewed in colour in the online version of the journal (http://www.blackwell-synergy.com).

Figure 3.

3-D and 2.5-D radial, transverse and vertical components of the radiation pattern in the xz-plane for a point source with a focal mechanism of dip = 45.0°, rake = 90.0° and strike = 90.0°.

Figure 4.

(a) A homogeneous half-space model and (b) location of desired site.

Figure 5.

3-D radiation pattern (vertical components) for (a) dip = 90.0°, rake = 0.0° and strike = 45.0°, (b) dip = 90.0°, rake = 0.0° and strike = 0.0°,(c) dip = 90.0°, rake = 0.0° and strike = 30.0°, and (d) dip = 90.0°, rake = 0.0° and strike = 330.0°

Figure 6.

(a) 2.5-D and (b) 3-D responses of the model (Fig. 4a) using a source with dip = 90.0°, rake = 0.0° and strike = 45.0°. (c) 2.5-D and (d) 3-D responses of the model (Fig. 4a) using a source with dip = 90.0°, rake = 0.0° and strike = 0.0°. (Note: the number to the right of each trace denotes the maximum amplitude in centimetres.)

5-D and 3-D Simulation of a Layered Model

2.5-D and 3-D responses of a five-layer geological model (Fig. 7) were computed using a double-couple point shear-dislocation source with a focal mechanism with a dip of 45.0°, a rake of 90.0° and a strike of 0.0°, for comparative study and validity of 2.5-D modelling. The model dimensions and source–receiver configuration are shown in Fig. 7, and the model parameters are listed in Table 2. The grid size, time step and dominant frequency were taken as 100.0 m, 0.005 s and 1.0 Hz, respectively. The distance between two consecutive receivers was 1.0 km, and the focal depth was 8.0 km. The moment per unit volume used in the simulation was 10¹¹ N m⁻².

Figs 8(a) and (b) show, respectively, the 3-D and 2.5-D ground displacement (cm) in the radial direction. Comparison of these two figures reveals that the ground displacement pattern is similar in both cases: it is low at the epicentre and then increases out to an offset of 4.0 km, after which it decreases. Furthermore, the generation and character of the multiples are also similar. Figs 8(c) and (d) show the 3-D and 2.5-D ground displacements in the transverse direction, respectively. The ground displacement in the transverse direction is much smaller than the displacement in the radial and the vertical directions.

Figure 7.

Layered geological model and the source–receiver configuration.

3-D and 2.5-D ground displacements in the vertical direction are shown in Figs 8(e) and (f), respectively. The ground displacement caused by the compressional wave is larger than that caused by the shear wave in both cases. Furthermore, the ground displacement caused by the compressional wave is larger in the 2.5-D responses in all cases as compared to the 3-D responses. The amplitudes of the primary P and S waves decrease with epicentral distance, but the amplitude of the mode-converted S wave increases with offset. Comparison of the 2.5-D responses with the 3-D responses reveals the accuracy of the 2.5-D modelling in predicting the various seismic phases and the amplitude response of a complicated geological model.

Table 2.

Parameters and the layer thickness of a geological model.

Figure 8.

3-D and 2.5-D responses of the model (Fig. 7) using a point source with a focal mechanism of dip = 90.0°, rake=0.0° and strike=0.0°. (Note: the number to the right of each trace denotes the maximum amplitude in centimetres.)

3-D and 2.5-D Amplitude Behaviour

3-D and 2.5-D snapshots were computed in the xz-plane at various times for a layered geological model (Fig. 7) for comparative study of the amplitude behaviour at various discontinuities. The model parameters, and the source and its position are the same as in the previous case. Fig. 9 shows the 3-D and 2.5-D snapshots of the radial component at 3.0, 5.0, 7.0 and 9.0 s. The 3-D and 2.5-D snapshots at time 3.0 s show a shear lobe and partially two compressional lobes. Tremendous amplitude amplification in the surficial layer is visible in both cases at time 5.0 s. The up-going primary S phase and the reflected Pp phase from the surface are clearly visible. Reflected Pp and Ss phases from the surface are also clearly visible in the snapshots at time 7.0 s. Similarly, 3-D and 2.5-D snapshots at time 9.0 s depict the reflected Pp and Ss phases, multiples and the new phases in the surficial layer. Comparative study of the snapshots reveals that the behaviour of wavefronts at the various discontinuities encountered in the path is more or less similar in the two cases. The decrease and increase of the wavelength in individual layers, depending on the velocity of the layer, can be easily inferred from the 3-D and 2.5-D snapshots.

Figure 9.

3-D and 2.5-D snapshots (radial component) at various times for a model (Fig. 7) with a point source with a focal mechanism of dip = 90.0°, rake=0.0° and strike=0.0°.

Figure 10.

Variation of the maximum amplitude (m) in the 3-D and 2.5-D wavefronts with time.

Figs 10(a) and (b) show the maximum negative and positive amplitudes present in the snapshots of the 3-D and 2.5-D radial component at various times. The analysis reveals that the difference in the 2.5-D and 3-D amplitude responses is very small except at the impedance contrasts. There is a rapid decrease in the amplitude of the wavefronts up to 3.0 s, and then amplification starts as the wave enters the second layer. The increase in the amplitude of the wavefront is due to the amplitude amplification caused by the decrease of shearing strength in the overlying layer. The amplitude amplification at time 6.0 s is due to a sudden decrease of shearing strength in the surficial layer, as well as possible interference of the outgoing and reflected wavefronts. Comparison of Figs 9(a) and (b) shows that the amplification pattern is different for the 2.5-D and 3-D cases. It can be easily inferred that the crest is more amplified in the 2.5-D simulation, and that the trough is more amplified in the 3-D simulation. The maximum positive and negative amplitudes in the 3-D and 2.5-D snapshots at various times are given in Table 3. The similarity in the amplitude decay in the wavefront of the 3-D and 2.5-D simulations reveals the spatial and temporal damped nature of the 2.5-D wave equation.

Table 3.

Maximum positive and negative amplitudes in the -D and 2.5-D snapshots of the radial component.

Conclusions

An algorithm was developed using Takenaka & Kennett's (1996a) 2.5-D elastodynamic wave equation. The main aim of this development was to present a new algorithm to predict the site-specific ground motion around the epicentre, taking into account the 3-D radiation pattern and the equivalent moment tensor source formulation in a 2-D environment. The similarity of the 2.5-D and 3-D amplitude behaviour and the generation of different seismic phases at discontinuities in both cases proves the accuracy and the validity of the 2.5-D wave equation. Furthermore, analysis of the nature of the amplitude decay with time in the 3-D and 2.5-D wavefronts reveals the spatial and temporal damped nature of the 2.5-D wave equation; however, it is somewhat smaller than in the 3-D case. The slowness vector (p_y) was kept in the supersonic region so that all the components of the inertia coefficient are positive. The requirements of the spatial and temporal sampling in 2.5-D simulations to avoid grid dispersion and instability are the same as for 2-D simulations.

It was concluded on the basis of iterative 3-D and 2.5-D numerical experiments that moment per unit volume should be used in the 2.5-D simulation instead of moment per unit area. Moreover, it is necessary to scale down the required moment per unit volume by a factor equal to the finite-difference grid size in the 2.5-D modelling. The similarity of the computed 3-D and 2.5-D radiation patterns for different focal mechanisms in the xz-plane reveals that the 3-D radiation pattern can be generated in the 2-D environment using the 2.5-D wave equation. It is concluded that site-specific ground motion can be predicted using 2.5-D modelling, taking into account the 3-D radiation pattern, by rotating the strike of the fault counter-clockwise by the same amount as the azimuth of the site with respect to the epicentre. This main conclusion was reached on the basis of a detailed analysis of the 3-D surface radiation patterns and the corresponding 3-D and 2.5-D responses for various source mechanisms. The development of 2.5-D modelling for site-specific ground motion prediction taking into account the 3-D radiation pattern and local site conditions will be most useful for the purpose of earthquake engineering in areas without ground motion recordings.

Acknowledgments

Financial assistance from the Indian National Science Academy (INSA), New Delhi is gratefully acknowledged.