Layered Earth Model Research Articles

Finite-Difference Modeling and Dispersion Analysis of High-Frequency Love Waves for Near-Surface Applications

Love-wave propagation has been a topic of interest to crustal, earthquake, and engineering seismologists for many years because it is independent of Poisson’s ratio and more sensitive to shear (S)-wave velocity changes and layer thickness changes than are Rayleigh waves. It is well known that Love-wave generation requires the existence of a low S-wave velocity layer in a multilayered earth model. In order to study numerically the propagation of Love waves in a layered earth model and dispersion characteristics for near-surface applications, we simulate high-frequency (>5 Hz) Love waves by the staggered-grid finite-difference (FD) method. The air–earth boundary (the shear stress above the free surface) is treated using the stress-imaging technique. We use a two-layer model to demonstrate the accuracy of the staggered-grid modeling scheme. We also simulate four-layer models including a low-velocity layer (LVL) or a high-velocity layer (HVL) to analyze dispersive energy characteristics for near-surface applications. Results demonstrate that: (1) the staggered-grid FD code and stress-imaging technique are suitable for treating the free-surface boundary conditions for Love-wave modeling, (2) Love-wave inversion should be treated with extra care when a LVL exists because of a lack of LVL information in dispersions aggravating uncertainties in the inversion procedure, and (3) energy of high modes in a low-frequency range is very weak, so that it is difficult to estimate the cutoff frequency accurately, and “mode-crossing” occurs between the second higher and third higher modes when a HVL exists.

EMDPLER: A F77 program for modeling the EM response of dipolar sources over the non-magnetic layer earth models

A program EMDPLER coded in FORTRAN 77 is presented for modeling the electromagnetic (EM) response of a wide class of dipolar sources (namely, vertical magnetic dipole, horizontal magnetic dipole and horizontal electric dipole) over the layered earth model in induction ( f≤50 kHz) as well as higher frequency (50 kHz≤ f≤1000 kHz) regions. The program is based on a semi-analytical numerical approach for evaluating infinite integrals occurring in the expressions of EM field components. It is capable of computing the EM response with or without the displacement current factor for any arbitrary position of source and/or receiver either in air or on the surface of the model due to all the conventional EM dipole sources where other available algorithms fail due to one or other reasons. The program is well documented, easy to understand and amenable to modification for the user’s specification. The validity and accuracy of the program is demonstrated by computing the EM response of dipolar sources over homogeneous and multi-layer earth models, and comparing the computed results with those of published results in the literature. Response curves depict their characteristic variations. The computed results match very well with the published results for the induction region ( f≤50 kHz) and are extension of their characteristic variation in the higher frequency (50 kHz≤ f≤1000 kHz) region. The matching of computed results with the published results illustrates the validity of the program.

Hydrogeophysical exploration of three-dimensional salinity anomalies with the time-domain electromagnetic method (TDEM)

The time-domain electromagnetic method (TDEM) is widely used in groundwater exploration and geological mapping applications. TDEM measures subsurface electrical conductivity, which is strongly correlated with groundwater salinity. TDEM offers a cheap and non-invasive option for mapping saltwater intrusion and groundwater salinization. Traditionally, TDEM data is interpreted using one-dimensional layered-earth models of the subsurface. However, most saltwater intrusion and groundwater salinization phenomena are characterized by three-dimensional anomalies. To fully exploit the information content of TDEM data in this context, three-dimensional modeling of the TDEM response is required. We present a finite-element solution for three-dimensional forward modeling of TDEM responses from arbitrary subsurface electrical conductivity distributions. The solution is benchmarked against standard layered-earth models and previously published three-dimensional forward TDEM modeling results. Concentration outputs from a groundwater flow and salinity transport model are converted to subsurface electrical conductivity using standard petrophysical relationships. TDEM responses over the resulting subsurface electrical conductivity distribution are generated using the three-dimensional TDEM forward model. The parameters of the hydrodynamic model are constrained by matching observed and simulated TDEM responses. As an application example, a field dataset of ground-based TDEM data from an island in the Okavango Delta is presented. Evaporative salt enrichment causes a strong salinity anomaly under the island. We show that the TDEM field data cannot be interpreted in terms of standard one-dimensional layered-earth TDEM models, because of the strongly three-dimensional nature of the salinity anomaly. Three-dimensional interpretation of the field data allows for detailed and consistent mapping of this anomaly and makes better use of the information contained in the TDEM field dataset.

