Prestack Parallel Modeling of Dispersive and Attenuative Medium

1 Institute of Geophysics, Department of Earth Sciences, National Central University, Chung-Li, Taiwan, ROC * Corresponding author address: Prof. How-Wei Chen, Institute of Geophysics, Department of Earth Sciences, National Central University, Chung-Li, Taiwan, ROC; E-mail: hwchen@earth.ncu.edu.tw This study presents an efficient parallelized staggered grid pseudospectral method for 2-D viscoacoustic seismic waveform modeling that runs in a highperformance multi-processor computer and an in-house developed PC cluster. Parallel simulation permits several processors to be used for solving a single large problem with a high computation to communication ratio. Thus, parallelizing the serial scheme effectively reduces the computation time. Computational results indicate a reasonably consistent parallel performance when using different FFTs in pseudospectral computations. Meanwhile, a virtually perfect linear speedup can be achieved in a distributed-memory multi-processor environment. Effectiveness of the proposed algorithm is demonstrated using synthetic examples by simulating multiple shot gathers consistent with field coordinates. For dispersive and attenuating media, the propagating wavefield possesses the observable differences in waveform, amplitude and travel-times. The resulting effects on seismic signals, such as the decreased amplitude because of intrinsic Q and temporal shift because of physical dispersion phenomena, can be analyzed quantitatively. Anelastic effects become more visible owing to cumulative propagation effects. Field data application is presented in simulating OBS wide-angle seismic marine data for deep crustal structure study. The fine details of deep crustal velocity and attenuation structures in the survey area can be resolved by comparing simulated waveforms with observed seismograms recorded at various distances. Parallel performance is analyzed through speedup and efficiency for a variety of computing platforms. Effective parallel implementation requires numerous independent CPU intensive sub-jobs with low latency and high bandwidth inter-processor communication. TAO, Vol. 17, No. 1, March 2006 204 (


INTRODUCTION
Realistic simulation of field data or pre-survey design frequently involves complex field acquisition parameters that require extensive computational resources and effective control of data flow.Viscoacoustic (Carcione 1988a, b;Chen 1996) and viscoelastic (Carcione 1988c, Carcione 1990;Robertsson et al. 1994) theories were previously implemented using pseudospectral (PS) or finite-difference (FD) schemes.Using the PS method for 3-D (Reshelf et al. 1988a, b;Chen and McMechan 1993;Furumura et al. 1998) or 2-D (Chen and McMechan 1992b;Chen 1996) seismic modeling and data processing (Chen and Chang 1999) is computationally and I/O intensive.Prestack modeling of seismic data for multiple shot and receiver geometry by executing a code sequentially on a uni-processor machine, makes it virtually impossible to achieve such goals.Therefore, efficient computer simulation of field surveys, and use of the simulated dataset for modeling, data processing, and analysis is essential.To compensate for long execution times (Fornberg 1987), parallel computing of wavefield simulations on a cluster of workstations (personal computers) or massive parallel processors (MPP) are alternative means of simulating full waveform data.
Effective parallel implementation requires numerous independent CPU intensive sub-jobs with low latency and high bandwidth inter-processor communication (Koshy et al. 1991).For a pre-defined large-scale model, the calculated synthetic responses for each shot point are independent of each other.Thus, the simulation of multiple shot seismic data has a useful coarse-grain level for parallelism.Versteeg et al. (1994) presented their parallel version of inversion method on a dedicated cluster of workstations.Liao and McMechan (1993) demonstrated their modeling scheme through a regular grid PS approximation with a machine dependent message passing interface on a specific Intel iPSC/860 machine.Hung and Forsyth (1998) simulate 3-D anisotropic response using a multidomain PS scheme.Herein, a numerical staggered grid PS method was employed and the algorithm was tested on various platforms, including a distributed-memory multicomputer system, a loosely coupled set of workstations and a cluster of PCs connected by a network.
In this study, I demonstrate the effectiveness of a wavebased parallel algorithm for simulation of full wavefield seismic data and also illustrate the diagnostic effects of anelasticity on synthetic seismograms.To achieve this, concurrent and scalable high performance computer simulations of 2-D anelastic responses are necessary for wide angle reflection and refractions.Finally, by simulating field data recorded on ocean bottom seismographs (OBS) acquired by the R/V Maurice Ewing during the TAICRUST survey in 1995, I demonstrate the potential feasibility of modeling of wide-aperture seismic data.

VISCOACOUSTIC MODELING
The key theory is based on Boltzmann's superposition of relaxation mechanisms to include viscoelasticity for various media.The loss mechanism, based on a standard linear solid model that includes a spring and a Kelvin element connected in series, is adapted for numerical simulation (Gross 1953;Christensen 1982).The equation of motion to simulate 2-D viscoacoustic responses in an isotropic linear viscoacoustic medium (Carcione et al. 1988a, b) is given by: where e denotes the dilatation, ρ represents the density, s x f is a source term given by the divergence of the i th component of body force per unit volume divided by the density, the dot above a variable denotes a time derivative, and p denotes the pressure.If the density ρ is constant, the equation (1a) then becomes: The general linear relation between pressure p and dilatation e in a viscoacoustic medium can be expressed as: in which Ψ c x z t ( , , ) represents the relaxation function (Liu et al. 1976; Carcione et al. 1988a, b) of the medium.The wave phenomenon associated with such a stress-strain relation is often expressed in the form of convolution integral based on Boltzmann's superposition principle.The kernel function of the integral is a memory function which describes the stress history dependence of strain.Thus, the Ψ c in (2) can be replaced by introducing memory variables.The relaxation function is: where stress ( τ σl ) and strain ( τ εl ) denote material relaxation times for each relaxation mecha- nism l.L denotes the number of relaxation mechanisms and M R represents the relaxed bulk modules of the material.Variables M R , τ εl and τ σl , are all spatial functions that characterize material properties.A material of this type is considered to have a memory because the current pressure depends on the full history of the stress and strain.Therefore, the relationship be-tween stress and strain is time dependent.The number of relaxation mechanisms and their corresponding relaxation times are chosen to provide the desired Q behavior as a function of frequency (ω ), where Q is expressed as the ratio of the real and imaginary parts of the complex modulus, that is: and the complex modules M c can be expressed as: If a set of strain and stress relaxation times is given, the related Q values can be obtained.The combination of relaxation times that produce any given Q( ) ω is clearly not unique; hence, many physical causes of viscoacoustic effects may behave similarly.
The phase velocity, defined as frequency devided by the real part of the complex wavenumber k M ω ρ 1 2 is the relaxed velocity of the medium.Source radiation simulation is generally expressed in terms of equivalent body forces (Burridge and Knopoff 1964).For a point pressure source approximation, with magnitude G acting at a point r s in the r direction, the body force term (Eq.1b) can be further derived from a scalar potential φ( , ) x z s s by differentiation.By means of the band-limited Gaussian-type distribution of Kosloff et al. (1984), φ( , ) x z s s is introduced over a small region of the grid given by: where r x z s s s ( , ) denotes the source location, and the δ( ) r r s − is an abbreviation for spatial Dirac delta function δ acting along x-and z-directions.Superposition of individual components and their derivatives allows for the generation of a variety of sources in elastic (Chen 1992a; Fig. 1) and viscoacoustic media (Chen 1996; Fig. 7).

MODEL CONSTRUCTION AND STAGGERED GRID COMPUTATION
Viscoacoustic theory holds that the model can be parameterized by appropriately combining strain ( τ εl ) and stress ( τ σl ) relaxation times, density, relaxed and unrelaxed moduli, and velocity.Attenuated and de-coupled propagating P-or SH-wavefields are synthesized by constructing P-and SH-velocity, ρ, Qp, and Qs models separately.A FD equation from wave equation ( 1) is derived based on numerical operators that approximate the partial derivatives both in time and space.Meanwhile, a staggered grid scheme removes non-physical scattered waves generated by "steps", caused by the discretization of the model (Chen 1996).That is, the pressure is defined at each sampling location, whereas pressure gradients are evaluated half-way between sampled positions (Fig. 1).Zeng and West (1996) proposed an alternative approach, using a spatial-averaging operator to reduce spurious diffractions in the calculated wavefields.The staggered grid PS method is highly desirable for modeling materials with moderate to highly contrasting Poisson's ratios, such as boundaries between fluids and solids (namely, land and ocean, or boreholes) or a boundary between a salt dome and surrounding country rock.According to this scheme, a widely varying Q can be modeled in a highly heterogeneous medium.
The solution of the wave equation [Eq.( 1)] is obtained by PS computation with the finitedifference leap-frog scheme.The even-based Fast Fourier transform (FFT) with a staggered grid is implemented to approximate spatial derivatives (Chen 1996).Meanwhile, temporal derivatives are approximated by second-order central differences.The leap-frog time marching scheme is applied such that an updated wavefield is obtained from wavefields at the previous two time steps (Chen 1996).Wraparound and edge reflections are suppressed through a composite absorbing boundary mechanism by systematically reducing the wave amplitude in a strip around the numerical mesh (Kosloff and Kosloff 1986;Sochacki et al. 1987, Chen andHuang 1998).A vanishing pressure, p = 0, is applied at the top of the model as the free-surface boundary condition.Thus, viscoacoustic wave propagation in a dispersive and attenuative medium is accurately described by implementing a staggered grid PS technique.

CONCURRENT FORWARD MODELING
Concurrent forward simulation in a distributed-memory multi-processor computer is implemented.During a seismic field experiment, data are collected by common-shot configurations along a predetermined seismic survey line.Thus, a two-dimensional flat earth that disperses and attenuates propagating waves is considered.The scenario is to mimic a wideangle refraction or a reflection survey for seismic exploration or to generate synthetic seismograms from various earthquake sources for earthquake research.The parallel strategy is based on partitioning the computation domain, and the field quantities are assigned to a number of processors.Each processor either shares the same earth model with each other or has its own, simulating different responses for different source locations.The characteristics of viscoacoustic simulation of seismic data are presented through two realistic applications.A synthetic example illustrates the main features of anelastic seismic responses in viscoacoustic wave simulation, while a field data application illustrates the feasibility of modeling wide-angle data and its potential for large scale imaging of deep crustal structure.

Strong Heterogenity and Complex Structure: a Synthetic Application
Modeling and imaging of complex structures around and beneath salt diapirs, domes, and sheet structures with the potential to trap significant hydrocarbon accumulations has been considered challenging.A high velocity, low density, high Q salt body flowing upwards under geopressure attributed to the weight of overlying sediments creates a highly complex geological structure.A salt diapir structure provides stratigraphic traps for hydrocarbon accumulation around the sides of a salt dome.Full-wavefield modeling for pre-survey design can provide the necessary acquisition parameters to ensure cost-effective field experiments that collect the most useful data necessary to reveal the complex geologic features around and under salt bodies.Thus, the proposed simulation strategy can reduce the risks associated with sub-salt drilling and help optimize well locations.
To illustrate how viscoacoustic media respond to external disturbances from various source types, and to simulate seismic responses concurrently, a realistic salt-tongue model that contains a highly heterogeneous, laterally varying velocity and Q model with a steep salt flank (Fig. 2) is used.Table 1 presents the specified quality factors derived from two stress and strain relaxation times (τ σl , τ εl ), relaxed moduli ( M R ), density ρ, and the acoustic velocity for each layer.The relaxed velocity and maximum phase velocity were derived from Eq. ( 6).For simplicity, the densities of the media are all the same (2.0 g cm −3 ) and the second stress (and strain) relaxation time is defined as one-tenth of the first one.For further simplicity, the boundaries of various Q structures are assumed to be coincident with the boundaries of velocity changes, though the Q model can be specified independently of the velocity model.
A complete set of survey data was simulated by a group of impact sources and its corresponding common-source gather.Table 2 lists the related parameters, including source and receiver locations, for conducting a seismic survey over a salt diapir structure.The radiation pattern of an impulsive source is mimicked by a vertically oriented single force (cf.Chen 1996; Fig. 7). Figure 3   shown in gray scale.The number indicated within each layer corresponds to each medium listed in Table 1.Each medium has different material property contenting fixed Qp and corresponding variable phase velocity structure.
Table 1.Physical material parameters used in the model.Two relaxation times for each medium; quality factor; unrelaxed modules; density and pure (acoustic) P-wave velocity for a salt-tongue model in Fig. 2 are calculated according to the theory of viscoelasticity.Phase velocity for attentive media can be derived from relaxation times.
attenuation curves ( Q −1 ) varies over a diverse range of frequency (ω ) values.The property of each layer is described by a relaxed modulus M R for velocity and by two relaxation mechanisms for attenuation that yield two minimum values in the plot of quality factor versus frequency.Snapshots of the acoustic (infinite Q, Fig. 4a) and viscoacoustic (variable Q, Fig. 4b) wavefields are compared at the same times.In a geological environment where the sediment velocity increases with depth, turning waves and salt bottom reflections can be recorded in a wide-aperture seismic experiment.In Figs.4a and b, noticeable salt bottom reflections (denoted BR) are present in both acoustic and viscoacoustic cases.The effectiveness of the composite absorbing boundary condition can be also visualized from snapshots.Those unreal reflections, Table 2. Survey parameters of a salt-tongue model.such as spurious reflections and wraparounds, from artificial boundary of the finite numerical space are sufficiently suppressed.Figures 5a and b illustrate the corresponding representative common-source gathers from finite aperture geophone arrays located at the surface of the model.The computed synthetic seismograms contain prominent reflections from both the sediments and salt flanks.A noticeable amplitude decay is more pronounced for the viscoacoustic wavefield simulation.
Figures 6a and b compare viscoacoustic and acoustic time histories, along with their corresponding amplitude spectra, at a fixed surface recording station (x = 9.0 km) above the salt, with shot location at 8.61 and 2.85 km, respectively.In both figures, both time histories contain the direct waves, and multiply reflected arrivals from the salt body and surrounding layered structures.For the shot located directly above the salt dome, the recorded data (Fig. 6a) Fig. 3.The Qp-ω curves obtained from model parameters presented in Table 1.
The variable Qp is specified within each layer denoted by the same number shown in Fig. 2.  2, Table 1 and Table 2, respectively.All panels have the same amplitude scaling factor.Salt bottom reflection (BR) may travel downward first then turn upward before being reflected at the salt flank (Fig. 2) and returning to the surface.exhibits relatively small amplitude differences between direct acoustic and viscoacoustic arrivals.However, time accentuates the difference in the latter seismic phases.For the shot located far from the salt dome that has the same receiver position (Fig. 6b), anelasticity exerts more significant effects on arrival time and amplitudes in the recorded time history.Because anelastic effects on wavefields are cumulative along the propagation journey, the resulting attenuation and dispersion, revealed in reduced amplitudes and time-shifted phase arrivals, are more visible for the synthetic records of Fig. 6b, in which the wave travels through the same structure but over a longer distance, than they are for the records of Fig. 6a.The corresponding amplitude spectra of the recorded waveforms also exhibit similar behavior; that is, the high frequency components are attenuated more significantly.As travel distance increases, the decay of low-frequency components and overall amplitude becomes more prominent.
The attenuation effects due to intrinsic Q are apparent in both the time and frequency content of the recorded signal.Comparison of the time records in Fig. 6 indicates that the viscoacoustic pulse arrives before the acoustic one.This early arrival is caused by physical dispersion, in which the unrelaxed phase velocity (therefore group velocity) exceeds the relaxed phase (acoustic) velocity for the entire frequency range.To clarify and verify this observation regarding early arrival, the theoretical difference in phase velocity between viscoacoustic and acoustic dispersion effects for each specified layer (Table 1) is calculated from Eq. ( 6) and compared in Fig. 7.Although the maximum difference in phase velocity is approximately 1.4% for a 13 Hz source signal, the cumulative effects on travel time and waveforms of the recorded signals are obvious.Thus, realistic description of attenuation effects for field data simulation is necessary.

Wide-Angle Deep Crustal Study: a Field Data Application
The ultimate goal of developing viscoacoustic wavefield modeling is the application of simulated waveforms to the acquisition and interpretation of field data.The deep seismic refraction and reflection profiles to be modeled herein were collected from the TAICRUST survey by the R/V Maurice Ewing off eastern and southern Taiwan in 1995.The TAICRUST marine seismic survey was aimed at imaging subduction-collision tectonic features and investigating the influence of the Taiwan arc-continent collision at the termination of the Ryukyu subduction system and the deformation of the Philippine Sea plate (Lallemand and Liu 1997).Line 1 (EW9509-01) was located west of Gagua Ridge and ran roughly perpendicular to the Ryukyu Trench.Eight common-receiver OBS profiles along Line 1 (Fig. 8), crossing the Ryukyu trench-arc-backarc system, were acquired in the area off eastern Taiwan.Line 1 comprises 9422 shots with average shot spacing of 34 m and total length of 325 km.Only a representative field record is used to demonstrate the feasibility of PS modeling.
To model wide-angle seismic data, parallel simulations of viscoacoustic wave propagation are performed in the velocity and attenuation model suitable for the study region.The results from previous works of ray modeling (Chen 1990) and traveltime inversion (Wang 1996) provide low-resolution but reasonable initial velocity distributions and those from spectral decay analysis (Chen 1993) give Qp structure distributions.Examination of shallow crustal structures from the near-trace seismic section (Schnurle et al. 1988) provides additional constraints for the inferred velocity discontinuity in the model.Meanwhile, the stacking velocities obtained from velocity analysis are employed to construct the near-surface velocity model.Table 3 lists the estimated intrinsic attenuation derived from strain and stress relaxation times.The 325 km long model was further divided into 8 sub-domains with model dimensions of 80 km × 40 km.Each sub-model is discretized into 2048 × 1024 grid points along the horizontal (position) and vertical (depth) directions.The source dominant frequency is set at 10 Hz to reveal the main features in the synthetic seismograms.After carefully considering the run time and stability conditions required for numerical simulation, the modeling is performed mainly on an IBM SP2 computer.For effective computer simulation, an explosive source is Fig. 7. Theoretical phase velocity difference predicted from the viscoelastic theory [Eq.( 6)].The numbers next to individual curves correspond to the media listed in Table I.The maximum difference in phase velocity for a given model is approximately 1.4% higher than that in the purely acoustic model at the dominant frequency of 13 Hz (dashed line).The effects on traveltime, waveform and amplitude in the recorded signal are obvious both in time and spectral analyses (Fig. 6) and in phase velocity-frequency analysis.
Fig. 8.The main subsurface structure features along profile Line 1 (EW9509-01).This profile runs N-S across the Ryukyu arc-trench and backarc system.Shallow (upper panel) and whole models (lower panel) show the complexity of the structure along the northwestern convergent region.Vertical exaggerations of 8.0 and 4.0 times are applied in the upper and lower panels, respectively.A portion of the synthetic seismic profiles in Fig. 9 for OBS 8 is marked by dashed lines.Stars denote OBS locations along line 1.The numbers denote the Q values in each medium (Table 3).
(OT: Okinawa Trough; RA: Ryukyu Arc; RAS: Ryukyu Arc Slope; NB: Nanao Basin; YR: Yaeyama Ridge; RT: Ryukyu Trench; TC: Taitung Canyon; HB: Huatung Basin).located at the center of the model and on the ocean floor.1988 receivers are located at the two grid points below the sea surface.
Complete sets of prestack synthetic seismograms are obtained by simultaneous computation over eight processors and compared with field recorded data.Both acoustic and viscoacoustic wavefield responses are constructed and compared to estimate the combined Table 3. Physical attenuation parameters used in the model.Stress (τ σl ) and strain ( τ εl ) relaxation times for each relaxation mechanism l. and their corresponding Q values of the material.
effects from attenuation and spreading loss.Because of excessive outputs, only a minimum amount of data for the interpretation of tectonic structure features is presented herein to illustrate the feasibility of prestack forward modeling.OBS station 8 located at the center of the Nanao Basin and associated with two other cross-line profiles is crucial for revealing a typical convergent margin structure in trench, accretionary wedge, forearc, arc, and backarc basins (Fig. 8).The Nanao Basin is known to be at the center of three sequential step-wise elevated forearc basin structures as the Ryukyu arc extends toward the Taiwan collision zone.Therefore, particular phases from OBS 8 must be structurally interpreted and identified for the entire survey.According to our experience, model building and updates as well as interpretation and delineation of subsurface interfaces are more time consuming than performing massive parallel computer simulation.
The OBS raw data from station 8 (on Nanao Basin) in Fig. 9a exhibit earlier arrivals on the north side, implying a shallower sea-floor and/or a higher apparent velocity north of the station.Faster and weaker refractions on the north side branch arise from the bathymetric and low transmission effects of the Ryukyu Arc Slope (RAS in Fig. 8).The seismic reflections of the forearc Nanao basin strata, up to 2.5 sec, or about 3 km thick, are strong and continuous.The representative viscoacoustic synthetic records of OBS 8 (Fig. 9b) exhibit similar seismic responses.Major features, including first arrivals, crustal refractions, head waves, inter-crustal reflections, reflections from the bottom of the sedimentary layer, diffractions, water-bottom and peg-leg multiples, are evident in both the field and synthetic profiles, indicating those features in Figs.9a and b can be helpful in tracing differences and similarities between synthetic and real data.Simultaneous modeling of arrivals recorded at all stations through waveform modeling and cross checking with ray tracing results can help provide reasonable constrains; besides which, consistencies in the travel time and relative amplitude variations between field and synthetic data further support the model.Certain discrepancies may be attributed to uncertain velocity, Q, and source signature.The computational wall clock time of a complete run over 8 processors on an IBM SP2 is approximately 27 hours.Viscoacoustic simulation usually takes no more than 10% more of elapsed time than acoustic computation.Therefore, concurrent forward modeling of large scale, wide-angle seismic field data is necessary and feasible.

PARALLELISM AND PERFORMANCE ANALYSIS
In full waveform modeling, a group of shots and associated parameters is distributed to each high-performance processor to simulate seismic responses concurrently.Our parallel implementation uses a master-multislave strategy.One processor is assigned as a master.The master spawns and directs slave processors that carry out most of the computations.The master processor also handles data input and output.The message passing routines provide communication among processors via a message passing interface called the Parallel Virtual Machine (PVM).Performance can be evaluated by measuring increases in speedup and efficiency (Appendix A).Implementation with a PVM software system simplifies network parallel computing without requiring detailed knowledge of how to work with a heterogeneous computer network.Because the system has distributed memory, the time required for inter-processor communication is significant.Owing to the nature of the FD leap-frog scheme used in wavefield simulation, seismograms are stored in a temporary array and written to disk at specified time intervals.The outputs, including seismograms and snapshots, are written to a shared disk via a high-performance switch (on IBM SP2) during the computation.The current approach facilitates a more realistic measurement of throughput performance (Appendix B) with any given system configuration, including I/O operations.
For 2-D seismic applications, the entire data set can reside in the CPU memory, simplifying the implementation of network parallel computing.The analysis of parallel efficiency herein indicates that an almost perfectly linear speedup can be expected under the dedicated mode (Fig. B-1a).The excellent performance-to-cost ratio makes parallel computing extremely attractive.Therefore, parallel computing with a distributed-memory multiprocessor system is a viable means of implementing 2-D prestack modeling.
The novel algorithm was also tested on network-connected workstations and on a cluster of Unix-based personal computers.However, according to those results, the computational capacity is limited by the local computational environment, including network bandwidth, wiring, heterogeneous CPU speed and time sharing (Fig. B-2).The requirements for costeffective approaches and practical applications remain the primary considerations in designing parallel algorithms for both 2-D and 3-D computation.Our parallel algorithm provides an alternative and economical approach that can be run on a PC or workstation cluster when no parallel computer is available.For optimum computational performance, a custom-made massively parallel computer should always be considered first.The optimum parallel paradigm is an approach that minimizes both time and cost requirements.

CONCLUSION
This study shows how numerical modeling of viscoacoustic seismic responses for 2-D media is developed, tested, and implemented in a network distributed parallel computing environment.Implementation of viscoelasticity provides a highly effective means of simulating realistic geological structures.The numerical scheme proposed herein includes the typical parameters of physical materials, such as density, velocity and intrinsic Q, to simulate attenuated, dispersive material.Considerable anelastic effects on reduced amplitude due to intrinsic Q structure and temporal shift because of physical material dispersion can be modeled and analyzed quantitatively.Staggered-grid viscoacoustic parallel computing becomes more realistic when a large-scale geological model is used and increasingly attractive when limited computer memory is available for waveform modeling.
Synthetic data examples demonstrate that the parallel implementation of prestack modeling of viscoacoustic wave propagation for various geophysical problems is feasible and highly effective.With the advent of multi-processor computer technology, most seismic forward and inverse processing can be easily decomposed into shot domain jobs and distributed amount several individual processors.The proposed algorithm permits the use of several nodes to solve a single large problem at high computation to communication ratios.This scheme is particularly useful and efficient when many processors share the same task, such as prestack seismic modeling of common source gathers.If the model is too large to fit into core memory, then domain-decomposition is an alternative approach for parallelism.
The capacity to perform efficient parallel computing for data interpretation and for subsurface structure delineation through modeling of wide-angle OBS seismic data is also presented.By integrating the combined technologies of 2-D seismic acquisition, 2-D seismic depth processing (Chen and Chang 1999), and an efficient wave-based parallel simulation strategy, we have been able to reconstruct many of the observable seismic signatures collected from the wide-angle refraction and reflection survey of the TAICRUST project off east Taiwan.Detailed discussion and interpretation of seismic field data in terms of their tectonic implications are beyond the scope of this study and will be presented elsewhere.Results presented herein provide a valuable reference for efforts underway to develop a parallel forward modeling and prestack seismic data processing system for geophysical surveys conducted in Taiwan.
Future investigations will explore the possibility of frequency-dependent Q behavior for diverse types of rocks.Developing a viscoelastic wavefield modeling algorithm that accounts for the propagation of coupled S-and surface waves is the next step in confirming the prediction of recorded seismic data.Extension to 3-D viscoacoustic and viscoelastic forward simulation is also necessary to identify the effects from 3-D structure.In terms of 3-D application, the volume of the dataset can be too large to accommodate available memory.Hence, a strategy that distributes data among several disks with parallel disk I/O control and high performance dynamic networking switch may be necessary.the single program multiple data (SPMD) model.In this approach, the application programmers are responsible for the choice of loops and distribution of loop iterations.Partitioning of data arrays over processors can be handled within the application.The parallel algorithm was tested on a cluster of various SUN workstations and an MPP machine.The later was an IBM 9076 SP2 computer, consisting of thirty-two RISC System/6000 processor nodes.
The performance can be assessed in terms of speedup and efficiency measurement.The speedup factor is defined as S T T serial parallel = / , where T serial denotes the run time for a single processor to complete a sequential task, and T parallel represents the run time for multi-proces- sors to complete the same task in parallel.For the convenience of analysis, we compare the parallel run time with one slave to parallel run time with multi-slaves as a measure of the speedup factor and the corresponding efficiency.The parallel run time is T N p ( ), where N p denotes the number of slave processors used.The speedup factor as a function of N p is given by S N T T N p p ( ) ( )/ ( ) = 1 and the efficiency (E) can be defined as E N S N p p ( ) / = . where E = 1 represents a perfect parallel performance, and 0 < E < 1 refers to an imperfectly parallel performance.
Table B-1.List of relevant FFT kernels used for PS wavefield simulation.(a) A comparison is made of average run times used from three different viscoacoustic wavefield simulators.The FFT kernels include Convex vector library routines ("Cldfft" and "C2dfft") and a canned routine "nlogn" from Robinson (1967).The nearly linear speedup and efficiency can be achieved as the number of processors is increased (dashed line and squares).The parallel performances drop (solid lines with cross and circle symbols) for more than four processors is mainly attributed to the narrow bandwidth of the computer network.
the narrow bandwidth of such computer networks.Because the calculated wavefield responses in our computations are output to a shared disk at regular time intervals, network communication is necessary.The efficiency performance in the lower panel indicates that inter-processor communication and scheduling can markedly reduce computation efficiency even though inter-processor communication has been reduced to a minimum.Perhaps hostess style parallelism rather than the host-slave style can be adopted for this purpose.Further investigation of inter-processor communication for real data applications and for various vector lengths is essential for large-scale problems.

Fig. 1 .
Fig. 1.A spatially staggered grid system to model acoustic and viscoacoustic waves in two-dimensional Cartesian coordinates.
depicts the Qp versusω curves that show the behavior of frequencydependent P-wave quality factors associated with individual velocity layers in the salt-tongue model for a specific number of relaxation mechanisms [L = 2 in Eqs (3) and (5)].A family of

Fig. 2 .
Fig. 2. Salt-tongue model.Variation of compressional wave-speed distributionshown in gray scale.The number indicated within each layer corresponds to each medium listed in Table1.Each medium has different material property contenting fixed Qp and corresponding variable phase velocity structure.

Fig. 4a .
Fig. 4a.Fixed-time wavefield snapshots generated from a salt-tongue model.Variations of seismic responses for purely acoustic (infinite Q) conditions are clearly discernible.From top to bottom, acoustic wave propagation through various media is shown at 0.75, 1.5, 2.25 and 3.0 seconds, compared with the viscoacoustic responses in Fig. 4b.Model and survey parameters are shown in Fig.2, Table1 and Table 2, respectively.All panels have the same amplitude scaling factor.Salt bottom reflection (BR) may travel downward first then turn upward before being reflected at the salt flank (Fig.2) and returning to the surface.

Fig. 5 .
Fig. 5. Synthetic common-shot gathers calculated from purely acoustic (5a) and viscoacoustic (5b) wave modeling (Fig. 2).Seismograms are recorded along a fixed cable length.A large survey aperture and appropriate recording time are the prerequisites for successful salt and subsalt imaging.Purely acoustic and viscoacoustic responses are displayed for comparison.

Fig. 6 .
Fig. 6.Comparison of attenuation and dispersion effects between viscoacousticand acoustic responses for a fixed channel at different shot locations.The temporal history and their corresponding amplitude spectra illustrate the fundamental difference in viscoacoustic (solid line) and acoustic (dashed line) responses for a source location at 8.61 (a) and 2.85 km (b), respectively.The receiver location is above the salt at 9.0 km.

Fig. 9 .
Fig. 9. Prestack modeling of seismic data collected at OBS 8 station along Line 1 (Fig.8).The processed field seismic data on the vertical component (9a) and synthetic viscoacoustic responses (9b) are plotted in reduced time format with a velocity of 8.0 km s 1 − .A different distance dependent, amplitude correction is applied to compensate the discrepancy resulting from different geometrical spreading between the real (3-D effects) and synthetic (2-D effects) data.

Fig. A- 1 .
Fig. A-1.The PVM system designed to permit a network of heterogeneous UNIX computers used as a single parallel computer (reproduced from Konuru, Otto and Walpole 1997).Local and remote message (Msg.)communications among processes (Tn) are controlled by the daemon (PVMD) and a run-time library (PVMLIB).
Fig. B-1.Timing information of sequential and parallel tasks run on an IBM 9076 SP2 computer consisting of thirty-two RISC System/6000 processors.Timing analysis results are compared with their average elapse time at one time step through a complete execution of the program.The timing information under dedicated mode (B-1a) and time-sharing mode through simultaneous tests (B-1b) are compared for different vector lengths.The results from simultaneous tests, which mimic multiple jobs running on the SP2, indicate the difference in throughput performance between original and vectorized code.

Fig. B- 2 .
Fig. and efficiency (lower) curves under dedicated and timesharing situations.The nearly linear speedup and efficiency can be achieved as the number of processors is increased (dashed line and squares).The parallel performances drop (solid lines with cross and circle symbols) for more than four processors is mainly attributed to the narrow bandwidth of the computer network.