SIMULATION OF HIGH-FREQUENCY WAVE
PROPAGATION IN COMPLEX CRUSTAL
WAVEGUIDES USING GENERALIZED
SCREEN PROPAGATORS
RU-SHANWU, XIAN-YUNWU ANDXIAO-BIXIE
Modeling and Imaging Laboratory, Institute of Geophysics and Planetary Physics, University of California, Santa Cruz, California, USA
ABSTRACT
In the crustal waveguide environment, the major part of wave energy is carried by forward propagating waves, including forward scattered waves. Therefore, neglect-ing backscattered waves in numerical modelneglect-ing will not modify the main features of regional waves in most cases. By neglecting backscattering in the theory, the wave modeling becomes a forward marching problem in which the next step of propaga-tion depends only on the present values of the wavefield in a transverse cross-secpropaga-tion and the heterogeneities between the present cross-section and the next one (wave-field extrapolation interval). The saving of computation time and computer memory is enormous. A half-space screen propagator (generalized screen propagator) has been developed to accommodate the free-surface boundary condition for modeling SH wave propagation in complex crustal waveguides. The SH screen propagator has also been extended to handle irregular surface topography using conformal or non-conformal topographic transforms. The screen propagator for modeling regional SH waves has been calibrated extensively against some full-wave methods, such as the wavenumber integration, finite-difference and boundary element methods, for differ-ent crustal models. Excelldiffer-ent agreemdiffer-ent with these full-wave methods demonstrated the validity and accuracy of the new one-way propagator method. For medium size problems, the screen-propagator method is 2–3 orders of magnitude faster than finite-difference methods. It has been used for the simulation of Lg propagation in crustal models with random heterogeneities and the related energy partition, atten-uation and blockage. It is found that the leakage attenatten-uation of Lg waves caused by large-angle forward scattering by random heterogeneities, which scatters the guided waves out of the trapped modes and leaking into the mantle, may contribute signif-icantly to Lg attenuation and blockage in some regions. In the case of P-SV elastic screen propagators, plane wave reflection calculations have been incorporated into the elastic screen method to handle the free surface. Body waves, including the reflected and converted waves, can be calculated by real wavenumber integration; while surface waves (Rayleigh waves) can be obtained with imaginary wavenumber integration. Numerical tests proved the validity of the theory and methods.
Keywords: Lg-wave, Crustal wave guide, One-way propagator, Seismic wave
scat-tering
© 2006 Elsevier Inc. All rights reserved ISSN: 0065-2687 DOI 10.1016/S0065-2687(06)48006-7 323
1. INTRODUCTION
High-frequency regional wave propagation in complex crustal waveguides is one of the most challenging problems in theoretical and computational seismol-ogy. A good understanding of propagation, scattering, attenuation and wave-type conversion of regional waves and the availability of analytical/numerical tools to simulate and analyze these phenomena for complex crustal structures, includ-ing rough surface, Moho topography and small-scale heterogeneities, are crucial to the applications of regional waves to various geophysical problems. Regional wave tomography for crustal structures, path correction for discrimination and yield estimation of low-yield nuclear tests, location determination of earthquakes, or underground explosions using regional phases are examples among the possi-ble applications. Nuclear explosion monitoring at regional distances is even more demanding for the simulation and analyzing tools. For this purpose, simulation algorithms are desirable for generating synthetic waveforms for high frequencies up to 25 Hz at distances greater than 1000 km.
Substantial efforts have been made in modeling regional wave propagation. Methods based on layered earth models, such as the reflectivity and mode sum-mation methods (e.g.,Bouchon et al., 1985; Kennett, 1989, 1990; Maupin, 1989; Baumgardt, 1990; Campillo, 1990; Campillo and Paul, 1992; Campillo et al., 1993; Gibson and Campillo, 1994) have very high efficiency and can be applied to relatively high frequencies, but they can be used only for very simplified cases with layered or smoothly varying layered models, or be applied to part of the wavefield. Modeling techniques that can treat realistic 3-D heterogeneous me-dia, rather than smoothly varying layered meme-dia, are needed to test and study many observations and hypotheses. New modeling methods are needed to handle sudden changes of crustal thickness, strong lateral variations and irregular 3-D heterogeneities. As pointed out byCampillo et al. (1993), actual Lg amplitudes are reduced more than 10 times for paths passing through an anomalous zone on the east side of the Alpine range, while the modeling results using existing methods (including the effect of known large-scale lateral structural variation) only account for 20–30% of the amplitude reduction. Other attenuation mecha-nisms such as the scattering and attenuation by small-scale heterogeneities must be taken into account.
Kennett (1984, 1998) and Maupin and Kennett (1987) developed a cou-pled mode method for calculating guided seismic waves in horizontally varying structures. The method works well for relatively low frequency waves in mod-erately heterogeneous models (for a summary of the coupling mode method, see Chapter 2 of this book by Maupin). However, the implementation of the method for high frequency 3-D models still requires formidable computational efforts.
Chen (1990, 1995)developed a global R/T (Reflection/Transmission) matrices method, for simulating the seismic wave excitation and propagation in an
arbi-trarily multi-layered medium with irregular interfaces, which can be regarded as an extension of the generalized R/T coefficients method (reflectivity method) for the horizontally layered case by incorporating the T-matrix approach (for a sum-mary, seeChapter 4of this book by Chen). Again the application of the method is limited to low frequencies and short propagation distances.
Cormier and Anderson (1996, 1997) applied elastic Born scattering (in the regime of Rayleigh scattering) to the locked-mode solution for plane layered media to calculate the effects of small-scale heterogeneities. However, the ap-proximation is limited to single scattering and is only good for heterogeneities with scales much smaller than the wavelength. Ray method has very limited suc-cess in modeling regional waves due to the chaotic behavior of rays caused by the multiple reflections from the free-surface and Moho. Keers et al. (1996a, 1996b) applied the Maslov integral method to avoid the caustics and pseudo-caustics (pseudo-caustics of plane waves) by working in the phase-space. However, when chaos develops in the ray system, more complicated caustics arise for which the Maslov method does not work. In addition, ray-tracing computation is very time-consuming in this case. An alternative and flexible approach using ray-tracing has been developed byKennett (1986),Bostock and Kennett (1990)andKennett
et al. (1990), in which ray diagrams are used to study Lg waves crossing struc-tural boundaries. The method agrees well with modal calculations and can be applied to surface topography, 3-D crustal structures and other cases. However, the method cannot provide information on wave phenomena for complicated waveguides.
Finite-difference methods (e.g.,Xie and Lay, 1994; App et al., 1996; Goldstein
et al., 1996, 1997, 1999; Husebye and Ruud, 1996; Jih, 1996; Nolte et al., 1996;
Jones et al., 1997; McLaughlin and Wilkins, 1997; Bradley and Jones, 1998, 1999; Xie et al., 2005) and pseudo-spectral methods (e.g.,Kosloff et al., 1990; Archambeau et al., 1996; Schatzman, 1996; Furumura and Kennett, 1997) are commonly-used numerical methods that have been extensively applied to re-gional wave propagation. Theoretically, these methods can deal with arbitrarily heterogeneous media. However, it is necessary to use very dense spatial sam-pling to avoid grid dispersion for long distance regional wave propagation (for grid dispersion problem, see Fornberg, 1987). The capability of the present-day computers usually restricts them to short propagation ranges and relatively low frequencies, which prevents them from being applied to more realistic cases.
The state-of-the-art of the traditional simulation techniques for regional waves has its application to relatively low frequencies and short propagation distances. Correspondingly, the volume heterogeneities and surface irregularity in the crustal models are limited to rather large scales. However, high-frequency regional waves up to 20 Hz or higher have been observed over different distances, ranging from a few hundred kilometers to more than one thousand kilometers (e.g., Ni et al., 1996; Herrmann et al., 1997; Lay et al., 1999). Since high-frequency waves
can be used for event locations with high accuracy, simulation and modeling of high-frequency regional wave propagation are very desirable for many applica-tions. For high-frequency wave propagation, scattering and attenuation, the role of small-scale heterogeneities and surface roughness are all important.
The existence of small-scale heterogeneities in the crust and the associated seismic wave scattering have been known among seismologists (e.g., Aki and Richards, 1980; Wu and Aki, 1988, 1989, 1990; Sato and Fehler, 1998). How-ever, the effects of these heterogeneities on guided wave (Lg) propagation in the crust have not been explored extensively. The reasons may be the following. First, the spectra, strength and distribution of the small-scale heterogeneities in differ-ent regions are not well-known. Very few data sets can be used to characterize the paths concerned. Second, there lack analytical and numerical tools to model or analyze their influence on the guided wave propagation. The theory of wave propagation in unbounded random media has been well developed. However, for waves in complex crustal waveguides with random heterogeneities, the theoretical difficulties are overwhelming, and no analytical tools are available for performing realistic calculations. Therefore, numerical methods for simulating regional wave propagation in complex waveguides with small-scale heterogeneities are highly desirable. It has become clear that small-scale heterogeneities are widely distrib-uted in tectonically active regions. Strong topographic variation is the manifesta-tion of tectonically active regions and often the indicamanifesta-tion of small-scale hetero-geneities.Figure 1gives a topographic profile (top panel) and its power spectrum (bottom panel) for a path crossing the Tibet region. The slope of the spectral roll-off is close to 1/k, a flicker noise spectrum, very rich in small-scale variations. This spectrum is similar to the observations of the sonic well-log in the KTB super-deep continental drilling well (Wu et al., 1994; Jones and Holliger, 1997; Goff and Holliger, 1999), where the spectrum also has a 1/k slope. Recently,Goff and Holliger (1999)explained the 1/k spectra as a combination of hierarchical, multi-scale heterogeneities. Overall, the 1/k spectra demonstrate the richness of small-scale heterogeneities.
Recently, the generalized screen method has been introduced into seismic wave simulations and applied to the problems of both exploration and theoretical seis-mology. The generalized screen method is based on the one-way wave equation and the one-return approximation. The one-way generalized screen propaga-tor (GSP) neglects backscattered waves, but correctly handles all the forward multiple-scattering effects, e.g., focusing/defocusing, diffraction, interference, and conversion between different wave types. The one-return approximation is also called the De Wolf approximation (De Wolf, 1971, 1985), which neglects the reverberation between screens and can simulate multiple-forescattering–single-backscattering (MFSB). Significant progress has been made on the development of an elastic complex screen (ECS) method for modeling elastic wave propaga-tion and scattering in arbitrarily complicated structures (Wu, 1994, 1996; Xie and Wu, 1995, 2001; Wild and Hudson, 1998; Wu and Wu, 1999). The method is two
FIG. 1. Topographic profile (top) and its power spectrum (bottom) for a path crossing the Tibet region.
to three orders of magnitude faster than the elastic finite-difference method for a medium sized 3-D problem. For detailed derivation and the physical meaning of the MFSB approximation and the generalized screen method, seeChapter 5 of this volume. The screen method has been successfully used in forward model-ing (Wu, 1994; Wu and Huang, 1995; Wu et al., 1995; Xie and Wu, 1995, 1996, 1999, 2001; Wu and Wu, 2001, 2005) and as backpropagators for seismic wave imaging/migration in both acoustic and elastic media (e.g.,Wu and Xie, 1994; Huang and Wu, 1996; Huang et al., 1999a, 1999b; Jin and Wu, 1999; Jin et al., 1999; Xie and Wu, 1998, 2005; Xie et al., 2000).
The rest of this chapter is divided into three sections. We first briefly describe the basic concept of the one-way propagator for simulating elastic wave propaga-tion in crustal wave guide. In the second part, we give a systematic review of the screen propagator for the SH wave propagation in complex crustal waveguides including irregular topography. In the last part, we present the P-SV Lg screen propagator.
2. A BRIEFDESCRIPTION OF THEGENERALIZEDSCREENPROPAGATOR FORGUIDEDWAVES
In the crustal wave guide environment, major part of wave energy is carried by forward propagating waves, including forward scattered waves. The Lg en-ergy, which is in the form of guided waves, is carried by forward propagating waves bouncing up and down between the free surface and major geophysical discontinuities such as the Moho and Conrad discontinuities. Beyond the critical reflection angle, these waves are systematically dominated by small-angle waves (relative to the main propagation direction) trapped in the crustal waveguides. Therefore, neglecting backscattered waves in the propagation during numerical modeling will not modify the main features of regional waves in most cases. With this approximation, the modeling method becomes a forward marching al-gorithm in which the next step of propagation depends only on the present values of the wavefield in a transverse cross-section and the heterogeneities between the present cross-section and the next one (within an extrapolation interval). To for-mulate the problem, we divide the crustal wave guide into a sequence of vertical slabs. The horizontal direction is chosen as the main propagation direction. The geometry of the model is shown inFig. 2a. Choosing one slab as the example, Fig. 2b shows the interaction between the incident waves and the slab. By intro-ducing the local Born approximation, both wavefields and the elastic parameters can be separated into two parts, the background values and the perturbations. The “thin-slab” must be thin enough to satisfy the local Born approximation: the scat-tered field due to the heterogeneities in the slab be much smaller than the incident field. The incident P- and S-waves uP0 and uS0 enter the slab from the vertical plane at x0. After the incident waves pass through the thin-slab between x0and x1, and interacting with the heterogeneities within it, there will be both incident
(a) (b)
FIG. 2. (a) Geometry using screen method to simulate Lg wave; (b) Sketch showing the interaction between the incident waves and a thin slab.
waves and different types of forward scattered waves at the exit plane at x1. The new P-wave uP = uP0 + UPP+ USP is composed of incident P-wave and scat-tered P-waves UPPand USP, respectively from incident P- and S-waves, and the new S-wave uS = uS0+ UPS+ USS is composed of incident S-wave and scat-tered S-waves UPSand USS, respectively. The propagation and scattering in the thin-slab can be formulated using the perturbation theory and calculated by two separated steps: (1) The interactions between the incoming waves and the hetero-geneities are conducted in the spatial domain, accounting for the scattering and the coupling between different wave types. (2) Plane wave propagation through the background medium is conducted in the wavenumber domain by simple phase-shift. In both domains, the calculations are local and highly efficient. There is no time-consuming spatial or wavenumber domain convolution involved. Forward and inverse fast Fourier transforms (FFT) shuttle the wavefield between the two domains. By iteratively using this process and making the output from one slab as the input of the next slab, the wavefield can be propagated through the entire model.
The conventional wavenumber integral method is for the horizontally layered model and the integral is along horizontal wavenumber kx. By contrast, the elastic
screen method for propagating guided waves in crustal environment uses vertical screens and the wavenumber integration is carried out along the vertical kzaxis.
Under this geometry, the postcritical reflections become small-angle events with respect to the main propagation direction (the x-direction), therefore, the gen-eralized screen propagator (GSP) methods based on small angle approximation is suitable for modeling Lg waves. Our discretized model is composed of vertical thin-slabs and therefore the wavenumber integral is along vertical wavenumber kz,
resulting in different features compared with the traditional wavenumber integra-tion.Figure 3is a sketch showing the difference between two integral axes. For the reflected P- and S-waves coupled at a horizontal free surface, both waves have the same horizontal wavenumber (Fig. 3a). Choosing the horizontal wavenum-ber as the integral variable, P- and S-waves have the same sampling points in the kx axis. The Rayleigh pole, which contributes to the generation of Rayleigh
wave, has a unique location on the kx axis (Figs. 3a and b). However,
choos-ing the vertical wavenumber kzas the integral variable, the P- and S-waves have
different sampling points. The Rayleigh pole in the real kx axis splits into two
points in the imaginary kzaxis (Fig. 3c). This makes resampling necessary when
calculating converted waves. Fast Fourier transform with regular sampling inter-val cannot be used in the case of vertical wavenumber integration for converted waves. All these introduce additional complexity in wavenumber integration and special treatment is required. However, the general principle of plane wave (in-cluding the inhomogeneous plane waves) superposition for representing point sources still holds.
GSP is accurate for small-angle propagation and scattering (near horizon-tal for the crushorizon-tal waveguide environment). A half-space screen propagator has
FIG. 3. Sketch showing different wavenumber domain integrals: (a) Dispersion relation in the kx−kz plane; (b) Traditional wavenumber integration in the complex kx plane; (c) Wavenumber integration in the complex kzplane for the generalized screen method using vertical screens in a half-space.
been introduced by Wu et al. (1996, 1997, 1998, 2000a, 2000b) to accom-modate the free-surface boundary condition and treat the SH wave propaga-tion in complex crustal waveguides. The new one-way method for modeling regional SH waves has been calibrated extensively with various full-wave meth-ods for different crustal models, such as the wavenumber integration method for flat structures and full-wave finite-difference method for heterogeneous crustal waveguides. Excellent agreement with these methods demonstrated the validity and accuracy of the new one-way method. For a model with propagation dis-tance of 250 km, dominant frequency at 0.5 Hz and with similar accuracy, the GSP method is about 300 times faster than the finite-difference method. The GSP method has been applied to the simulation of Lg propagation in random media for calculating the related energy partition and attenuation (Wu et al., 2000b). It is found that the leakage attenuation of Lg waves caused by forward large-angle scattering from random heterogeneities, which scatters the guided waves out of the trapped modes resulting in energy leaking into the mantle, may contribute significantly to Lg attenuation and blockage in some regions. The apparent Q for leakage attenuation, as a function of normalized scale length
ka of the random heterogeneities, agrees well with the scattering theory. Later,
the SH screen propagator is extended to the case of irregular surface topogra-phy by conformal or non-conformal topographic transforms (Wu et al., 1999; Wu and Wu, 2001). In the conformal transform method, the coordinate system is rotated according to the local topographic slope, and the mirror image method can be applied to the local plane surface; the non-conformal method is a surface flattening transform which turns the free surface topography into modified vol-ume perturbations of elastic parameters. The former method is suitable to deal with smoothly varying topography, while the latter method can treat rough but moderate topography.
In the P-SV wave case, the derivation and application of one-way GSP screen propagators are much more complicated. Unlike for SH waves, the mirror image method of generating the half-space Green’s function cannot be used to account for the effect of the free surface. Plane wave reflection calculations are incorpo-rated into the elastic screen method (Wu et al., 2000c). Body waves, including the reflected and converted waves, can be calculated by real wavenumber inte-gration; while surface waves (Rayleigh waves) can be obtained with imaginary wavenumber integration. Numerical tests show good agreement with the theory. More work has to be conducted for the coupling among the body waves, guided wave and surface wave caused by lateral heterogeneities and irregular topography.
3. SH WAVECASE 3.1. Half-Space Screen Propagator
For a 2-D SH problem, only the y-component of the displacement field, noted as u, exists. With the perturbation method, the medium and the wave field are decomposed into
ρ = ρ0+ δρ, µ= µ0+ δµ, u= u0+ U,
where ρ0and µ0are the density and shear rigidity of the background medium, δρ and δµ are their corresponding perturbations, u0is the primary field and U is the scattered field. The SH wave equation in the frequency domain can be written as
(1)
µ0∇2u+ ω2ρ0u= −
ω2δρu+ ∇ · δµ∇u,
where∇ is the 2-D gradient operator and “·” stands for inner product. Equation(1) is a scalar Helmholtz equation. With a half-space scalar Green’s function gh, the
scattered field U can be written as (2) U (r1)= k2 v d2r gh(r1; r)ερ(r)u(r)− 1 k2∇g h(r 1; r) · εµ(r)∇u(r) .
Under the forward-scattering approximation, or more generally the multiple-forescattering–single-backscattering (MFSB) (De Wolf, 1971, 1985; Wu and Huang, 1995; Wu, 1996), the total field and Green’s function under the integration in the above equation can be replaced by their forward-scattering approximated counterparts, and the field can be calculated by a one-way marching algorithm along the x-direction using a dual domain technique (seeChapter 5of this book). Note that the half-space Green’s function must be used here to account for the free surface effect.
For each step of the marching algorithm under the forward-scattering approx-imation, the total field at x1 is calculated as the sum of the primary field which is the field propagating in the half-space from x to x1, and the scattered field caused by the heterogeneities in the thin-slab between xand x1. The thickness of the slab should be made thin enough to ensure the validity of the local Born approximation. The Green’s function in the homogeneous half-space can be ob-tained by the mirror image method. The stress should vanish at the free surface
z= 0. Therefore we have (3) g0h(r1; r) = g0(r1; r) + g0 r1; r∗ ,
where g0is the infinite homogeneous Green’s function and r∗is the mirror image point of r with respect to the free surface.
The free space Green’s function in wavenumber domain is (Wu, 1996) (4) g0(x1, Kz; x, z) = 1 2γe iγ|x1−x|e−iKzz with (5) γ = k2− K2 z. Therefore, (6) g0h(x1, Kz; x, z) = 1 2γe iγ (x1−x)2 cos(K zz).
In a similar way we can obtain
(7) ∇gh 0 = 1 2e iγ (x1−x) ˆe
x2 cos(Kzz)− ˆez(Kz/γ )2i sin(Kzz)
,
whereˆexandˆezare the unit vectors in the x- and z-directions, respectively.
Taking Fourier transform of Eq.(2)along z1for the case of a thin-slab perpen-dicular to the x-direction, and substitute the half-space Green’s functions into it, the scattered field by the thin-slab can be calculated by (for details seeWu et al.,
2000a) U (x1, Kz)= Uρ(x1, Kz)+ Uµ(x1, Kz), Uρ(x1, Kz)= ik x1 x dx eiγ (x1−x)C k γερ(z)u0(z) , (8) Uµ(x1, Kz)= ik x1 x dx eiγ (x1−x) × C εµ(z) ¯∂xu0(z) − iS Kz γ εµ(z) ¯∂zu0(z) , where ερ(r)= δρ(r) ρ0 , εµ(r)= δµ(r) µ0 , and γ = k2− K2
z is the propagating wavenumber in the x-direction, Kz is
transverse wavenumber along the z-axis, and ¯∂x=
1
ik ∂ ∂x
are dimensionless partial derivatives. In the above equations, C[f (z)] and S[f (z)] are the cosine and sine transforms:
C f (z)= ∞ 0 dz 2 cos(Kzz)f (z), (9) S f (z)= ∞ 0 dz 2 sin(Kzz)f (z),
and u0, ¯∂xu0and ¯∂zu0can be calculated by u0(x, z)= 1 2π ∞ −∞dK zeiK zzeiγ(x−x)u 0(x,Kz) (10) = C−1 eiγ(x−x)u0(x, Kz) , ¯∂xu0(x, z)= C−1 eiγ(x−x)γ ku0(x , K z) , (11) ¯∂zu0(x, z)= iS−1 eiγ(x−x)K z k u0(x , K z) .
The above equations are the general wide-angle formulation. When the energy of crustal guided waves is carried mainly by small-angle waves (with respect to the horizontal direction), the phase-screen approximation can be invoked to sim-plify the theory and calculations. Summing up the primary and scattered fields
and invoking the Rytov transform results in the dual-domain expression of phase-screen propagator (12) u(x1, Kz)≈ eiγ (x1−x ) C eikSs(z)u 0(x,z) ,
where eikSs(z)is the phase delay operator with
(13) Ss(z)= 1 2 x1 x dx ερ(x, z)− εµ(x, z) ≈ x ¯εs(z),
where ¯εs(z)is the average S-wave slowness perturbation over the thin-slab at
depth z, (14) ¯εs(z)= 1 x1− x x1 x dxs(x, z)− s0 s0 ,
with s(x, z) = 1/v(x, z) and x = (x1− x)is the thin-slab thickness. Equa-tion(12)is the SH phase-screen propagator for the half space. It has a similar form as the whole space propagator with the Fourier transform replaced by a co-sine transform.
The phase-screen propagator has long been used in ocean acoustics to simulate long range acoustic wave propagation in the heterogeneous ocean due to internal waves. Most work in the literature deals with the stochastic treatment of waves in random media. For an introduction and brief summary of the work in that field, the reader is referred toFlatté et al. (1979). However, in this work we will use the half-space screen propagator for deterministic modeling in heterogeneous crustal waveguides.
3.2. Treatment of the Moho Discontinuity
The Moho discontinuity can be treated in two ways. One way is to put the impedance boundary conditions in the formulation, the other is to treat the pa-rameter changes as perturbations and therefore incorporate the discontinuity into the screen interaction. In this paper, we adopt the latter approach because of its flexibility in treating irregular Moho discontinuity. The validity of the perturbation approach for the Moho discontinuity is verified by the comparison with wavenum-ber integration and full-wave finite-difference algorithms. Since for guided waves, or crustal waves with critical or post-critical reflections, the related mantle waves are nearly horizontal, the screen approximation is quite accurate in this case. The excellent agreements of the method with the wavenumber integration for flat Moho, and with the finite-difference method for irregular Moho demonstrate the validity of this approach (Wu et al., 2000a).
Figure 4compares the reflection coefficients of the Moho discontinuity calcu-lated using the theoretical equation (dotted line) and using phase-screen method (solid line). A constant velocity crust model (vc = 3.5 km/s, ρc = 2.8 g/cm3,
FIG. 4. Comparison of reflection coefficients at the Moho discontinuity. Dotted line denotes result from theoretical equation and solid line denotes result calculated using screen method. A constant velocity crust model is used in the calculation and the source is located 30 km above the Moho dis-continuity.
vm = 4.5 km/s, ρm = 3.1 g/cm3) is used in the calculation and the source is
located 30 km above the Moho. The curve from the screen method is obtained by using the root-mean-square (RMS) of the reflected waveforms. We see that the process of critical reflection is well matched, except that the transition from pre-critical to pre-critical calculated by the screen method is not as sharp as the theoretical curve. This may be caused by the small phase error in the phase-screen approxi-mation. In addition, the reflection for wide-angle incidence, especially for nearly vertical incidence, is not well modeled by the phase-screen method. This error results from the small angle approximation used in the screen formulation. How-ever, this limitation occurs only at short range, well before the critical distance (around 80 km in this case). As can be verified using numerical simulations, the screen method behaves exceptionally well beyond the critical distance, making it a good candidate for guided wave simulation.
3.3. Numerical Verifications and Simulation Examples
In this section we show some examples demonstrating the validity of the method and its potential applications to various problems of regional wave prop-agation. First, we show a comparison between the screen method and a full-wave finite-difference method for a heterogeneous crustal model. Shown inFig. 5a is a wave guide model with a crust necking.Figure 5b shows the synthetic seis-mograms along a vertical profile at an epicentral distance of 250 km. The thin lines are from the finite-difference method and the thick lines are from the screen method. The source is located at a depth of 2 km and the source time function has a dominant frequency of 0.5 Hz.Figure 5demonstrates excellent agreement
(a)
(b)
FIG. 5. Comparison of synthetic seismograms along a vertical profile at a distance of 250 km. Shown in (a) is the velocity model with a laterally varying crustal wave guide. Shown in (b) are synthetic seismograms calculated using the screen method (thick lines) and a finite-difference method (thin lines). The source depth is 2 km and the source time function is a Gaussian derivative with a dominant frequency of 0.5 Hz.
between the two methods. For this example, the GSP method is about 300 times faster than the finite-difference method. Note that the grid spacing used in the FD simulation was 3–4 times smaller than the stability requirement in order to reduce the numerical dispersion. Other comparisons with wavenumber integration and finite-difference methods can be found inWu et al. (2000a, 2000b).
The importance of small-scale random heterogeneities to seismic wave prop-agation is well known. There are extensive publications on this subject in seis-mology. However, due to the complexity of the problem, the role of random heterogeneities in Lg excitation, propagation, attenuation and blockage is still unclear. For waves in complex crustal waveguides with random heterogeneities, there are still no analytical tools available for performing realistic calculations. Numerical simulation is a useful alternative to the theory. Some finite-difference simulations have been conducted (e.g.,Frankel and Clayton, 1986; Frankel, 1989; Xie and Lay, 1994; Jih, 1996). Due to the limit of the computing power, the wave-propagation distances in the finite-difference simulations are relatively short.
FIG. 6. A heterogeneous crustal model representing a mountain root with small-scale random het-erogeneities (top panel). The comparisons between synthetic seismograms with and without random heterogeneities are shown on the middle and bottom panels, respectively.
Liu and Wu (1994)have done some numerical simulations using the phase-screen method, but the models simulated are limited to unbounded media. The devel-opment of the half-space GSP method enables us to simulate long distances, high-frequency wave propagation in complex crustal waveguides. We present here two numerical examples to demonstrate the capability of the GSP method.
Figure 6shows a heterogeneous crustal model representing a “mountain root” with small-scale random heterogeneities. The top panel is the velocity model, and the comparisons between synthetic seismograms with and without random heterogeneities are shown on the middle and bottom panels, respectively. The het-erogeneities have an exponential correlation function, with the scale length ax = az = 1.6 km (in horizontal and vertical directions, respectively). The RMS
ve-locity perturbation is 5%. The dominant frequency of the source time function is 2 Hz.Figures 7a and b show a comparison of wavefield snapshots between mod-els with and without random heterogeneities. We see that random heterogeneities drastically increase the complexity of the wavefield and the energy leakage to the upper mantle.
FIG. 7. Comparison between snapshots for waves passing through a “mountain root” with or with-out random heterogeneities, shown on A and B, respectively.
3.4. Application to Energy Partition and Attenuation in Crustal Waveguide with Random Heterogeneities
In heterogeneous crustal wave guides, the upper boundary is the free surface, which is a perfect reflector. The lower boundary of the wave guide is the Moho discontinuity. For waves incident on the Moho discontinuity, a part of the en-ergy will leak into the upper mantle. However, for waves incident on the Moho with post-critical angles, total reflections occur and all the energy is reflected and trapped in the waveguide. Generally speaking, the guided wave energy can be expressed as (15) Eg= C Kz<Kc u(Kz) 2 dKz,
where C is a constant, Kzis the wavenumber in the z-direction, namely the
trans-verse wavenumber, and Kcis the critical wavenumber. Scattering processes can
redistribute the energy in wavenumber domain, causing the leak of trapped energy into the upper mantle. In addition to the leakage loss, the guided waves suffer also the anelastic loss and backscattering loss. Assuming a homogeneous mantle and neglecting reverberation in the x-direction, the energy balance after propagating a short distance dx in the x-direction is
(16)
Eg(x+ dx) = Eg(x)− Ea(x)− Eb(x)− El(x),
where Egis the energy of guided crustal waves; Ea, energy lost due to absorption
(anelastic loss); Eb, energy lost due to backscattering by random heterogeneities; El, energy lost due to leakage to the mantle caused by heterogeneities. In terms
of different attenuation coefficients, it can be written as
(17) dEg
dx = −[ηa+ ηb+ ηl]Eg(x)= −ηgEg(x),
where ηa = (Ea/Eg)/dx, ηb = (Eb/Eg)/dx, and ηl = (El/Eg)/dx are the
apparent attenuation coefficients for guided crustal waves. Equivalently,
(18)
ηg= kQ−1g = k
Q−1a + Q−1b + Q−1l ,
where Q’s are the corresponding apparent quality factors.
The leakage loss is the scattering loss due to the redistribution of Lg angular spectra. It is caused dominantly by large-angle forward scattering and therefore, it is several orders of magnitude larger than the backscattering loss, i.e. ηl ηb.
In the following, we will concentrate on the analysis of leakage loss of guided waves. For the leakage analysis, the angular spectrum representation or the energy distribution versus propagation angle (or vertical slowness) will be very useful and can show clearly which part of the energy would be trapped in the wave guide and which part of the energy would leak into the mantle.
In first-order approximation, the anelastic (intrinsic) attenuation is additive to the leakage loss, so that we can calculate and analyze these two attenuation mech-anisms separately. For the Lg RMS amplitude attenuation, one more attenuation mechanism is involved:
(19)
bg= kQ−1g = k
Q−1a + Q−1b + Q−1l + Q−1d
where Q−1d is the equivalent Q of diffusion loss, which represents the amplitude decrease of Lg due to the transfer of coherent energy into incoherent energy (Lg coda) by random heterogeneities.
Shown in Fig. 8is a comparison between angular spectra from a waveguide model with 5% RMS random velocity perturbation in the crust and a reference flat crust without velocity perturbation. The perturbation has an exponential cor-relation function with horizontal and vertical characteristic scales (corcor-relation lengths) of 5.0 and 3.0 km, respectively. From the top panel to the bottom panel
FIG. 8. Energy distribution for different crustal models. From top to bottom are: waveguide model with 5% RMS velocity perturbations in the crust; energy angular spectra versus distance for a flat crust; energy angular spectra versus distance for random crust; and relative energy attenuations versus distance. The dotted line is for the flat crust model and the solid line is for the random crust model.
inFig. 8, they are a random velocity model, energy distribution for homogeneous crust, energy distribution for random crust; and energy attenuation curves versus distances, respectively. For energy distributions, the vertical coordinate is the
nor-malized vertical slowness Kz/k, corresponding to the cosine of incident angles
(or sine of the grazing angles). Note that zero vertical slowness means horizontal propagation. The frequency range is between 0.6 and 1.9 Hz. For the flat crust model, there is a considerable portion of energy with large vertical slowness (or steep angles) at the initial stage. After multiple reflections, energy with larger ver-tical slowness is depleted due to the leakage to the mantle, leaving the energy with small vertical slowness, i.e., the guided waves, propagating in the waveguide. For the model with random velocity perturbations, the distinct feature is the contin-uous energy repartition, moving from small (grazing) angle waves to large-angle waves due to scattering by small scale heterogeneities. The energy propagating with large angles tends to leak into the mantle and causes Lg-wave energy at-tenuation. The bottom panel ofFig. 8is the wave energy attenuation versus the distance. The energy is calculated from synthetic seismograms on the free sur-face. The dotted line is for the reference (homogeneous) crust model. It can be seen that for this case, after passing 100 km or more, the energy is basically kept constant, which means that the trapped mode has been formed. The solid line is for the random waveguide. Due to the scattering, the energy is decreasing with distance.
Figure 9 gives the attenuation curves for different characteristic scales. The upper panel is the attenuation curve of total energy, which is the energy con-tained in the entire seismogram recorded on the surface. The thin solid line is for ka = 1, the thick solid line is for ka = 10, and the dashed line is for the reference (homogeneous) model. We see that for the reference model, the total energy remains constant beyond critical distance, which serves as a check-ing point for the numerical simulations. The middle panel gives the coherent Lg energy, which is calculated using waves within the Lg window (group ve-locity between 3.7 km/s and 3.2 km/s) versus distance. Again, the thin, thick and dashed lines are for ka = 1, ka = 10 and the reference model, respec-tively. In both measurements, the cases with ka= 1 are always associated with stronger attenuation than ka = 10 cases. We also see that the coherent Lg en-ergy corresponding to the peak amplitude suffers more attenuations than the total energy. This is due to the extra attenuation, i.e. the diffusion loss which scatters the waves out of the Lg window and transfers them into incoherent waves (Lg coda). However, in these numerical simulations, there is no intrinsic attenuation, and leakage attenuation dominates. The difference between the coherent energy attenuation and the total energy attenuation is relatively small. In the bottom panel ofFig. 9, we plot the curve of apparent inverse quality factor for leakage atten-uation Q−1l versus the normalized scale length (ka) of random heterogeneities, where k = 2π/λ with λ being the wavelength of the dominant frequency, and
a the correlation length. Since no intrinsic (anelastic) attenuation exists in the model and no backscattering is involved, the attenuation is solely caused by the leakage loss due to scattering. From the curve we see that Q−1l reaches its peak at ka ≈ 1.5–2.0 and keeps flat until ka ≈ 8.0. This is a feature of large-angle
FIG. 9. Total energy attenuation (top panel), and windowed Lg energy attenuation between group velocities 3.1 km/s and 3.7 km/s (middle panel) versus distance for ka= 1 ( thin lines) and ka = 10 (thick lines). The bottom panel shows the equivalent Q−1for leakage attenuation versus the normal-ized scale length ka. The dashed line is for the reference model of a homogeneous crust.
forward-scattering dominance. For backscattering, the maximum scattering Q−1 is around ka ≈ 1.0 and decreases rapidly at ka > 1.0 for exponential correla-tion funccorrela-tions; while for large-angle forescattering, the plateau is quite wide after
ka = 1.0 (Wu, 1982; Frankel and Clayton, 1986). The numerical simulations agree well with the scattering theory. The values of the equivalent Q (300–900 for f0= 1 Hz) are comparable with some observations (Xie and Mitchell, 1991;
Xie, 1993). This suggests that the leakage attenuation caused by small-scale ran-dom heterogeneities may be responsible and even the ran-dominant mechanism for some observed Lg attenuations and blockages.
3.5. SH-Waves in Crustal Waveguides with Irregular Surface Topography Theoretical studies and observations show that surface topography is one of the important factors affecting Lg wave propagation. For example, irregular surface can cause anomalous variation of Lg amplitude along the propagation path (Sills, 1978; Geli et al., 1988; Bouchon and Barker, 1996). Methodologi-cally, range-independent boundary conditions for flat surface must be replaced by range-dependent boundary condition for an irregular surface. In the case of sur-face topography, the global mirror symmetry no longer exists. To use the GSP method for solving range-dependent boundary condition problems, both confor-mal and non-conforconfor-mal coordinate transforms were incorporated into the GSP method and their relative merits and accuracies were analyzed (Wu et al., 1999; Wu and Wu, 2001). The following is a summary of these two approaches.
3.5.1. Conformal Coordinate Transform Method for Smoothly Varying Topography
For a flat free surface, Wu et al. (2000a) derived a half-space GSP solution for Lg wave propagation. In the case of irregular topography, the global mirror symmetry for the problem no longer exists. However, taking a local plane-surface approximation for the topography, we can modify the mirror wavefield method to a local mirror wavefield method and apply the corresponding coordinate trans-form to obtain a GSP solution for the irregular topography.
Figure 10shows the geometry of the derivation. Assume that u+0(x, z)is the wavefield on the half-screen S+in the lower half-space. To calculate the wave-field in the next screen with the existence of a locally dipping surface, we first obtain the mirror wavefield u−0(˜x,˜z) on the half-screen S− in the upper half-space. The total wavefield in the next screen is composed of contributions from incident waves u+0(x, z)and u−0(˜x,˜z) plus the scattered field which is generated
by the local heterogeneities in the thin-slab. The effects of the heterogeneities and the topography can be calculated separately for each step in the GSP method. The effect of the slant free-surface can be incorporated into the propagation integral. Assume ut(x, z) is the total field including the scattering effect of the volume
inte-FIG. 10. Geometry of the conformal coordinate transform. gral u(x1, z1)= S ds g(x, z; x1, z1) ∂ut(x, z) ∂n − ∂g(x, z; x1, z1) ∂n ut(x, z) (20) = S− ds{· · ·} + S+ ds{· · ·},
where g(·) is the Green’s function for the full space with the background velocity,
S = S+ + S− is the integration surface composed of lower and upper half-surfaces S+ and S−, respectively. The Rayleigh integral can be used to replace the Kirchhoff integral for each half surface integral. For the lower half-surface the contribution of S+is u+t (x1, z1)= −2 ∞ 0 dz u+t (x, z) ∂g(x, z; x1, z1) ∂n (21) = 1 2π dKT eiKTz1u+t (x1, KT), where (22) u+t (x1, KT)= eiγ (x1−x) ∞ 0 dz1u+t (x, z1)e−iKTz1.
Here u+t (x, z)is the total wavefield composed of incident field u+0(x, z)and the scattered field U+(x, z)caused by the heterogeneities within the slab (seeWu, 1994; Wu et al., 2000a). If we put the slab entrance at x = xand the wavefield on the entrance surface S+as u+t (x, z), then
(23)
where U+(x, z)= k2 x1 x dx e−iγ (x1−x) ∞ 0 dz g(x, z; x1, z1)ερ(x, z)u0(x, z) (24) − 1 k2∇g(x, z; x1, z1)· εµ(x, z)∇u0(x, z) .
For the bent upper half surface, we perform a coordinate transform by clock-wise rotation of 2θ to a new coordinate system (˜x, ˜z). Taking the downward direction as positive z-direction and the rotation angle from x to z as positive, the relation connecting the two systems is
˜x = x cos 2θ + z sin 2θ,
(25) ˜z = −x sin 2θ + z cos 2θ.
In the new system, the surface S−is parallel to the˜z-axis, so that
(26) u−t (˜x1, ˜KT)= ei˜γ( ˜x1− ˜x) 0 −∞d˜z u− t (˜x,˜z)e−i ˜KT˜z ,
where u−t (˜x,−˜z)= u+t (x, z). The field in the space domain can be obtained
by synthesizing the contributions from all plane waves
(27) u−t (˜x1,˜z) = d ˜KTei˜γ( ˜x1− ˜x ) ei ˜KT˜zu− t (˜x, ˜KT), where (28) u−t (˜x, ˜KT)= 0 −∞d˜z u− t (˜x,˜z)e−i ˜KT˜z .
Transform back to the original coordinate system, resulting in
u−t (x1, z1)= d ˜KT exp i (˜γ cos 2θ − ˜KTsin 2θ )(x1− x) (29) + ( ˜γ sin 2θ + ˜KTcos 2θ )z1 u−t (˜x, ˜KT).
We see that in the original coordinate system, the effective transversal and propa-gating wavenumbers are
KT = ˜γ sin 2θ + ˜KT cos 2θ,
(30)
γ = ˜γ cos 2θ − ˜KTsin 2θ.
If we transform the ( ˜KT,˜γ) system into (KT, γ ),
(31)
u−t (x1, z1)=
dKT u−t (˜x, KTcos 2θ− γ sin 2θ)eiγ (x1−x
)
The total field is a summation of the contributions from both u+t (x1, z1) and u−t (x1, z1) u(x1, z1)= dKTeiγ (x1−x ) eiKTz1 (32) × u+t (x, KT)+ u−t (˜x, KT cos 2θ− γ sin 2θ) .
When small-angle waves prevail such as in the case of Lg propagation, the spectral interpolation in Eq.(32)can be avoided and replaced by operations in the space domain using a narrow-angle approximation. From(23), it can be seen that to calculate the reflection response we need to find the spectral component
u+t (− ˜KT). We try to obtain the approximate space-domain operations corre-sponding to the wavenumber-domain interpolation. We know that
(33) u+t (− ˜KT)= ∞ 0 dz ei(−KTcos 2θ+γ sin 2θ)zu+ t (z).
With narrow-angle approximation, γ ≈ k, therefore,
(34) u+t (− ˜KT)= ∞ 0 dzeiKTz 1 cos 2θe ik(tan 2θ )zu+ t z cos 2θ ,
where θ is the dipping angle of the free surface at x = x. We see that the wavenumber-domain interpolation is transformed into a space-domain operation which is a modulation plus a coordinate stretching. For a flat surface, Eq.(32) reduces to the original half-space GSP method (Wu et al., 2000a).
Shown inFig. 11are the synthetic seismograms obtained using the conformal screen method for a Gaussian hill model (Fig. 11a). The Gaussian hill is repre-sented by h(x)= −h0exp[−(x−x0)2/2σ2] with x0= 62.25 km, h0= 4 km, and σ = 9.129 km. Synthetic seismograms calculated with a more accurate
bound-ary element method (Fu and Wu, 2001) are also given as a reference. The solid lines are from the screen method and the dashed lines are from boundary element method. The comparison indicates that the screen method gives a satisfactory re-sult. It correctly modeled waveforms between distance 60 and 70 km, where two reflections from the convex free surface interfere with each other and generate complex waveforms. Note that the coordinate stretch z/cos2θ increases very fast with large angle θ , the conformal screen method works only for smoothly varying topography.
3.5.2. Non-Conformal Coordinate Transform Method for Rough Topography
Another alternative approach is to incorporate surface flattening transform into the GSP method. The transform converts surface height perturbations into mod-ified volume perturbations. In this way the range-dependent boundary condition becomes a stress release boundary condition on a flat surface in the new coor-dinate system where the half-space GSP method is applicable. The transform is
(a)
(b)
FIG. 11. (a) Velocity model with a Gaussian hill topography and (b) synthetic seismograms cal-culated from this model. For the calculation, dx = dz = 0.25 km and dt = 0.05 s. The source is located at a depth of 32 km and the dominant frequency of source time function is 3 Hz. Receivers are on the free surface. The solid lines are synthetic seismograms calculated using the screen method with a conformal transform, and the dashed lines are synthetic seismograms calculated with boundary element method (Fu and Wu, 2001).
defined as (Beillis and Tappert, 1979)
(35) χ = x,
ζ = z − h(x),
where h(x) is the height function of free surface. Equation(35)shows that the transform gives only a shift to depth variable z, i.e., depth measurement starts from the free surface. Thus, it is a non-conformal transform. Using the above transform and perturbation theory, the original half-space screen propagator be-comes (Wu and Wu, 2001)
ˆu(χ1, kς)= eiγ χC eik0Ss(ζ )C−1 ˆu 0(χ, kζ) (36) − Z(χ1)ˆµ(ζ) µ0 S−1 kζˆu0(χ, kζ) ,
where C and C−1are the forward and inverse cosine transforms, and S−1is the inverse sine transform defined by Eq.(9), µ0 is shear modulus of background medium, SS is the relative slowness perturbation of the thin-slab and given by
Eq.(13), and
(37)
Z(χ1)= h(x1)− h(x).
Equation(36)is expressed in the new coordinates (χ , ζ ). It is clear that the second term in the bracket in Eq.(36)comes from the roughness of topography, which is proportional to the height difference of the adjacent two screens for each for-ward step. For the upgoing slope, Z(χ ) < 0, the field scattered by topography is in-phase with the background field and strengthens the background field, while for downgoing slope Z(χ ) > 0, the field scattered by topography is out-phase with the background field and weakens the background field. Equation(36) is computationally efficient, in which all calculations can be done by FFT.
Shown in Fig. 12 are synthetic seismograms calculated using the non-conformal screen method for the Gaussian hill model shown inFig. 11. The solid lines are from the screen method and the dashed lines are from the bound-ary element method. The excellent agreement between the two methods is clearly seen except at the vicinity of the hill top where a small discrepancy exists both in wave shapes and amplitudes. The error can be reduced by using a smaller step
FIG. 12. Synthetic seismograms for a Gaussian hill model (Fig. 11a). The solid lines are calculated using the screen method with a non-conformal transform and the dashed lines are calculated using the boundary element method. The parameters for the calculation are the same as inFig. 11.
(a)
(b)
(c)
FIG. 13. (a) A crustal model with a rough random surface. The correlation length is 2.5 km, RMS perturbation is 0.6 km. (b) Synthetic seismograms, and (c) energy distribution versus horizontal distance. (b) and (c) show a comparison between the non-conformal screen method and BE method for a crustal waveguide with a rough random surface. The thick smoothly varying curve in (c) is calculated with finite difference method for a uniform waveguide. The source is located at the depth of 8 km, the dominant frequency of the source time function is 1 Hz.
length x. For forward marching algorithms, the step length x can be adjusted according to the roughness of topography. The more severe the topography is, the finer the step length x should be. Therefore, the non-conformal screen method can handle more severe topography than the conformal screen method.Figure 13a shows a crustal model with a rough random surface used for testing feasibil-ity and accuracy of the non-conformal screen method. The correlation length is 2.5 km, RMS height fluctuation is 0.6 km.Figures 13b and c show a comparison of synthetic seismograms calculated by the non-conformal screen method and BE method, and the corresponding energy attenuation curves, respectively. The thick smoothly varying curve inFig. 13c is the energy distribution for a similar waveguide but with a flat free surface. We see that the presence of a rough random surface makes the waveforms and attenuation curves more complicated. Except for large-angle Moho reflections, the results of the screen method agree well with those of the BE method. However, for this example, the screen method took about 35 minutes, while the BE method took about 72 hours.
Figure 14shows an investigation of the combining effect of rough topogra-phy and volume heterogeneity on Lg wave propagation using the non-conformal screen method. The rough topography is the same as shown inFig. 13. The het-erogeneities are velocity variations only. The correlation lengths are 6 km in range and 4 km in depth, RMS velocity fluctuations are 5% and 10%, respectively. The thickly dashed line calculated by finite-difference method for a uniform crustal waveguide is used as a reference. We see fromFig. 14that random heterogeneities combined with rough topography drastically increase the attenuation of high fre-quency Lg waves. This example shows that the non-conformal screen method can handle the effects of both volume heterogeneities and moderately rough topogra-phy on Lg wave propagation at long distances and high frequencies.
4. P-SV CASE
To introduce the P-SV elastic screen propagator for a flat free surface, the basic idea is to incorporate plane wave reflection calculation into elastic screen method (Wu, 1994; Wu et al., 2000c). The half space is extended upward in vertical di-rection from free surface. The extended part has the parameters of background medium and will be used to keep records of upgoing waves which can be used for the calculation of reflected/converted waves by the free surface. The incident P and S waves at vertical profile x= xcan be decomposed into a superposition of plane waves uP0(Kz, x)and uS0(Kz, x). The propagating waves are represented
by the real vertical wavenumbers, and the imaginary vertical wavenumbers corre-spond to the surface waves (inhomogeneous waves). Reflection at the free surface can be calculated at each step, and the total field including the reflected waves will interact with the heterogeneities. We will first treat the propagating waves (homogeneous waves) and then discuss the calculation of the fundamental mode Rayleigh wave as an example of the surface wave modeling.
TION OF HIGH-FREQ UENCY REGION A L W A V E 351
FIG. 14. Lg wave attenuation versus horizontal distances. A random medium whose correlation lengths are 6 km in range and 4 km in depth, and RMS velocity fluctuations are 5% and 10%, respectively. The source is located at a depth of 8 km, the dominant frequency (f0)of source time function is 2 Hz. In the figure, “rough” means the crust with rough topography, “ho” and “het” denote homogeneous and heterogeneous crustal models, respectively.
Applying the reflection coefficients, the free surface reflected P and S waves due to incident P wave can be expressed by
(38) uPP(x, z)= eiγα(x−x) dKzuP0(Kz, x)PPˆa1e−iKzz, (39) uPS(x, z)= eiγα(x−x) dKzuP0(Kz, x)PSˆa2e−iK ∗ zz, where γα = k2
α− Kz2is the propagating wavenumber for P waves (here in the x-direction) and Kz∗ =
k2β− k2
α+ Kz2 is the transverse wavenumber of
con-verted S waves determined by Snell’s law. Unit vectorsˆa1= (γα,−Kz)/kα and
ˆa2= (Kz∗, γα)/kβ·|uP0(Kz, x)| is the scalar spectrum of the incident P wave with
a transverse wavenumber Kz (here in z-direction). The reflected P and S waves
due to the incident plane S wave can be obtained by
(40) uSP(x, z)= eiγβ(x−x) dKzuS0(Kz, x)SPˆa3e−iK ∗ z z, (41) uSS(x, z)= eiγβ(x−x) dKzuS0(Kz, x)SSˆa4e−iKzz, where γβ = k2β− K2
z is the propagating wavenumber for S waves (here in the x-direction) and Kz ∗ =
k2
α− kβ2+ Kz2 is the transverse wavenumber of the
reflected P wave. Unit vectors ˆa3 = (γβ,−Kz ∗)/kα and ˆa4 = (Kz, γβ)/kβ.
|uS
0(Kz, x1)| is the scalar spectrum of the incident S wave with a transverse wavenumber Kz. PP, PS, SP and SS in Eqs.(38)–(41)are reflection coefficients
of different wave types at the free surface (Aki and Richards, 1980).Figure 15 is an example of those reflection coefficients versus horizontal slowness (ray pa-rameter p). InFig. 15, pA corresponds to P slowness (inverse velocity) and pB
to S slowness. For p < pA, P and S waves are both homogeneous waves, their
transverse wavenumbers are real. For p > pB, P and S waves are both
inhomo-geneous waves, their transverse wavenumbers are imaginary. For pA< p < pB,
P wave is inhomogeneous while S wave is homogeneous. A Rayleigh pole is lo-cated in the region of p > pB. In general, we can calculate all reflected waves
using Eqs.(38)–(41), once the incident fields|uP0| and |uS0| are known. However, numerically, it is more convenient to separate the calculation of Eqs.(38)–(41) into homogeneous and inhomogeneous waves, respectively.
For homogeneous waves, Eqs. (38) and (41) (common-type) can be imple-mented by FFT. However, the reflected waves of converted-type cannot be directly implemented by FFT because the nonlinear relationship exists between Kz and Kz∗ for P–S conversion (or Kz and Kz ∗ for S–P conversion). Although we can
obtain uniform samples with respect to Kzand Kz∗(or Kzand Kz ∗) by complex
FIG. 15. The free surface reflection coefficients (in logarithmic scale) versus horizontal slowness. The P and S wave velocities for the elastic half-space are α= 5 km/s and β = 3.5 km/s. The pAand
pSdenote P and S slownesses.
the noise due to the interpolation is so strong that the accumulated errors increase very fast for multiple step propagation. In our study, the direct summations over the incident waves (p < pA for P incidence or p < pB for S incidence) are
performed for calculating the converted reflections. Figure 16shows synthetic seismograms calculated with Eqs.(38)–(41) for an elastic half-space with only homogeneous waves. The results calculated with wavenumber integration (WI) method (dashed lines) are also shown as references. Since the source is deep compared with the propagation distance, Rayleigh wave is very weak in the ex-act solution.Figure 16a shows the vertical component of the displacement and Fig. 16b shows the horizontal component. From Fig. 16 we see that the cal-culations of the reflection and conversion by the free surface are in excellent agreement with the theory.Figure 17 shows synthetic seismograms for Flora– Asnes crustal model (see,Fig. 18) using elastic screen method. A double-couple source is located at a depth of 16 km and has a dominant frequency of 2 Hz. We see that both P and S waves are well excited.Figure 18is the corresponding snapshots. FromFigs. 17 and 18, the short-period phases Pn, Sn, Lg, etc., can be identified. For the elastic screen method at its current stage, only real trans-verse wavenumbers are used in FFT, which can only handle propagating waves (homogeneous waves).
(a)
(b)
FIG. 16. Synthetic seismograms calculated by the elastic screen method (solid lines) and wavenumber integration method (dashed lines) for an elastic halfspace. Only homogeneous waves are included in the results of elastic screen method. (a) Shows the vertical components of displace-ment; (b) Shows the horizontal components. A point explosion source is located at the depth of 16 km and the dominant frequency of source time function is 1 Hz. The first 4 receivers are placed along the free surface separated from the source by 100–124 km, and the last 5 receivers are placed in a vertical profile at an epicenter distance of 132 km and with depths ranging from 0–32 km.
For inhomogeneous waves, their transverse wavenumbers are imaginary so that Eqs.(38)–(41)cannot be calculated by FFT. However, the imaginary trans-verse wavenumber makes the propagation of inhomogeneous waves simple. The phase advance takes place only along the horizontal direction. It can be easily incorporated into the screen method, once the spectra of inhomogeneous waves
(a)
(b)
FIG. 17. Synthetic seismograms for Flora–Asnes crustal model (seeFig. 18) using P-SV elastic screen method. Only homogeneous wave are involved. (a) Shows the vertical components of displace-ment and (b) shows the horizontal components. A double-couple source is located at the depth of 16 km and has a dominant frequency of 2 Hz. Receivers are on the surface.
are known. Another important feature of inhomogeneous waves is the exponen-tial decay only in the direction perpendicular to propagation direction. Then the spectra of inhomogeneous waves can be calculated with Laplace transform. Fig-ure 19shows an example of such a treatment for Rayleigh wave propagating in homogeneous elastic half-space. The source is located at a depth of 2 km and has a dominant frequency of 0.5 Hz. The vertical receiver array is located at a distance of 100 km.Figure 19a shows the horizontal component of Rayleigh wave synthetic seismograms andFig. 19b shows the vertical component of syn-thetic seismograms. The solid lines are exact solution. The agreement between the screen calculation and the exact solution is excellent. The interaction between in-homogeneous waves and heterogeneities and the conversion between body waves and surface wave are still on-going research.
FIG. 18. Snapshots (horizontal component of displacement) for Flora–Asnes crustal model us-ing P-SV elastic screen method. A double-couple source is located at the depth of 16 km and has a dominant frequency of 2 Hz. The thicknesses of layers (from top to bottom) are 1 km, 14 km, 22 km and infinity, respectively. Their velocity and density parameters are α1 = 5.2 km/s, β1= 3 km/s, ρ1 = 2.6 g/cm3; α2 = 6.0 km/s, β2 = 3.46 km/s, ρ2 = 2.8 g/cm3; α3 = 6.51 km/s, β3 = 3.76 km/s, ρ3 = 3 g/cm3; α4 = 8.05 km/s, β4 = 4.65 km/s, ρ4 = 3.3 g/cm3. The ma-jor phases are labeled in the figure.
5. CONCLUSION
In the crustal waveguide environment, the major part of wave energy is car-ried by forward propagating waves, including forward scattered waves. Therefore, the neglect of backscattered waves in the modeling can still simulate the main features of regional waves in most cases. By neglecting backscattering in the theory, the method becomes a forward marching algorithm. A half-space screen propagator (generalized screen propagator) has been developed to accommodate the free-surface boundary condition and treat the SH wave propagation in com-plex crustal waveguides. The SH screen propagator has also been extended to the case of irregular surface topography by conformal or non-conformal topo-graphic transforms. For medium sized problems, the screen-propagator method is two to three orders of magnitude faster than the finite-difference meth-ods.
FIG. 19. Comparison of synthetic Rayleigh wave calculated using the screen method (dotted lines) with those calculated using the exact solution (solid lines). The source is located at a depth of 2 km and has a dominant frequency of 0.5 Hz. (a) Shows the horizontal components of displacement of Rayleigh wave and (b) shows the vertical components. The half-space parameters are α= 6 km/s and β= 3.5 km/s.
In the case of P-SV elastic screen propagators, plane wave reflection calcula-tions have been incorporated into the elastic screen method. Body waves including the reflected and converted waves can be calculated by real wavenumber inte-gration, while surface waves (Rayleigh waves) can be obtained with imaginary wavenumber integration. Numerical tests show good agreement with the the-ory.
From the theoretical developments and numerical tests of both SH and P-SV screen-propagators, we see that the one-way screen propagator approach for re-gional wave simulation is a viable approach and the savings in computation