Fast, laterally smooth inversion of airborne time‐domain electromagnetic data

ABSTRACTWe present a fully developed, fast approximate method for 1D inversion of time‐domain electromagnetic data. The method is applied to a helicopterborne transient electromagnetic data set from the Toolibin Lake area of Western Australia using the lateral parameter correlation method to ensure lateral smoothness of the inverted models.The method is based on fast approximate forward computation of transient electromagnetic step responses and their derivatives with respect to the model parameters of a 1D model. The inversion is carried out with multi‐layer models in an iterative, constrained least‐squares inversion formulation including explicit formulation of the model regularization through a model covariance matrix. The method is 50 times faster than conventional inversion for a layered earth model and produces model sections of concatenated 1D models and contoured maps of mean conductivity in depth intervals almost indistinguishable from those of conventional inversion.To ensure lateral smoothness of the model sections and to avoid spurious artefacts in the mean conductivity maps, the inversion is integrated with the lateral parameter correlation method. In this way, well determined parameters are allowed to influence the more poorly determined parameters in the survey area.Applied to the Toolibin data set, the inversion produces model sections and conductivity maps that reveal the distribution of conductivity in the area and thereby the distribution of salinity. This information is crucial for any remediation effort aimed at alleviating the salinization problems.

Erratum to Improvements on Computation of Phase Velocities of Rayleigh Waves Based on the Generalized R/T Coefficient Method

Among many improvements on phase velocity computation of Rayleigh waves in a layered Earth model, the methods based on reflection and transmission (R/T) coefficients produce stable and accurate phase velocities for both low and high frequencies and are appropriate for viscoelastic and anisotropic media. The generalized R/T coefficient methods are not the most efficient algorithms. This paper presents a new, more efficient algorithm called the fast-generalized R/T coefficient method to calculate the phase velocity of surface waves for a layered Earth model. The improvements include (1) computation of the generalized R/T coefficients without calculation of the modified R/T coefficients, and (2) presentation of an analytic solution for the inverse of the 4×4 layer matrix E . Compared with the traditional generalized R/T coefficient method, the fast-generalized R/T coefficient method, when applied on Rayleigh waves, significantly improves the speed of computation, cutting the computational time at least by half while keeping the stability and accuracy of the traditional method.

Fast approximate 1D inversion of frequency domain electromagnetic data

ABSTRACTWe present a fast approximate method for 1D inversion of frequency domain data and apply it to frequency domain helicopter‐borne data from the Bookpurnong area of the Murray River, South Australia. The method is based on fast approximate forward computation of transient electromagnetic step responses and their derivatives with respect to the model parameters of a 1D model, with the frequency domain responses and derivatives then found through Fourier transformation of the time‐domain counterparts. The inversion is carried out with multi‐layer models in an iterative, constrained least‐squares inversion scheme including explicit formulation of the model regularization through a model covariance matrix. The method is 30 times faster than conventional full inversion for a layered earth model and produces model sections of concatenated 1D models and contoured maps of mean conductivity in elevation intervals almost indistinguishable from those of a conventional full inversion. In a theoretical forward and inverse modelling study, the fast approximate and conventional computation methods are compared demonstrating the applicability of the approximate method and its limitations.Applied to the Bookpurnong RESOLVE® FDHEM data set from South Australia, the inversion produces model sections and conductivity maps that reveal the distribution of conductivity in the area and thereby the distribution of salinity. This information is crucial for any remediation effort aimed at alleviating the salinization of the river and the degradation of floodplain vegetation and associated ecosystems.

