S31B-1039 0800h
Towards Computing Full 3D Seismic Sensitivity: The Axisymmetric Spectral Element Method
Finite frequency tomography has recently provided detailed images of the Earth's deep interior. However, the Fr\'{e}chet sensitivity kernels used in these inversions are calculated using ray theory and can therefore not account for D''-diffracted phases or any caustics in the wavefield, as e.g. occurring in phases used to map boundary layer topography. Our objective is to compile an extensive set of full sensitivity kernels based on seismic forward modeling to allow for inversion of any seismic phase. The sensitivity of the wavefield due to a scatterer off the theoretical ray path is generally determined by the convolution of the source-to-scatterer response with, using reciprocity, the receiver-to-scatterer response. Thus, exact kernels require the knowledge of the Green's function for the full moment tensor (i.e., source) and body forces (i.e., receiver components) throughout the model space and time. We develop an axisymmetric spectral element method for elastodynamics to serve this purpose. The axisymmetric approach takes advantage of the fact that kernels are computed upon a spherically symmetric Earth model. In this reduced dimension formulation, all moment tensor elements and single forces can be included and collectively unfold in six different 2D problems to be solved separately. The efficient simulations on a 2D mesh then allow for currently unattainable high resolution at low hardware requirements. The displacement field $\mathbf{u}$ for the 3D sphere can be expressed as $\mathbf{u}(\mathbf{x},\mathrm{t})$=$\mathbf{u}(\mathbf{x}_{\phi=0},\mathrm{t})\mathrm{f}(\phi)$, where $\phi$=$0$ represents the 2D computational domain and $\mathrm{f}(\phi)$ are trigonometric functions. Here, we describe the variational formalism for the full multipole source system and validate its implementation against normal mode solutions for the solid sphere. The global mesh includes several conforming coarsening levels to minimize grid spacing variations. In an effort of algorithmic optimization, the discretization is acquired on the basis of matrix-matrix operations and performed via unit-stride cache access.
S31B-1040 0800h
Finite-Frequency Sensitivity Kernels for Boundary Topography Perturbations
Seismic traveltimes are sensitive not only to three-dimensional wavespeed variations, but also to the topography of internal discontinuity surfaces within the earth. This sensitivity has been exploited in a number of tomographic investigations, beginning with the largely unsuccessful attempts to infer the long-wavelength topography of the core-mantle boundary using the ISC-archived traveltimes of {\textit PcP\/} and {\textit PKP\/} waves. More recent studies have focused upon the topography of the major upper-mantle discontinuities and the thickness of the transition zone between 410 and 660 km depth, using the traveltimes of {\textit SdS\/} and {\textit PdP\/} bottomside reflections and {\textit Pds\/} conversions. All of these studies are based upon infinite-frequency ray theory, which maps a traveltime perturbation to the topographic perturbation at a single point on the source-to-receiver geometrical ray. A number of workers have questioned the applicability of ray theory to finite-frequency {\textit SdS\/} traveltimes, because the X-shaped Fresnel zones of off-ray sensitivity are particularly extensive for these minimax-time waves. In this paper, we use the Born approximation in conjunction with body-wave ray theory to derive a finite-frequency, boundary-topography Fr\'{e}chet kernel that fully accounts for the off-ray sensitivity of a transmitted, reflected or converted wave.
S31B-1041 0800h
Tomography, Adjoint Methods, Time-Reversal, and Banana-Doughnut Kernels
We demonstrate that Fr\'{e}chet derivatives for tomographic inversions may be obtained based upon just two calculations for each earthquake: one calculation for the current model and a second, `adjoint', calculation that uses time-reversed signals at the receivers as simultaneous, fictitious sources. For a given model~$m$, we consider objective functions $\chi(m)$ that minimize differences between waveforms, traveltimes, or amplitudes. We show that the Fr\'{e}chet derivatives of such objective functions may be written in the generic form $\delta\chi=\int_VK_m(\mathbf{x})\,\delta\ln m(\mathbf{x})\,d^3\mathbf{x}$, where $\delta\ln m=\delta m/m$ denotes the relative model perturbation. The volumetric kernel $K_m$ is defined throughout the model volume $V$ and is determined by time-integrated products between spatial and temporal derivatives of the regular displacement field $\mathbf{s}$ and the adjoint displacement field $\mathbf{s}^\dagger$ obtained by using time-reversed signals at the receivers as simultaneous sources. In waveform tomography the time-reversed signal consists of differences between the data and the synthetics, in traveltime tomography it is determined by synthetic velocities, and in amplitude tomography it is controlled by synthetic displacements. For each event, the construction of the kernel $K_m$ requires one forward calculation for the regular field $\mathbf{s}$ and one adjoint calculation involving the fields $\mathbf{s}$ and $\mathbf{s}^\dagger$. For multiple events the kernels are simply summed. The final summed kernel is controlled by the distribution of events and stations and thus determines image resolution. In the case of traveltime tomography, the kernels $K_m$ are weighted combinations of banana-doughnut kernels. We demonstrate also how amplitude anomalies may be inverted for lateral variations in elastic and anelastic structure. The theory is illustrated based upon 2D spectral-element simulations.
S31B-1042 0800h
Using Adjoint Methods to Construct 3-D Banana-Doughnut Kernels
We use adjoint methods popular in climate and ocean dynamics to calculate Fr\'{e}chet derivatives for tomographic inversions. The Fr\'{e}chet derivative of an objective function $\chi(m)$, where $m$ denotes the Earth model, may be written in the generic form $\delta\chi=\int K_m(\mathbf{x})\,\delta\ln m(\mathbf{x})\,d^3\mathbf{x}$, where $\delta\ln m=\delta m/m$ denotes the relative model perturbation. We construct the 3-D finite-frequency `banana-donut' kernel $K_m$ by simultaneously computing the so-called `adjoint' wave field $\mathbf{s}^\dagger$ forward in time and reconstructing the regular wave field $\mathbf{s}$ backward in time. The adjoint wave field is produced by using time-reversed signals at the receivers as fictitious, simultaneous sources, while the regular wave field is reconstructed on the fly by propagating the last frame of the wave field saved by a previous forward simulation backward in time. The approach is based upon the spectral-element method, and only two simulations are needed to produce density, shear-wave speed, and compressional-wave speed sensitivity kernels. This method is applied to a complicated 3-D southern California model. Various density, shear-wave speed, and compressional-wave speed sensitivity kernels are presented for different phases in the seismograms. These kernels illustrate the sensitivity of the observations to the structural parameters, and form the basis for fully 3-D tomographic inversions to update the structural model.
S31B-1043 0800h
Tomography of the global Earth based on the Spetral Elements Method
We present a new tomographic method based on the non linear least squares inversion of seismograms using the spectral elements method (SEM). The SEM is used for the forward modeling and to compute partial derivatives of seismograms with respect to the model parameters. The main idea of the method is to use a special data reduction scheme to overcome the prohibitive numerical cost of such a inversion. The SEM allows us to trigger several sources at the same time within one simulation with no incremental numerical cost. Doing so, the resulting synthetic seismograms are the sum of seismograms due to each individual source for a common receiver and a common origin time, with no possibility to separate them afterward. These summed synthetics are not directly comparable to data, but using the linearity of the problem with respect to the seismic sources, we can sum all data for a common station and a common zero time, and we perform the same operation on synthetics. Using this data reduction scheme, we can then model the whole dataset using a single SEM run, rather than a number of runs equal to the number of events considered, allowing this type of inversion to be feasible on a reasonable size computer. We will present different tests that show the feasibility of the method. It appears that this approach can work owing to the combination of two factors: the off-path sensitivity of the long period waveforms and the presence of multiple-scattering, which compensate for the loss of information in the summation process. We will discuss the advantages and drawbacks of such a scheme
S31B-1044 0800h
Seismic Applications of the Fast Marching Method: Traveltime Prediction in Complex Layered Media and 3-D Tomographic Imaging
The accurate prediction of seismic traveltimes is required in many areas of seismology, including the processing of seismic reflection profiles, earthquake location, and seismic tomography at a variety of scales. We present two seismic applications of a recently developed grid based numerical scheme for tracking the evolution of monotonically advancing interfaces, via finite difference solution of the eikonal equation, known as the fast marching method (FMM). Like most other grid-based eikonal schemes, FMM can only locate the first-arrival phase in continuous media; however, unlike other schemes, it combines unconditional stability with rapid computation. The first application of FMM that we consider focuses on the prediction of multiple reflection and refraction phases in complex 2-D layered media. By treating each layer that the wavefront enters as a separate computational domain, we show that sequential application of FMM can be used to track phases comprising any number of reflection and transmission branches in media of arbitrary complexity. We also show that the use of local grid refinement in the source neighbourhood, where wavefront curvature is high, significantly improves the accuracy of the scheme with little extra computational expense. The second application of FMM that we present is in the context of 3-D teleseismic tomography. We show that FMM can rapidly and robustly calculate two point traveltimes from an impinging teleseismic wavefront to a receiver array located on the surface, despite the presence of significant lateral variations in wavespeed in the intervening crust and upper mantle. Combined with a rapid subspace inversion method, the new FMM based tomographic scheme is shown to be extremely efficient and robust. For example, an iterative non-linear teleseismic tomography problem involving nearly 10,000 unknowns and 5,938 ray paths is solved in less than 10 minutes on a 1.6 GHz Opteron PC running Linux.
S31B-1045 0800h
Travel Times of Rays and Waves in Strongly Heterogeneous Media
Numerical validation experiments demonstrate that the finite-frequency travel times predicted by Born-Fr\'{e}chet kernel theory agree well with those measured by cross-correlation of synthetic seismograms ({\it Hung et al.}, 2001; {\it Baig et al.}, 2003). However, the agreement wanes with the enhanced strength of velocity heterogeneity. When the non-linear effect of high-order perturbations becomes important, the genuine path trajectory of a ray or wave that gives the minimum travel time differs from the stationary unperturbed one assumed in both linearized kernel and ray theories. In addition to the predictions from linearized ray theory and Fr\'{e}chet kernel theory, we extend the study of {\it Baig et al.} (2003) by exploring the ground-truth travel times in 3-D acoustic media of large velocity perturbations ($\varepsilon\ge$4%) predicted by generic ray theory based on a ray-bending method. The ground-truth time shifts are determined by cross-correlation of pairs of synthetic waveforms in one homogeneous medium and the other with Gaussian random wavespeed perturbations. The results reveal that (1) the kernel-predicted travel times are consistently better than those from ray theory in all circumstances for small perturbations ($\varepsilon\le$2%). Nevertheless, the predictions from both linearized theories are significantly deviated from the ground-truth ones for $\varepsilon\ge$4% and the heterogeneity scale $a<\sqrt{{\lambda}L}$, where $\lambda$ and $L$ are the characteristic wavelength and propagation distance, respectively. (2) Generic ray theory is superior to finite-frequency kernel theory for $\varepsilon\ge$3%, as long as $a$ is greater than $\sim{0.4}\sqrt{{\lambda}L}$. (3) While $a\le\sim{0.4}\sqrt{{\lambda}L}$, the effect of wavefront healing becomes pronounced and thus high-frequency approximation of generic ray theory begin to overestimate the travel times for all the heterogeneity strengths. In short, Fr\'{e}chet kernel theory taking finite-frequency diffractive effects into account can interpret seismic travel times very well and is tenable to resolve small-scale, weakly heterogeneous structures. Whereas generic ray theory including detoured path trajectories provides better approximations to high-frequency travel times in strongly but large-scale heterogeneous media.
S31B-1046 0800h
How Well Do Tomographic Images Explain Surface Waveforms?
Surface wave tomography based on phase measurements has achieved great successes in the last decades. However, in order to interpret the fine details of tomographic models, it is necessary to go beyond classical resolution checks, using tests with data not employed in the construction of the images. Moreover, it is desirable to apply different forward modelling techniques to assess inaccuracies due to simplified theoretical formalisms. Seismic wave amplitudes potentially provide complementary constraints on the Earth's structure and should be used to check the accuracy of existing models. We have implemented surface wave full ray theory (FRT) to investigate observed minor and major arc three-component surface waveforms, giving special emphasis to amplitudes. Calculations are performed upon mantle model S20RTS and for three sets of published phase velocity maps, combined with global crust model CRUST2.0. The synthetic waveforms calculated fit the data better than the great circle approximation (GCA). To understand the reasons for this improvement we study amplitude and phase anomalies and how they are controlled by the source, path and receiver terms. We have conducted a synthetic numerical study of phase and amplitude anomalies due to each of the terms, as predicted by the different phase velocity maps. The resulting phase perturbation is mostly controlled by path effects. The influence of the local structure at the receiver is generally negligible compared with the other terms, except for the amplitudes of the horizontal components of Rayleigh waves. However, for Rayleigh waves, theoretical amplitude anomalies due to the local structure at the source and the deviation in takeoff azimuth of the ray leaving the source are more important than previously thought. The resulting perturbation is as important as that produced by focusing and defocusing effects along the propagation path. We subsequently present results of a statistical experiment to assess the impact of these effects on observed waveforms. The improved waveform fit of full ray theory synthetics is mainly caused by more accurate phase predictions due to path corrections. The predicted amplitude anomalies correlate with observations. Nevertheless, the observed amplitude signal is not fully explained, so that on average, full ray theory does not improve the amplitude fit compared with that for a spherical Earth model. Finally, we investigate the influence of uncertainties in the source mechanism and location on the waves, by implementing an iterative inversion procedure similar to the Harvard-CMT method to determine earthquake source parameters. We use the forward modelling techniques FRT and GCA upon the laterally varying models. In all cases using the new source locations reduce the misfit, but the considered laterally varying models still do not improve the average amplitude match compared with a spherical Earth calculation. The phase velocity maps considered in this study cannot explain observed surface wave amplitudes alone, even after correcting for the earthquakes source location and mechanism. This probably reflects unmodelled effects, such as lateral variations in attenuation, azimuthal anisotropy or, possibly, stations miscalibrations.
S31B-1047 0800h
Envelope Synthesis In Random Media - Radiative Transfer Versus Finite Difference Modeling
The analysis of the coda portion of seismograms is an effective strategy to investigate the heterogeneous structure of the Earth at small scales. Usually the shape of seismogram envelopes at high frequencies are studied. A powerful method to synthesize envelopes is based on the radiative transfer theory, which describes energy transport through a scattering medium. The radiative transfer equations can conveniently be solved by a Monte Carlo simulation of random walks of energy particles through such a medium. Between single scattering events each particle moves through the background medium along ray paths. The probability of a scattering event is determined by the mean free path length depending on the total scattering coefficient of the medium. Monte Carlo simulations have so far mostly assumed isotropic scattering and acoustic approximations, as well as isotropic source radiation. Here we present an extension of this method to the full elastic case including P and S waves, and for angular dependent scattering coefficients according to the Born approximation. In order to validate this procedure, the results of the simulations are compared to envelopes obtained from full wave field modeling in 2D employing a finite difference method. Envelope shapes agree remarkably well for both short and long lapse times and for a broad range of scattering parameters. This leads to the conclusion that the use of Born scattering coefficients does not pose severe limits to the validity range of Monte Carlo method. From the comparison between elastic and acoustic simulations it becomes apparent that wave type conversions should not be neglected, especially when both P and S coda are interpreted simultaneously. Additionally, the influence of density fluctuations on envelope shapes has also been studied. It appears that the amount of density variations has a large effect on the level of the late coda only, thus showing a possibility to discriminate between velocity and density fluctuations.
S31B-1048 0800h
Polarization of Rayleigh Waves as a Function of Depth for Different Poisson's Ratio
\par In this work we study the polarization of Rayleigh waves in relation to depth for different Poisson's ratio. Rayleigh waves have the following features: the amplitude decreases rapidly with depth and the velocity is smaller than the velocity of body waves (P and S). Rayleigh waves are of great importance both in earthquake seismology and exploration seismology. Perhaps the best example is the presence of a Rayleigh wave called ground roll in the data collected in the field. In general, ground roll is considered noise in the seismic data, so that must be removed during the seismic processing stage, or at least have its effect attenuated. \par In the general theory of Rayleigh waves appears a cubic equation with a very special role. The roots of this cubic equation can be related to the Poisson's ratio $\nu$. In general one assumes the condition of a Poisson solid, where $\nu$ = 0.25. Few researchers payed attention to other values of $\nu$. In this work we calculated the velocity values of Rayleigh waves within a wide range of Poisson's ratio, from 0.00 to 0.50, where the interval is $\Delta \nu$ = 0.01. The choice of this range is due to the fact that zero is the minimum theoretical value for Poisson's ratio and 0.50 is the maximum one. \par We show the relation between $\nu$ and three quantities: $c/c_2$, $c/c_1$, and $c_2/c_1$, where $c_1$ is the velocity of the P wave, $c_2$ is velocity of the S wave, and $c$ is the velocity of Rayleigh wave. The first quantity may be denoted by $\gamma = c/c_2$. For each Poisson's ratio there are three roots, since the equation is cubic. For some values of Poisson's ratio although the velocity of the Rayleigh wave exist, that is, it is real, the surface wave itself does not exist. \par For each Poisson's ratio we calculated the path of particles confirming the elliptic behavior. The particle motion is retrograde when the displacements $u$ and $v$ are positive, and direct when $u$ is negative and $v$ is positive. \par In the case of Poisson solid ($\nu$ = 0.25), the real root of the cubic equation is $\gamma^2$ = 0.8453, or $\gamma$ = 0.9194, which means that the surface wave velocity $c$ is 0.9194 times the value of S wave velocity $c_2$. It is this value generally mentioned in textbooks.
S31B-1049 0800h
Accuracy & Computational Considerations for Wide--Angle One--way Seismic Propagators and Multiple Scattering by Invariant Embedding
Pseudodifferential operators (PSDOs) yield in principle exact one--way seismic wave equations, which are attractive both conceptually and for their promise of computational efficiency. The one--way operators can be extended to include multiple--scattering effects, again in principle exactly. In practice approximations must be made and, as an example, the variable--wavespeed Helmholtz equation for scalar waves in two space dimensions is here factorized to give the one--way wave equation. This simple case permits clear identification of a sequence of physically reasonable approximations to be used when the mathematically exact PSDO one--way equation is implemented on a computer. As intuition suggests, these approximations hinge on the medium gradients in the direction transverse to the main propagation direction. A key point is that narrow--angle approximations are to be avoided in the interests of accuracy. Another key consideration stems from the fact that the so--called ``standard--ordering'' PSDO indicates how lateral interpolation of the velocity structure can significantly reduce computational costs associated with the Fourier or plane--wave synthesis lying at the heart of the calculations. The decision on whether a slow or a fast Fourier transform code should be used rests upon how many lateral model parameters are truly distinct. A third important point is that the PSDO theory shows what approximations are necessary in order to generate an exponential one--way propagator for the laterally varying case, representing the intuitive extension of classical integral--transform solutions for a laterally homogeneous medium. This exponential propagator suggests the use of larger discrete step sizes, and it can also be used to approach phase--screen like approximations (though the latter are not the main interest here). Numerical comparisons with finite--difference solutions will be presented in order to assess the approximations being made and to gain an understanding of computation time differences. The ideas described extend to the three--dimensional, generally anisotropic case and to multiple scattering by invariant embedding.
S31B-1050 0800h
Filtering Applications on Seismogram Analysis
Useful data extraction from registered signals is an everyday problem. Seismogram filtering can be used to enhance the earthquake signal, especially, knowing that period range is strongly amplified by the seismograph itself. The goal of this project was to estimate tradeoffs between noise reduction and the sharpness of the filtered earthquake signal. The different filter bands were applied on the seismograms from different epicentral distances, including both local and regional earthquakes, in order to see the quality of the data. During the seismogram analysis we had two approaches: reading the wave phase (arrival time, amplitude, period) and spectral analysis, since frequency filtering is essential. The well-known set of Fourier orthogonal basis functions was used for spectral analysis. The main purpose of filtering was to remove the noise from the seismograms and to get as much possible information from the record. Practical application was obtained using custom made "Analysis" software. We processed 5 earthquake signals; four local and one distant earthquake. The data came from Djerdap Seismological Station, Serbia. Filtering was based on a band-pass filter with frequency cut-offs from 0.2Hz to 10Hz. Every record was filtered with five filters (10-1Hz,10-2.5Hz, 2.5-1Hz, 1-0.2Hz). The manipulation with different filter bands allowed us to define precisely the beginning of the earthquakes and first arrivals of different wave phases. The best result for the 45km epicentral distance earthquake was obtained using window 0.1-0.5sec, where very clear Pg and Sg arrivals could be distinguished. The distant earthquake showed the advantages of filtering seismograms to enhance the long period energy. The lowest frequency band-pass filter (1-0.2Hz) was the most successful. The results of the experiment were with very clear Pg and Sg arrivals when the right bands were used for filtering. This improved the quality of the original signals by emphasizing the temporal and frequency characteristics of interest.
S31B-1051 0800h
Time-Frequency analysis of seismic signals using the Wigner-Ville Spectrum
Fourier analysis and the spectral characterization of time series have played an important role in modern seismological studies. The frequency contents of a stationary time series are generally expressed by the use of the power spectrum. In seismology, we deal with earthquake records that are non-stationary transient time sequences. The usual spectral analysis fails to provide an acceptable estimate of the frequency contents of the signal with time because temporal information is lost and thus it is unsuitable for earthquake records. Different approaches have been proposed in the last 25 years to obtain an appropriate time-frequency representation, with the restriction imposed by the uncertainty principle, precluding one to obtain arbitrarily fine resolution in time and frequency. The Wigner-Ville distribution (or its stochastic counterpart the Wigner-Ville spectrum) has many of the desired properties for a power distribution: reduces to the PSD for stationary signals, its real-valued everywhere and is capable of correctly reflecting input-output relationship for linear filters. It has also some drawbacks, mainly the positivity of all its values cannot be guaranteed, which would be desirable property in terms of a density function. We apply the Wigner-Ville spectrum estimation to earthquakes ($M = 1-5$) recorded by the ANZA seismic network in Southern California and discuss the results. A measure of the degree of non-stationarity of the stochastic process is estimated, to test if the assumption of locally stationarity of the earthquake record is acceptable.
S31B-1052 0800h
CIP-MOC Modeling of Seismic Wave Propagation in Elastic Media
In many fields such as hydrodynamics and MHD, the CIP method, an upwind difference hyperbolic equation solver, has widely been employed for advection calculation. The CIP scheme was constructed considering that an advected property and its spatial derivative follow same advection equation. This effects low numerical dispersion and relaxed CFL condition in the advection calculation. In the present work, we developed a CIP-MOC (CIP with method of characteristics) scheme for seismic wave propagation in 3D elastic heterogeneous media with flat free surface. 3D elastic wave equations in velocity-stress formulation and their spatial derivatives, as well, are converted into sets of 1D advection equations and non-advection equations for each direction (x,y,z in Cartesian coodinate system) with the method of characteristics. Since the Riemann invariant of each advection equation consists of stress and velocity, updatings of velocity and stress are simultaneous and a collocated grid system is employed. A free surface is modeled as a zero-stress surface. A reflection free boundary is installed by considering no incident wave comes from outside of the boundary. A double coupled seismic point source is introduced as external point stresses. Overall scheme is made up of multiphases employing time-splitting and directional-splitting techniques. Each time step is composed of three directional updating phases each for wave propagation in x, y and z direction. Each directional updating phase is made up of advection phase and non-advection phase. In the advection phase, advection equations are solved with the CIP method. In the non-advection phases, non-advetion equations and boundary conditions are evaluated with central finite differences. We conducted CIP-MOC seismic wave propagation simulations in a half-space, layered and fully heterogeneous media for embedded point source. By comparing our products with those produced with discrete wavenumber method and finite difference method, we found that the newly developed scheme successfully calculated heterogeneous seismic wave propagations.
S31B-1053 0800h
Estimating Scattering Strength from the Transition to Equipartitioning
A pulse launched in a strongly heterogeneous medium first travels as a coherent packet of energy. Suppose the packet travels from one receiver to another. Then cross-correlating the signal at the two receivers will show a (causal) correlation peaked at the time associated with the group velocity of the packet. As time progresses (i.e., after a few mean free times) scattering causes the pulse to spread out and its energy to propagate diffusively. This means that there will be some energy propagating from the second receiver back to the first, resulting in a build-up of acausal correlations. Eventually, after many mean free times, there is as much energy propagating from one receiver to the next as there is back. In this "equipartitioning" regime, after suitable spatial averaging, the cross correlation of the diffuse wave speckle, leads to a completely symmetric result, corresponding the the sum of the causal and the anti-causal Green Function of the average medium. We will show results from a laboratory experiment tracking this transition in a granular rock sample. The amount of time required for this transition depends on the scattering properties of the medium. Through the symmetry of the two peaks it is possible to extract the mean-free path, which is a measure of the scattering strength of the sample.
S31B-1054 0800h
Modelling Borehole Wave Propagation using a 3D Variable-Grid Finite-Difference Scheme
Generating accurate synthetic seismograms for wave propagation in a three-dimensional fluid-filled borehole is difficult using conventional uniform-grid finite-difference schemes. The low-velocity modes generated in the borehole require very fine sampling to avoid numerical dispersion. The need for accurate discretisation of the borehole and casing also calls for fine grid spacing. In uniform-grid finite-difference schemes this fine grid spacing must be used throughout the model, which restricts such models to small dimensions, unrealistic velocity structures, or long wavelengths due to computer memory limitations. Using a variable-grid finite-difference scheme improves computational efficiency by partially avoiding oversampling of the wavefield in high velocity materials. The grid spacing is tailored to the velocity model, allowing the use of fine grid spacing inside the borehole and coarser grid spacing in the rock formation. We have implemented the variable-grid scheme of Pitarka (1999) in a parallel, 3D anisotropic finite-difference code, allowing large models with complex velocity structure to be tackled. First, the accuracy of this scheme for simple plane-layered geometries is verified by comparison with results obtained using the pseudo spectral method. The suitability of the scheme for acoustic logging problems is then assessed: finite-difference results for monopole and dipole sources in an open borehole are compared with analytical results obtained using the real axis integration method. Finally, the method is applied to an investigation of the modes generated in a borehole penetrating a fractured zone. Data from single well seismic surveys has suggested that this method holds promise for locating gas-filled fracture zones (Majer et al, 1997). An anomalous event observed in single-well data from the Newberry well is thought to indicate the presence of a fractured zone at some distance from the well. We use the variable-grid finite-difference scheme to investigate whether a fractured zone may act as a type of low-velocity wave guide which traps energy, giving rise to the observed high amplitude low-velocity arrivals.
S31B-1055 0800h
VSP and Crosswell Seismic Monitoring of a $CO_2$ Injection
The U.S. Department of Energy is conducting a small scale $CO_2$ injection to test issues related to the geologic sequestration of $CO_2$. The injection is taking place in southeast Texas in the Frio Formation. As part of this injection test, Lawrence Berkeley National Laboratory is conducting time-lapse seismic monitoring. Both VSP and crosswell surveys have been acquired as pre-injection baselines with a repeat post-injection survey planned. The VSP used an 80 level, 3-component geophone string with 7.5 m sensor spacing and explosive sources. Eight source shot points were acquired. The sensor string locations were interleaved resulting in 1.5 to 7.5 m spacings. The shotpoints were offset 100 to 1500 m from the sensor well on 4 azimuths. The location of the shotpoints was designed to monitor the estimated $CO_2$ plume location and to provide structural information at the injection site. VSP data reveal good quality direct and reflected events. Upgoing and downgoing wavefields will be shown. Data repeatability is investigated because up to 10 shot holes were used at each shot point. Initial results show good repeatability for separate explosive shots. The crosswell survey also used the 80 level, 3-component geophone string along with an orbital vibrator borehole source. The orbital vibrator is capable of generating both P- and S-wave direct arrivals with frequency content of about 70 - 350 Hz. For the crosswell survey, both source and sensor spacing was 1.5 m. The crosswell survey was conducted using the planned injection well (for sensors) and a nearby monitoring well (for source) which is about 30 m offset. Crosswell source locations spanned about 75 m, centered on the injection interval. The crosswell sensors spanned about 300 m., also centered on the injection interval. The planned injection interval is the 6-7 m thick, upper C sand in the Frio formation which is at a depth of about 1500 m. Initial analysis of the crosswell data shows good quality P- and S-wave direct arrivals. Tomographic imaging is planned for both P- and S-waves. Current results, probably including data from post-injection surveys, will be presented.
S31B-1056 0800h
Numerical Demonstration of A New Method of Computing Wave Fields as A Realization of An Isolated Linear Dynamic System Excited Periodically
A strong demand exists for the computation code of wave field in an arbitrary medium excited by a sinusoidal force for the system designing in ACROSS technology and also for the data analysis. Numerical examination is made on the new method of wave field computation proposed by Kumazawa et al., (2003). The theory is expected to be applied commonly to both elastic and electromagnetic waves in any heterogeneous medium with anisotropy and frequency dependent property. Wave field $w(t,x)$ is described by an inhomogeneous differential equation of the second order with respect to time and space, $DW = e$, where $D$ ($= p_1(\partial /\partial t)^2+\nabla (p_2 \nabla)$ ) is a dynamic operator ($p_1$ & $p_2$: two material parameters, $\nabla$: nabla) and $e$ is an excitation as an inhomogeneous term. Whereas the spatial discontinuity in $p_1$ & $p_2$ does not allow the space-differentiation in ordinary sense, we do differentiate all parameters and variables in wave equation by introducing hyper-function in their description in space. Then we can rewrite the wave equation in a local form of $D(\omega , \kappa, x)w(x)=e(x)$, where both $\omega$ and $\kappa$ originate from time and space derivatives, respectively. The consequence is the derivation of a set of two linear dynamic system equations of global nature; $D(\omega , \kappa',\kappa)w(\omega ,\kappa)$ = $e(\omega , \kappa')$ and $w(\omega , \kappa ) = R(\omega , \kappa ,\kappa') e(\omega , \kappa')$, which are mutually inverse with $DR=1$ at a specified frequency $\omega$, where $\kappa$ and $\kappa'$ are wavenumber vectors of wave field and excitation, respectively. The $R(\omega , \kappa ,\kappa')$ is the `frequency-wavenumber response characteristics', relating the wavenumber spectrum of excitation $e(\omega , \kappa')$ to that of the wave field $w(\omega , \kappa)$ linearly. Validity of this theory is numerically demonstrated with a simple model in 1D medium. There are two paths to derive $R(\omega , \kappa ,\kappa')$ numerically from $D(\omega , \kappa , x)$ : Path-1 via $D(\omega , \kappa',\kappa)$ involving large matrix inverse leads to more accurate result, whereas path-2 via $1/D(\omega , \kappa , x)$ without large computational load gives only approximate solution to data.
S31B-1057 0800h
Rotational Motions from Teleseismic Events - Modelling and Observations
Currently only ring lasers technology is capable of recording rotational motions resulting from earthquakes with a sensitivity and frequency band that are interesting for broadband seismology. One of those instruments is located at the Geodetic observatory in Wettzell/Germany. Here we present theoretical studies of rotational motions simulated with different Earth models and comparisons with several observations at the Wettzell ring laser. The 3-D global simulations were performed with the Spectral Element Method (Komatitsch and Tromp 2002a,b), that was modified to also allow the output of rotational seismograms. The Earth models used in these simulations range from simple radially symmetric ones, such as PREM, to more complex models including 3D velocity structures, attenuation and geometric effects like topography and bathymetry. Thus, by comparison of the theoretical rotation rates with the ring laser data we show how the results converge to the observed rotation rates when using more realistic Earth models. In a second step we compare rotation rates to the transverse component of translational acceleration both obtained from simulations with 3D velocity structures in crust and mantle. As expected from theory - under the assumption of plane wave propagation - those two signals should be in phase and scale linearly with the phase velocity. Using this relation, it is possible to determine the local phase velocity of transverse signals from collocated measurments of rotations and transverse accelerations. We compare the estimated phase velocities with those observed in a temporary seismic array installed around the ring laser.
S31B-1058 0800h
Comparison of ring laser rotational measurements with rotations derived from seismic array data: the M6.4 Morocco earthquake of February 24, 2004
Recently, ring laser technology has provided the first consistent observations of rotational motions induced by distant large earthquakes. Consistent in this context implies that the observed waveforms and amplitudes are in large part compatible with the theoretical prediction that - assuming plane wave propagation - transverse acceleration and rotation rate should be in phase and their ratio identical to twice the horizontal phase velocity. The ring laser installed at the Fundamentalstation Wettzell, South East of Germany, is recording the vertical component of rotation rate. Theoretically, this vertical component of rotation rate is a linear combination of the space derivatives of the horizontal components of translations. This suggests that - at least in principle - rotation could be determined from seismic array measurements using "finite differencing". For this purpose, we installed a cross-shaped array of nine (LE3D/5s-Mars Lite) seismometers around the ring laser instrument in Wettzell. The experiment was running for four months, from December 2003 to March 2004, recording several large earthquakes. The data set is complemented by the GRSN broadband instrument installed at Wettzell. Space derivatives are calculated by applying three different methods: A derivative weighted constant based on polynomial inversion, derivative calculation using triangular points and geodetic inversion of the surface strain and stress tensor as function of time from observed displacements. These procedures were applied to the M6.4 Morocco earthquake of February 24, 2004. We note that this is the first experiment where - in addition to seismic array data - direct measurements of a collocated rotational sensor are available.
S31B-1059 0800h
Poroelastic Wave Propagation With a 3D Velocity-Stress-Pressure Finite-Difference Algorithm
Seismic wave propagation within a three-dimensional, heterogeneous, isotropic poroelastic medium is numerically simulated with an explicit, time-domain, finite-difference algorithm. A system of thirteen, coupled, first-order, partial differential equations is solved for the particle velocity vector components, the stress tensor components, and the pressure associated with solid and fluid constituents of the two-phase continuum. These thirteen dependent variables are stored on staggered temporal and spatial grids, analogous to the scheme utilized for solution of the conventional velocity-stress system of isotropic elastodynamics. Centered finite-difference operators possess 2nd-order accuracy in time and 4th-order accuracy in space. Seismological utility is enhanced by an optional stress-free boundary condition applied on a horizontal plane representing the earth's surface. Absorbing boundary conditions are imposed on the flanks of the 3D spatial grid via a simple wavefield amplitude taper approach. A massively parallel computational implementation, utilizing the spatial domain decomposition strategy, allows investigation of large-scale earth models and/or broadband wave propagation within reasonable execution times. Initial algorithm testing indicates that a point force density and/or moment density source activated within a poroelastic medium generates diverging fast and slow P waves (and possibly an S-wave)in accord with Biot theory. Solid and fluid particle velocities are in-phase for the fast P-wave, whereas they are out-of-phase for the slow P-wave. Conversions between all wave types occur during reflection and transmission at interfaces. Thus, although the slow P-wave is regarded as difficult to detect experimentally, its presence is strongly manifest within the complex of waves generated at a lithologic or fluid boundary. Very fine spatial and temporal gridding are required for high-fidelity representation of the slow P-wave, without inducing excessive numerical dispersion. However, this wave attenuates extremely rapidly when appreciable fluid viscosity (even for water) is assigned. Sandia National Laboratories is a multiprogram science and engineering facility operated by Sandia Corporation, a Lockheed-Martin company, for the United States Department of Energy under contract DE-AC04-94AL85000.
S31B-1060 0800h
Reflection-Image-Spectroscopy (RIS) of the South American subduction zone
A new method for extracting additional information from deep seismic reflection images will be presented in this paper. In complicated evironments it is not always possible to deliver one uniquely valid seismic image. Instead, the image will differ significantly when focusing over different frequency ranges. For example enhanced concentrations of scatterers may be masked in one frequency band and become visible in another frequency band. We therefore apply 3D prestack Kirchhoff depth migration to data containing the full frequency range on one hand and different narrow frequency ranges on the other hand. We call this approach "Reflection Image Spectroscopy" (RIS). It provides frequency dependent images which enable for instance the characterization of the medium in terms of scatterer concentration and length scale. Furthermore, the analysis of these images allows to differentiate between small-scale structures in the high-frequency band and large-scale structures in the low-frequency band. This technique was applied to an onshore deep seismic reflection line (ANCORP96) across the Chilean subduction one. The dominant features in the full-frequency image are the mid-crustal Quebrada Blanca Bright Spot (QBBS) and the subducted oceanic crust (Nazca reflector). The latter appears as a thick reflector with no further interpretable internal structure. However, the narrow-frequency images yield more information (see figure below). At a depth of 70 km the Nazca reflector changes its appearance. The clearly separated double reflection zone continues downward but a wedge-shaped body attached to its top can be resolved. We interpret this body as a possible fluid trap above the slab. Serpentinization within this part of the mantle wedge may lead to the increased reflectivity. Further up towards the QBBS a strongly heterogeneous zone can be observed. This zone coincides with a tomographic low-velocity anomaly and can be explained as the location of possible fluid ascend paths, eventually along pre-existing faults. The QBBS therefore seems to be directly linked to the downgoing plate. Another seismic reflection line (PRECORP95) in that area shows similar features and supports our interpretations.
S31B-1061 0800h
One-way Wavefield Extrapolation in Riemannian Coordinates
Images of crustal and lithospheric-scale structure can be created through seismic wavefield extrapolation, which extends surface-recorded data to depth through application of a wave-equation operator. The one-way extrapolation operators are derived from the acoustic wave-equation dispersion relation, and are usually defined in a Cartesian coordinate system with a vertical extrapolation axis. However, in this formulation there is an assumption that a vertical extrapolation axis can handle all wavefield propagation effects. In particular, there are numerous seismic imaging experiments where challenges such as incorporating rough topography, modeling overturning waves, and accounting for spherical Earth effects make extrapolation in Cartesian coordinates impractical. Riemannian wavefield extrapolation (RWE) (Sava and Fomel, 2004) is a generalization of the downward continuation concept to coordinate systems that incorporate geometrical or wave propagation effects. For example, geometrical effects such as topography can be directly incorporated into the coordinate system, and wavefield extrapolation commenced directly from relief (Shragge and Sava, AGU 2004). Another example is modeling path propagation effects directly in the coordinate system enabling accurate one-way wavefield extrapolation, even in situations where wavefields overturn. The method is implemented with a mixed domain solution involving both Fourier and space-domain finite differences. The algorithm is modestly more costly than wavefield extrapolation performed on a Cartesian computational mesh.
S31B-1062 0800h
Reservoir Imaging Using Frequency-Dependent Seismic Attributes
A recently developed elastic fluid-flow model accounting for the physical phenomena associated with frequency-dependent anomalous reflectivity at low seismic frequencies is discussed. Besides theoretical calculations, observations from field and laboratory data indicate that the reflectivity is strongly related to the reservoir flow properties. An asymptotic scaling of the frequency-dependent component of the reflection coefficient with respect to the frequency of the signal and the reservoir fluid mobility has been obtained. This scaling incorporates the reservoir fluid mobility. The porous medium tortuosity factor from Biot's equations has been linked to the relaxation time in dynamic Darcy's law.
S31B-1063 0800h
Instantaneous Polarization Attributes in the Time-frequency Domain: Application to Wave Field Separation
In this contribution we propose a method of wave field separation from multicomponent data sets based on the continuous wavelet transform (CWT). We present different approaches for obtaining the time-frequency dependent instantaneous polarization attributes for multicomponent data sets (2C, 3C or more). Using these attributes, we show how to construct filters tailored to separate (filter) different wave types followed by an inverse wavelet transform to obtain the desired wave type in the time domain. The proposed methods are applied on synthetic and experimental data for illustration.
S31B-1064 0800h
Incorporating topography into crustal-scale wave-equation imaging through conformal mapping
Migration of seismic land data acquired over topography presents significant challenges for seismic imaging technology. One technique commonly used to alleviate the deleterious effects of topography is to include a wavefield datuming step in the pre-migration processing flow. Locations exhibiting strong, near-surface lateral velocity contrast accompanying topography, though, may have significant wavefield triplication that leads to non-optimal datuming procedures; especially if Kirchhoff-based methods are used. Accordingly, the use of novel wave-equation approaches that naturally handle wavefield-multipathing should lead in these situations to significantly improved datuming or migration procedures. Conformal mapping is a technique used widely in applied physics and engineering fields to facilitate numerical solution of boundary value problems involving surfaces exhibiting complex geometry. The general reason for applying this procedure is to transform solution domains marked by complicated boundaries, such as topography, to domains that are more regular, such as Cartesian. The conformal mapping transforms may then be used to define wavefield-extrapolation equations appropriate for using in the new domain. In this paper, we demonstrate that the conformal mapping transform coupled with the Riemannian wavefield extrapolation (Sava, 2004) generate orthogonal meshes, herein termed topographic-coordinate systems, and governing equations appropriate for incorporating topography directly into wave-equation migration. This premise is illustrated with poststack and prestack crustal-scale sized synthetic data sets with a topographic surface taken from the Foothills region of British Columbia, Canada.
http://sepwww.stanford.edu/sep/students/jeff.html
S31B-1065 0800h
MINEMOTION3D: A new set of Programs for Predicting Ground Motion From Explosions in Complex 3D Media
Predicting ground motion from complicated mining explosions is important for mines developing new blasting programs in regions where vibrations must be kept below certain levels. Additionally, predicting ground motion from mining explosions in complex 3D media is important for moment estimation for nuclear test treaty monitoring. Both problems have been addressed under the development of a new series of numerical prediction programs called MINEMOTION3D including 1) Generalized Fourier Methods to generate Green's functions in 3D media for a moment tensor source implementation and 2) MineSeis3D, a program that simulates seismograms for delay-fired mining explosions with a linear relationship between signals from small size individual shots. To test the programs, local recordings (5 - 23 km) of three production shots at a mine in northern Minnesota were compared to synthetic waveforms in 3D media. A non-zero value of the moment tensor component M12 was considered, to introduce a horizontal spall component into the waveform synthesis when the Green's functions were generated for each model. Methods using seismic noise crosscorrelation for improved inter-element subsurface structure estimation were also evaluated. Comparison of the observed and synthetic waveforms shows promising results. The shape and arrival times of the normalized synthetic and observed waveforms are similar for most of the stations. The synthetic and observed waveform amplitude fit is best for the vertical components in the mean 3D model and worst for the transversal components. The observed effect of spall on the waveform spectra was weak in the case of fragmentation delay fired commercial explosions. Commercial applications of the code could provide data needed for designing explosions which do not exceed ground vibration requirements posed by the U.S. Department of the Interior, Office of Surface Mining.
S31B-1066 0800h
Explosion Shear Wave Generation and Scattering
We use observations of explosion-generated Lg together with three separate types of numerical models to determine how underground nuclear explosions generate shear wave phases. This question is fundamental to how Lg phases are interpreted for use in explosion yield estimation and earthquake/explosion discrimination. A simple point explosion in a uniform medium generates no shear waves, so the Lg phase is generated entirely by non-spherical components of the source and conversions through reflections and scattering. Our results indicate that the most important sources of high frequency explosion shear waves are P to S conversions at the free surface and S waves generated directly by a realistic distributed explosion source including nonlinear effects due to the free surface and gravity. In addition, Rg scattering may contribute to lower frequency Lg. Near source S is observed on both radial and tangential component records from a diverse set of explosion data. The data sets include 1) Degelen Mountain explosions recorded at distances less than 100 km and corresponding recordings at Borovoye (BOR) at 650 km; 2) recordings from Russian deep seismic sounding experiments; 3) Nevada Test Site (NTS) explosion sources including the Nonproliferation Experiment (NPE) and nuclear tests covering a range of source depths and media properties. We model the overburied NPE, and underburied and overburied Degelen explosions, using point sources and two-dimensional nonlinear finite difference calculations to quantify the source effects. We use energy conservation to determine an upper bound on Rg to Lg scattering. Results indicate that Rg to Lg scattering may be important at frequencies less than 1 Hz, and in Lg coda, but is less than Lg generated directly by the explosion at higher frequencies. We use 2D and 3D finite difference calculations, using the known topography and velocity structure at Degelen Mt. and lateral heterogeneities within the crust, to estimate the effect of surface and crustal scattering on Lg. Results indicate that topographic scattering can significantly disrupt the surface pS scattered phase, so more shear wave energy is trapped in the crust. The topography also appears to be the most significant scatterer of Rg, and has a profound effect on its dispersion.
S31B-1067 0800h
A new method of normal mode calculation
Numerical method of calculating normal modes of the Earth is classical problem in seismology and many methods have been developed, e.g. MINEOS, OBANI (Woodhouse 1989) and DISPER80 (Saito 1989). These methods are based on direct numerical integration of the governing differential equations (Runge-Kutta method). The methods are efficient for search for (real) eigenvalues of non dissipative cases but have some trouble to seek complex eigenvalues of dissipative or leaky oscillations. On the other hand, in solar seismology, the Henyey type relaxation method is used to calculate acoustic and gravity modes of the sun and stars (e.g. Unno et al. 1989). The latter method is generally more stable for eigenvalue problems of a self-gravitating gaseous body where the density and pressure vary exponentially near its surface. The efficiency of the method is, however, depends on initial guess of eigenvalues and eigenfunctions. Here we propose a new method to overcome the deficiency of the above methods. This method is formulated with the aid of the concept of the Henyey type relaxation method but does not needs initial guess of eigenfuctions and solves the problem more directly like the Runge-Hutta method. Some mumerical tests show good convergence of complex eigenvalues in our method. So the proposed method is effective to calculate solid modes, acoustic and gravity modes and ocean modes interfered with the other parts of the earth and planets.
S31B-1068 0800h
Spatial variations of excitation amplitudes of Earth's background free oscillations
The phenomenon of Earth's background free oscillations has now been confirmed firmly [{\it Kobayashi and Nishida,} 1998; {\it Nawa et al.,} 1998; {\it Suda et al.,} 1998]. Their excitation amplitudes and statistical features indicate that the excitation source is persistent disturbance at or just above the whole Earth's surface. The observed amplitudes also show clear annual variations and acoustic resonance between the solid Earth and the atmosphere, suggesting that the most likely excitation source is atmospheric disturbance although other possibilities, such as oceanic disturbance, cannot be ruled out. Fukao {\it et al.} [2002] calculated the synthetic spectra of Earth's background free oscillations, assuming that the atmospheric turbulence generates random pressure forces acting uniformly on the surface of the solid Earth. The calculated spectra were consistent with the observed ones. In order to constrain more narrowly the excitation mechanism of Earth's background free oscillations, we investigated the spatial distribution of the excited amplitudes more closely. We first calculated the cross-spectrum of background free oscillations for every pair of 55 stations from 3 to 6 mHz. We then stacked these cross-spectra for each of the surface points with intervals of 10 degrees in latitude and longitude after a correction for relative phase shift of Rayleigh wave propagation from every surface point to the two paired stations. The cross-spectral amplitudes show a clear spatial pattern dominant in harmonic degree 2. The amplitude maxima are on the Pacific ocean and its antipode, whereas the minima are on North America and its antipode. This result suggests that neither purely uniform atmospheric disturbance nor purely oceanic disturbance explains the observed amplitude pattern of Earth's background free oscillations.
S31B-1069 0800h
Computation of Normal Modes for a spherical Planet with a non-adiabatic atmosphere : application to the Earth and Venus
Thanks to technological advances over the past fifteen years the atmosphere is now a new medium for seismological investigation. Surface waves emitted after large earthquakes are known to induce, by dynamic coupling, atmospheric infrasonic and gravity waves propagating upward through the atmosphere. Those waves have been detected recently at ionospheric heights using a variety of techniques, such as Doppler or GPS ionospheric sounding. The theoretical description of the coupling processes and propagation in the atmosphere already take into account the viscosity of the atmosphere. However, for acoustic waves in atmosphere, the effects of thermal-conduction can reach 30% of the total dissipation. We describe a new theory taking into account the thermal conduction effects and which allows a more precise modelling of the propagation and provide the temperature fluctuation associated to the seismic waves at high altitudes. This theory is use to model waves and signals for different earthquakes, and the modelled data are compared to the observations and results discussed. Modelling is also performed in the case of Venus, and the perspectives of detection by ESA's Venus Express Spacecraft are specified.
S31B-1070 0800h
Tracing back ionospheric perturbations detected by GPS monitoring in the epicenter area after large earthquakes to their source mechanisms.
Dense Global Positioning System (GPS) continuous networks allow to detect and image small-scale perturbations of the ionosphere, through measurement of the Total Electron Content (TEC) along satellite-receiver rays. We present observations of ionospheric perturbations observed above the epicenter area during two large earthquakes, the M=8.2 Tokachi-Oki earthquake near Hokkaido (09/22/2003) and the M=6.5 San Simeon earthquake in Central California. In both cases significant variations of TEC have been observed, propagating up to several hundred kilometers from the epicenter, and starting approximately 15 minutes after the event. Very similar in nature, but with different amplitudes, those two observations can be interpreted as the signature in the ionosphere of an acoustic wave induced in the atmosphere by the ground displacement. Indeed, we compared those signals with results from finite-source inversions performed for those two events. Most similarities and differences between those two cases can be traced back to the seismic sources.