S34B-01
Computation of Travel Time Through 3D Velocity Models for Applications in Real-Time, Global Seismic Event Monitoring
Three dimensional velocity models of the Earth have been little used by real-time global monitoring agencies despite the expectation that these models might improve the accuracy and reduce the uncertainty of the seismic event locations they calculate. There are many reasons for this reluctance to adopt 3D models, including 1) uncertainty that adoption of 3D models will in fact significantly improve locations, 2) questions about how to quantify the uncertainty of the travel time predictions, 3) uncertainty as to how to assess the fidelity of computed travel times relative to the input velocity model, 4) questions about the computational architecture most appropriate for calculating predicted travel times, and 5) concern that the computational cost of calculating travel time predictions in a real time monitoring environment will be prohibitive. In this paper we begin to address the implications of using 3D velocity models in real-time global monitoring environments by addressing the last 3 items in the list above, which focus on the computational aspects of using 3D velocity models for travel time prediction. There are three fundamental approaches to computing travel times through 3D velocity models: 1) fix a source location at some position in the Earth model and compute travel times to all nodes in a 3D grid of nodes surrounding the source locations by solving the eikonal equation, 2) fix the locations of a single source and single receiver within the 3D velocity model and find the ray path(s) that honor Snell's Law in between (boundary value problem; ray bending), and 3) fix a single source location within the model and iteratively modify an initial estimate of the ray parameter searching for a ray that arrives at the receiver (initial value problem; ray shooting). To assess the computational issues with the use of these types of travel time calculators we have implemented the Fast Marching Method of de Kool, et al (2006), which is an eikonal solver, and the pseudo-bending algorithm of Um and Thurber (1987). In this paper, we compare the relative merits of these approaches in the context of their use in a real-time global monitoring environment.
S34B-02
Extended Diffraction Tomography
We present the development of extended diffraction tomography, a new approach to the solution of the linear seismic waveform inversion problem for exploration-style acquisition geometries. This method has several appealing features, such as the use of arbitrary depth-dependent reference models and the decomposition of the full 2D or 3D inverse problem into a large number of independent 1D problems. This decomposition makes the method naturally highly parallelizable. Careful implementation yields significant robustness with respect to noise. Several synthetic examples are shown which characterize the benefits of our method and demonstrate the usefulness of choosing realistic 1D reference media.
S34B-03
Applying Waveform Tomography to Refraction Seismic Data - Inversion Strategies and Resolving Power
A number of studies on waveform inversion of synthetic refraction seismic data have proven its capability to resolve complex subsurface structure at sub-wavelength scale. However, these data often have an unrealistic bandwidth, were generated with unrealistic models (e.g. lacking attenuation), assuming optimal geometry (e.g. no topography), and unrealistic signal-to-noise ratio (mostly ∞). This study investigates the practical limitations of resolution when applying frequency domain full waveform inversion to exploration scale surface refraction data, and the possibilities to push these limits through appropriate inversion strategies. Instead of computing and inverting more realistic synthetic data --- certainly a valid approach in order to bridge the gap to application --- actual field data from two different surveys were inverted, which are similar in acquisition geometry but differ strongly in terms geological heterogeneity. Exploring the parameter space of these inversions allows for evaluating the importance of the various factors that affect the results (the background model, the data properties and preconditioning, the forward modeling, the inversion strategy…). The most crucial problem in this context is the variability of observed signal amplitudes: Near-surface layers often exhibit strong and strongly varying attenuation due to weathering, different levels of lithification, different porosity and the like, which may, in addition to receiver coupling, alter the signal amplitude on the order of magnitudes. When attempting to minimize an objective function that has been posed as the (squared) residual of recorded and computed seismograms --- the most common approach in waveform inversion --- the use of true amplitudes is most likely to fail in the face of non-negligible, but unknown, attenuation variations. Finding a strategy to address this problem is a requirement for a successful application. Attempting to simultaneously reconstruct an attenuation model (along with velocities) does not alleviate this problem. Using normalized amplitudes to reconstruct velocities and logarithmic amplitudes to subsequently image attenuation is a much more promising approach. As a consequence of amplitude normalization, the introduction of weighting factors becomes essential in order to mitigate the impact of noise in the data. Another problem that effectively limits the resolution for one of the surveys is the lack of a surface in the model despite strong elevation variations.
S34B-04
High-resolution P and S-wave Velocity Structures from Elastic Full Waveform Inversion of Multi-Component Ocean Bottom Cable Seismic Data
Full waveform inversion is becoming a realistic option with the advent of modern computing facilities, both in global and exploration seismology. Over the last ten years, we have developed a series of elastic full waveform inversion algorithm and have applied to a variety of acquisition geometry. The forward modelling is based on the finite difference approximation to the full elastic wave equation in the time domain, which can incorporate converted waves, refraction, and attenuation. The inversion algorithm is based on the minimisation of observed data with synthetic data in a least-squares sense, and requires a cross-correlation of the back propagation of residual with forward propagated wavefield in a background media. Starting with the background velocity obtained using travel time inversion, we first invert wide-angle and low frequency data, which provides medium wavelength velocity structure, and then invert near offset and high frequencies that leads to high-resolution P- and S-wave velocity structure. We first invert vertical component data to obtain short wavelength P- and S-wave velocities, which are constrained by amplitude versus offset behaviour of the P-P reflection, and then invert horizontal component data to obtain very-high resolution S-wave velocity structure, which is constrained by P-S reflection. Finally, we invert all the data simultaneously to have consistency over the data and model space. We found that the high-resolution S-wave velocity image is far superior than the P-wave velocity image and provides information that may not be present in the P-wave velocity image. Combined P and S-wave velocity structure could be used to quantify sub-surface lithology and fluid saturation and pressure. In this presentation we will highlight the challenges faced during the development of our waveform inversion and their implication for the global seismology problems.
S34B-05
On Characterization of Elasticity Parameters in Context of Measurement Errors
In this presentation, we discuss the one-to-one relation between the elasticity parameters and the traveltime and polarization of a propagating signal in the context of the measurement errors. The one-to-one relationship between seismic measurements and a model postulated in the realm of the constitutive equation of an elastic continuum provides the link between the observational and theoretical aspects of seismic tomography [1]. The existence of this link encourages us to develop methods of inferring the elasticity parameters from measurements. However, a consideration of required accuracy and the analysis of error sensitivity suggest that the pragmatic application of this one-to-one relationship might be a difficult task indeed [4]. There are eight symmetry classes of an elastic continuum whose properties are contained in the density-scaled elasticity tensor [6]. Given this tensor in an arbitrary coordinate system, we can identify to which symmetry class it belongs, as well as obtain the orientation of its symmetry axes and planes, and hence the elasticity parameters in a natural coordinate system [2]. To obtain the tensor to be studied, we consider either ray velocities and polarizations [1] or wavefront slownesses and polarizations [5]. For the former, we assume that the medium is homogeneous in order to invoke the straightness of rays to calculate ray velocity given the source and receiver position; for the latter, we assume that the medium is homogeneous in at least one direction in order to invoke the ray parameter. In spite of the limitations due to homogeneities, both approaches are sensitive to measurement errors, which are not negligible. In view of these observational concerns [4], we consider several weaker objectives based on the theoretical formulation. Rather than distinguishing among eight symmetry classes and obtaining the corresponding elasticity parameters, we might be able to distinguish among a few groups that contain several classes within them and are characterized by ranges of parameters. Such an approach takes advantage of similarities among several eigenproperties that distinguish a given group from the others. Furthermore, we might not require to measure the traveltime of the three waves --- the quasishear wave being more difficult to observe. Also, we might not require to measure polarizations, which, in general, exhibit a larger measurement error than do the traveltimes. (To obtain a complete elasticity tensor we need both polarizations and traveltimes for the three waves [3].) 1. Bóna, A., Bucataru, I., Slawinski, M.A. (2007) Elasticity parameters from traveltime and polarization measurements. Journal of Applied Geophysics (accepted) 2. Bóna, A., Bucataru, I., Slawinski, M.A. (2007) Coordinate-free characterization of elasticity tensor. Journal of Elasticity 87(2-3), 109--132 3. Bóna, A., Bucataru, I., Slawinski, M.A. (2007) Material symmetries versus wavefront symmetries. Q. Jl Mech. appl. Math 60(2), 73--8 4. Bóna, A., Slawinski, M.A. (2007) Comparison of two inversions for elasticity tensor. Journal of Applied Geophysics (submitted) 5. Dewangan, P., Grechka, V. (2003) Inversion of multicomponent, multiazimuth, walkaway VSP data for the stiffness tensor. Geophysics 68(3), 1022--1031 6. Ting, T.C.T. (2003) Generalized Cowin-Mehrabadi theorems and a direct proof that the number of linear elastic symmetries is eight. Internat. J. of Solids and Structures 40, 7129--7142
S34B-06
New 3D Vs model of Europe and the Mediterranean Basin
We invert a dense dataset of teleseismic and regional phase velocity observations for global 3-dimensional radially anisotropic shear velocity structure of the upper mantle. Our algorithm finds the least squares solution of the linear inverse problem by Cholesky factorization of the ATA matrix. We parameterize the model with a series of latitudinal, longitudinal, and radial splines. Multiple resolution parameterization allows us to discretize Europe and the Mediterranean with a higher density of splines than elsewhere, thereby utilizing the dense data coverage within the region. By using a this parameterization, we decrease the risk of erroneously mapping external anomalies into our high-density region. We calculate our sensitivity kernels iteratively, starting with 1D kernels based on a combined Crust 2.1 and PREM starting model. The velocity model resulting from inversion with the first kernels is then used to recalculate kernels for the subsequent inversion. We will present the resulting SV and SH models. Key advancements over previous models include increased resolution and improved imaging of key tectonic features such as continuous fast velocities below the Aegean Arc and Alpine Orogeny.
S34B-07
LOCAL EARTHQUAKES TOMOGRAPHY IN THE SOUTHERN TYRRHENIAN REGION (ITALY): GEOPHYSICAL AND PETROLOGICAL INFERENCES ON SUBDUCTING LITHOSPHERE
The Calabrian Arc, Southern Italy, is characterised by the subduction of the Ionian lithosphere -since Middle Miocene- beneath the Tyrrhenian basin. The related Benioff zone is seismically active to a depth > 500 km. The tomoDD code [Zhang and Thurber, 2003] was adopted to perform the tomography, using a set of 2463 earthquakes located in the window 14°30' E - 17°E and 37°N - 41°N, and recorded by seismic networks of the INGV in the period 1981-2005. Several inversions were performed using different selections of absolute and differential data obtained varying the maximum RMS and the threshold of the inter-event distance. Various synthetic and experimental tests were executed to evaluate the resolution and stability of the tomographic inversion. The inversions carried out for the synthetic and the restoration-resolution test [Zhao et al., 1992] were repeated several times with the same procedure used in the inversion of experimental data. The lack of bias in the models, related to the different grid- node positions, was tested performing inversions rotating, translating and deforming the original grid. To evaluate the dependence on the initial model, several inversions were also done using different 1D and 3D models simulating slab features. Finally, 35 models resulting from the inversions were synthesized in an average model obtained by interpolating each velocity model into a fixed grid. Each velocity value interpolated was weighted with a corresponding DWS (Derivative Weight Sum) resulting thus a Weighted Average Velocity model. The highly resolved sections through the average Vp, Vs and Vp/Vs models allowed us to image several relevant features of the structure of the subducting Ionian slab and of the Southern Tyrrhenian mantle: -the hypocenters are localized in the NW dipping fast area (Vp>8.2 km/s), 50-60 km thick, most likely composed of: (i) 10 km of eclogite, the former oceanic Ionian crust, and (ii) 50 km of anhydrous harzburgite, i.e. the Ionian litospheric mantle. Just below, an aseismic low Vp zone (6.6 - 7.7 km/s) 20-25 km thick, is assigned to the partially hydrated (serpentinized) harzburgite. The relation between the decrease of Vp with increasing serpentinization in peridotites [Christensen, 2004] suggests that a Vp of 7.0 km/s can be achieved with a 30-40 vol % of serpentinization. The serpentinized harzburgite, which should coincide with the inner (i.e. colder) portion of the suducting slab, disappears at a depth of 230-250 km, closely corresponding to the experimentally determined maximum pressure stability of antigorite-chlorite assemblages in hydrous peridotites [ca. 8.0 GPa, Schmidt and Poli, 1998; Fumagalli and Poli, 2005]. The vanishing of the low-velocity region with increasing depth could thus be ascribed to the dehydration of the peridotite-serpentinite to less hydrous high pressure phases (e.g. the phase A) , whose seismic characteristics are akin to anhydrous lherzolite [Hacker et al., 2003]. Some other interesting features imaged in the tomography are instead related to the roots of the volcanism of the area (Aeolian islands): two vertically elongated low-velocity areas (Vp ≤ 7.0 km/s) and high Vp/Vs ratios (>1.85) characterize the mantle domains beneath Stromboli and Marsili volcanoes, reaching a maximum depth of 180 km. We relate these low-Vp, Vs and high Vp/Vs bodies to accumulation of significant amounts of mantle partial melts.
S34B-08
SAsia3D: A New Crustal and Upper Mantle P- and S- Velocity Model in Central and Southern Asia from Joint Body- and Surface-Wave Inversion
Accurate travel-time and amplitude predictions for regional seismic phases are essential for locating and characterizing small seismic events with the accuracy needed for explosion monitoring decisions. Parameter estimates calculated through 3D Earth models have the best chance of achieving acceptable prediction errors, if the models are constrained by sufficient data. With this motivation, we have developed and applied a joint body- wave/surface-wave inversion method to produce a new 3D P and S velocity model (SAsia3D) for the crust and upper mantle to a depth of 400 km in the region of central and southern Asia between 10-50° N and 40- 110° E. The method uses Pn and Pg arrival times to determine the P velocity structure and Rayleigh-wave group velocities in the period range 10-150 s to constrain the S velocity structure and depth to Moho. The body- wave and surface-wave inverse problems are coupled through an assumed correlation coefficient between P and S velocity perturbations and the imposition of bounds on the velocities and Poisson's ratio as a function of depth. Both body-wave and surface-wave forward modeling are performed in 3D models with the aid of finite-difference numerical raytracing to calculate body-wave raypaths and 2D raytracing to calculate non-great circle surface-wave paths. Nonlinearity is addressed by iterating the inversion method with updated raypaths. The regional P-wave arrival-time observations used to obtain SAsia3D were collected from the Engdahl, van der Hilst and Buland (1998; EHB) bulletin, restricted to well-located earthquakes in the years 1988-2004. The group- velocity measurements were provided by groups at the University of Colorado and Lawrence Livermore National Laboratory. Our initial model for the inversion procedure was taken as a hybrid of the CRUST2.0 3D model (Bassin et al., 2000) and the upper mantle portion of the global 1D AK135 model (Kennett et al., 1995). SAsia3D was obtained with four iterations of our technique, achieving a fifty percent variance reduction for both the body- wave and surface-wave data. Relative velocity (in particular the P velocity) variations with respect to the AK135 mantle model are in good agreement with previous studies and reflect the major tectonic features across southern and central Asia. For example, at 250 km depth, P variations across new model range from -2.0 to +2.3%, while S variations vary from -2.1% to +1.5%. We have also noted intriguing differences between the P and S velocity models in regions of significant tectonic activity, such as the Tibetan Plateau and South Caspian Basin. Some of these variations may prove useful in explaining the tectonic evolution of the Indo-Asian collision zone. We have also completed a number of validation exercises to demonstrate the accuracy of SAsia3D in regional seismic event location. Most notably, SAsia3D performs well when both regional P and S phase arrivals are included in the location. The regional P/S location obtained with SAsia3D are frequently superior to the locations obtained with a large set of teleseismic and regional P arrivals and the AK135 reference model.