Initial results from spatially averaged coherency, frequency-wavenumber, and horizontal to vertical spectrum ratio microtremor survey methods for site hazard study at Launceston, Tasmania

The Tamar rift valley runs through the City of Launceston, Tasmania. Damage has occurred to city buildings due to earthquake activity in Bass Strait. The presence of the ancient valley, the Tamar valley, in-filled with soft sediments that vary rapidly in thickness from 0 to 250 m over a few hundreds metres, is thought to induce a 2D resonance pattern, amplifying the surface motions over the valley and in Launceston. Spatially averaged coherency (SPAC), frequency-wavenumber (FK) and horizontal to vertical spectrum ratio (HVSR) microtremor survey methods are combined to identify and characterise site effects over the Tamar valley.Passive seismic array measurements acquired at seven selected sites were analysed with SPAC to estimate shear wave velocity (slowness) depth profiles. SPAC was then combined with HVSR to improve the resolution of these profiles in the sediments to an approximate depth of 125 m. Results show that sediments thicknesses vary significantly throughout Launceston. The top layer is composed of as much as 20 m of very soft Quaternary alluvial sediments with a velocity from 50 m/s to 125 m/s. Shear-wave velocities in the deeper Tertiary sediment fill of the Tamar valley, with thicknesses from 0 to 250 m vary from 400 m/s to 750 m/s.Results obtained using SPAC are presented at two selected sites (GUN and KPK) that agree well with dispersion curves interpreted with FK analysis. FK interpretation is, however, limited to a narrower range of frequencies than SPAC and seems to overestimate the shear wave velocity at lower frequencies. Observed HVSR are also compared with the results obtained by SPAC, assuming a layered earth model, and provide additional constraints on the shear wave slowness profiles at these sites. The combined SPAC and HVSR analysis confirms the hypothesis of a layered geology at the GUN site and indicates the presence of a 2D resonance pattern across the Tamar valley at the KPK site.

Approximation to Cutoffs of Higher Modes of Rayleigh Waves for a Layered Earth Model

A cutoff defines the long-period termination of a Rayleigh-wave higher mode and, therefore is a key characteristic of higher mode energy relationship to several material properties of the subsurface. Cutoffs have been used to estimate the shear-wave velocity of an underlying half space of a layered earth model. In this study, we describe a method that replaces the multilayer earth model with a single surface layer overlying the half-space model, accomplished by harmonic averaging of velocities and arithmetic averaging of densities. Using numerical comparisons with theoretical models validates the single-layer approximation. Accuracy of this single-layer approximation is best defined by values of the calculated error in the frequency and phase velocity estimate at a cutoff. Our proposed method is intuitively explained using ray theory. Numerical results indicate that a cutoffs frequency is controlled by the averaged elastic properties within the passing depth of Rayleigh waves and the shear-wave velocity of the underlying half space.

Fast approximate inversion of SkyTem airborne electromagnetic data

Transient electromagnetic (TEM) soundings have become one of the standard methods of environmental geophysics (Fittermann and Stewart, 1986; Buselli et al., 1990; Hoekstra and Blohm, 1990; Christensen and Sørensen, 1998; Auken et al., 2006; Lane and Pracilio, 2000). Over the past decades, airborne TEM methods have found widespread use in hydrogeophysical investigations, making it possible to cover large areas in a cost-effective way. Several helicopterborne TEM systems have been developed and they now represent the most up-to-date method of airborne hydrogeophysical investigations, collecting 100,000s of data sets. We present a fast approximate method for one-dimensional inversion of time domain electromagnetic data and apply it to SkyTEM helicopterborne data from the Toolibin Lake area of West Australia. The method is based on fast approximate forward computation of transient electromagnetic step responses and their derivatives with respect to the model parameters of a 1D model. The inversion is carried out with multi-layer models in a state-of-the-art formulation of a least-squares iterative inversion scheme including explicit formulation of the model regularization through a model covariance matrix. The method is 50 times faster than conventional inversion for a layered earth model and produces model sections of concatenated 1D models and contoured maps of mean conductivity in elevation intervals almost indistinguishable from those of conventional inversion. To ensure lateral smoothness of the model sections and to avoid spurious artifacts in the mean conductivity maps, the Lateral Parameter Correlation method is applied to the results of the individual inversion models. In this way, well determined parameters are allowed to influence the more poorly determined parameters in the survey area. Applied to the SkyTEM Toolibin data set, the inversion produces model sections and conductivity maps that reveal the distribution of conductivity in the area and thereby the distribution of salinity. This information is crucial for any remediation effort aimed at alleviating the salinity problems. The results are in accordance with previous results from other investigators.

Fast approximate inversion of FDHEM data

Despite advances in computer speed and code efficiency, the most often employed interpretation method for frequency domain helicopterborne electromagnetic (FDHEM) data is still inversion with one-dimensional (1D) models. In areas where the lateral rate of change of conductivity is small, inversion with 1D models is justified and can be sufficient. These conditions are often found in hydrogeophysical investigations in a sedimentary environment. Even for 1D inversion, the task of inverting a survey volume of 100,000s of data points can be time-consuming, and more or less sophisticated approximate methods have been developed. These are often included in software packages used for data and model display and handling of geographical information, e.g. EMflow (Macnae et al., 1998), EMax (Fullagar and Reid, 2001), etc. A fairly comprehensive comparison between different inversion and transformation techniques used for FDHEM data is published in Sattel (2005). We present a fast approximate method for one-dimensional (1D) inversion of frequency domain data and apply it to frequency domain helicopterborne data from the Bookpurnong area of the Murray River, South Australia. The method is based on fast approximate forward computation of transient electromagnetic step responses and their derivatives with respect to the model parameters of a 1D model, and the frequency domain responses and derivatives are found through Fourier transformation of the time domain counterparts using Fast Hankel Transform filters. The inversion is carried out with multi-layer models in a state-of-the-art formulation of a least-squares iterative inversion scheme including explicit formulation of the model regularisation through a model covariance matrix. The approximate and ordinary 1D inversion approaches thus share the inversion formulation, the difference lying only in the forward mapping. The method is 30 times faster than ordinary inversion for a layered earth model and produces model sections of concatenated 1D models and contoured maps of mean conductivity in elevation intervals almost indistinguishable from those of ordinary inversion.

Coseismic and postseismic deformation associated with the 2003 Chengkung, Taiwan, earthquake

SUMMARY We use GPS-derived coseismic and postseismic displacements of the 2003 Mw 6.8 Chengkung, Taiwan, earthquake to examine seismogenic behaviour of the southern Longitudinal Valley. We invert for fault slip on the Chihshang fault, a segment of the Longitudinal Valley fault, based on a simplified layered earth model and a listric-shape fault geometry. Our model infers that maximum coseismic slip of about 0.72 m occurred at a depth of 15 km. Conversely, postseismic slip of about 0.1 m over a 157 d period mainly occurred at depths less than 10 km. We find an anticorrelation between coseismic and postseismic slip distributions. However, the depth profiles of seismicity before and after the main shock are very similar and both show the peak seismicity at 15–25 km depth, consistent with the depth range of large coseismic slip. The potency distribution at depth reveals a coseismic slip deficit at shallow depths, which is compensated by postseismic and interseismic creep associated with the Lichi Melange. Aseismic slip is an important component of the slip budget on the Chihshang fault. According to the historic earthquake records, interseismic slip rate and coseismic slip on the Chihshang fault, the recurrence interval of this type of event is between 12 and 36 yr. We use a rate-dependent friction model of afterslip to infer the frictional parameter, dτss/dlnV= 0.03–0.5 MPa, at shallow depths on the Chihshang fault. Assuming effective normal stress of 100 MPa at 5 km depth, the rheological parameter, , is about 3 × 10−4–5 × 10−3, which is in good agreement with laboratory experiments.

On the Mathematical Analysis of an Elastic-gravitational Layered Earth Model for Magmatic Intrusion: The Stationary Case

In the early eighties Rundle (1980, 1981a,b, 1982) developed the techniques needed for calculations of displacements and gravity changes due to internal sources of strain in layered linear elastic-gravitational media. The approximation of the solution for the half space was obtained by using the propagator matrix technique. The Earth model considered is elastic-gravitational, composed of several homogeneous layers overlying a bottom half space. Two dislocation sources can be considered, representing magma intrusions and faults. In recent decades theoretical and computational extensions of that model have been developed by Rundle and co-workers (e.g., Fernandez and Rundle, 1994a,b; Fernandez et al., 1997, 2005a; Tiampo et al., 2004; Charco et al., 2006, 2007a,b). The source can be located at any depth in the media. In this work we prove that the perturbed equations representing the elastic-gravitational deformation problem, with the natural boundary and transmission conditions, leads to a well-posed problem even for varied domains and general data. We present constructive proof of the existence and we show the uniqueness and the continuous dependence with respect to the data of weak solutions of the coupled elastic-gravitational field equations.

Data-resolution Matrix and Model-resolution Matrix for Rayleigh-wave Inversion Using a Damped Least-squares Method

Inversion of multimode surface-wave data is of increasing interest in the near-surface geophysics community. For a given near-surface geophysical problem, it is essential to understand how well the data, calculated according to a layered-earth model, might match the observed data. A data-resolution matrix is a function of the data kernel (determined by a geophysical model and a priori information applied to the problem), not the data. A data-resolution matrix of high-frequency (≥2 Hz) Rayleigh-wave phase velocities, therefore, offers a quantitative tool for designing field surveys and predicting the match between calculated and observed data. We employed a data-resolution matrix to select data that would be well predicted and we find that there are advantages of incorporating higher modes in inversion. The resulting discussion using the data-resolution matrix provides insight into the process of inverting Rayleigh-wave phase velocities with higher-mode data to estimate S-wave velocity structure. Discussion also suggested that each near-surface geophysical target can only be resolved using Rayleigh-wave phase velocities within specific frequency ranges, and higher-mode data are normally more accurately predicted than fundamental-mode data because of restrictions on the data kernel for the inversion system. We used synthetic and real-world examples to demonstrate that selected data with the data-resolution matrix can provide better inversion results and to explain with the data-resolution matrix why incorporating higher-mode data in inversion can provide better results. We also calculated model-resolution matrices in these examples to show the potential of increasing model resolution with selected surface-wave data.

A study of near source effects in array-based (SPAC) microtremor surveys

Array-based methods for exploiting ambient seismic noise are receiving increasing attention in the literature, particularly for use in providing shear wave velocity profiles to assist with simulation of site-specific earthquake responses. Many of these microtremor methods—such as the spatial autocorrelation technique (SPAC)—operate under the fundamental assumption of plane wave propagation of surface waves. However, when there are significant microtremor sources in close proximity to an array, wave front curvature becomes significant, and the plane wave assumption becomes a tenuous one. This paper explains the use of a simplified geometry-based approach to examine the effect of source distance on SPAC spectra for hexagonal (six-station) and triangular (three-station) arrays. The results suggest that near source effects are generally minimal provided that most of the sources are located a distance equal to at least two array radii from the centre of the array, although this distance increases when attenuation is considered. Results of a field trial involving the deliberate use of near-sources in proximity to a hexagonal (six-station) array support the findings of the theoretical model. Near source effects generally result in the measured SPAC curve being offset to lower frequencies for wavenumbers less than the first secondary maximum of the SPAC curve (kr < 7.0; λ > 0.86R), resulting in underestimates of shear velocity values in the inverted layered earth model. For the field example, shear velocity was underestimated by up to approximately 17 per cent. Studies of near source effects in a number of theoretical source distributions for hexagonal arrays revealed that cyclic variation in the imaginary SPAC coefficients are indicative of near source effects, although field examples indicate this effect is less significant than the geometric modelling suggests. Modelling of linear microtremor source distributions (intended to represent traffic travelling along a road) resulted in less severe near source effects than full (360°) azimuthal source coverage—a result confirmed by field data. Model and field studies also indicated that in situations of a dominant, close microtremor source of limited azimuthal extent, use of a hexagonal (six-station) array significantly mitigates near source effects compared with a similarly aligned triangular (three-station) array.

Applicability of artificial neural networks for obtaining velocity models from synthetic seismic data

Seismic velocity analysis is a crucial part of seismic data processing and interpretation which has been practiced using different methods. In contrast to time consuming and complicated numerical methods, artificial neural networks (ANNs) are found to be of potential applicability. ANN ability to establish a relationship between an input and output space is considered to be appropriate for mapping seismic velocity corresponding to travel times picked from seismograms. Accordingly a preliminary attempt is made to evaluate the applicability of ANNs to determine velocity and dips of dipping layered earth models corresponding to travel time data. The study is based on synthetic data generated using inverse modeling approach for three earth models. The models include a three-layer structure with same dips and same directions, a three-layer model with different dips and same directions, as well as a two-layer model with different dips and directions. An ANN structure is designed in three layers, namely, input, output, and hidden ones. The training and testing process of the ANN is successfully accomplished using the synthetic data. The evaluation of the applicability of the trained ANN to unknown data sets indicates that the ANN can satisfactorily compute velocity and dips corresponding to travel times. The error intervals between the desired and calculated velocity and dips are shown to be acceptably small in all cases. The applicability of the trained ANN in extrapolating is also evaluated using a number of data outside of the range already known to ANN. The results indicate that the trained ANN acceptably approximates the velocity and dips. Furthermore, the trained ANN is also evaluated in terms of capability of handling deficiency in input data where acceptable results were also achieved in velocity and dip calculations. Generally, this study shows that velocity analysis using ANNs can promisingly tackle the challenge of retrieving an initial velocity model from the travel time hyperbolas of seismic data.

Amplitude and phase correction of helicopter EM data

We aim to develop a quantitative method for recalibration of historic helicopter electromagnetic data sets. Recent research has shown that frequency-domain helicopter electromagnetic data collected over a conductive half-space such as calm seawater can be used to correct system calibration errors. However, most historic surveys consist only of data collected over land, where the conductive half-space assumption is rarely justified. We estimate the required recalibration parameters by analyzing systematic misfits in the inversion of statistically chosen measures of historic data. Our method requires the identification, within the survey area, of a zone of conductive responses that are reasonably uniform. From this zone, a set of altitude-corrected median responses are estimated. These are inverted using geologically specifiedconstraints to obtain a best-fit layered earth model. Systematic inconsistencies between the median measured altitude and the inverted depth to surface are attributed to altitude error. Remaining frequency-dependent fitting errors are assumed to be the calibration errors. We tested the method with partial success on helicopter electromagnetic data sets collected over uniform deep sediments where seawater data were also available and two different inland surveys over multiple lithologies in one general area. At high frequencies, our method works reliably. Recalibration of low-frequency data is not possible if the area used as a reference consists of moderate or poor conductors. In this case, data amplitudes are small and are greatly affected by imperfect drift and magnetic susceptibility corrections. Historic helicopter electromagnetic data may require amplitude rescaling up to 20%–30%, with phase shifts of up to 3°.

A limited-memory quasi-Newton inversion for 1D magnetotellurics

We apply a limited-memory quasi-Newton (QN) method to the 1D magnetotelluric (MT) inverse problem. Using this method we invert a realistic synthetic MT impedance data set calculated for a layered earth model. The calculation of gradients based on the adjoint method speeds up the inverse problem solution many times. In addition, regularization stabilizes the QN inversion result and a few correction pairs are sufficient to produce reasonable results. Comparison with the L-BFGS-B algorithm shows similar convergence rates. This study is a first step towards the solution of large-scale electromagnetic problems, with a full treatment of the 3D conductivity structure of the earth.

The rotational and gravitational signature of the December 26, 2004 Sumatran earthquake

Besides generating seismic waves, which eventually dissipate, an earthquake also generates a static displacement field everywhere within the Earth. This global displacement field rearranges the Earth’s mass thereby causing the Earth’s rotation and gravitational field to change. The size of this change depends upon the magnitude, focal mechanism, and location of the earthquake. The Sumatran earthquake of December 26, 2004 is the largest earthquake to have occurred since the 1960 Chilean earthquake. Using a spherical, layered Earth model, the coseismic effect of the Sumatran earthquake upon the Earth’s length-of-day, polar motion, and low-degree harmonic coefficients of the gravitational field are computed. Using a model of the earthquake source that is composed of five subevents having a total moment-magnitude Mw of 9.3, it is found that this earthquake should have caused the length-of-day to decrease by 6.8 microseconds, the position of the Earth’s generalized figure axis to shift 2.32 milliarcseconds towards 127° E longitude, the Earth’s oblateness J2 to decrease by 2.37 × 10−11 and the Earth’s pear-shapedness J3 to decrease by 0.63 × 10−11. The predicted change in the length-of-day, position of the generalized figure axis, and J3 are probably not detectable by current measurement systems. But the predicted change in oblateness is perhaps detectable if other effects, such as those of the atmosphere, oceans, and continental water storage, can be adequately removed from the observations.

Removing free-surface multiples from teleseismic transmission and constructed reflection responses using reciprocity and the inverse scattering series

We develop a new way to remove free-surface multiples from teleseismic P- transmission and constructed reflection responses. We consider two types of teleseismic waves with the presence of the free surface: One is the recorded waves under the real transmission geometry; the other is the constructed waves under a virtual reflection geometry. The theory presented is limited to 1D plane wave acoustic media, but this approximation is reasonable for the teleseismic P-wave problem resulting from the steep emergence angle of the wavefield. Using one-way wavefield reciprocity, we show how the teleseismic reflection responses can be reconstructed from the teleseismic transmission responses. We use the inverse scattering series to remove free-surface multiples from the original transmission data and from the reconstructed reflection response. We derive an alternative algorithm for reconstructing the reflection response from the transmission data that is obtained by taking the difference between the teleseismic transmission waves before and after free-surface multiple removal. Numerical tests with 1D acoustic layered earth models demonstrate the validity of the theory we develop. Noise test shows that the algorithm can work with S/N ratio as low as 5 compared to actual data with S/N ratio from 30 to 50. Testing with elastic synthetic data indicates that the acoustic algorithm is still effective for small incidence angles of typical teleseismic wavefields.

Fast computation of Hankel Transform using orthonormal exponential approximation of complex kernel function

The computation of electromagnetic (EM) fields, for 1-D layered earth model, requires evaluation of Hankel Transform (HT) of the EM kernel function. The digital filtering is the most widely used technique to evaluate HT integrals. However, it has some obvious shortcomings. We present an alternative scheme, based on an orthonormal exponential approximation of the kernel function, for evaluating HT integrals. This approximation of the kernel function was chosen because the analytical solution of HT of an exponential function is readily available in literature. This expansion reduces the integral to a simple algebraic sum. The implementation of such a scheme requires that the weights and the exponents of the exponential function be estimated. The exponents were estimated through a guided search algorithm while the weights were obtained using Marquardt matrix inversion method. The algorithm was tested on analytical HT pairs available in literature. The results are compared with those obtained using the digital filtering technique with Anderson filters. The field curves for four types (A-, K-, H-and Q-type) of 3-layer earth models are generated using the present scheme and compared with the corresponding curves obtained using the Anderson sc heme. It is concluded that the present scheme is more accurate than the Anderson scheme