Transformation relations for second-order derivatives of travel time in anisotropic media

In the computation of paraxial travel times and Gaussian beams, the basic role is played by the second-order derivatives of the travel-time field at the reference ray. These derivatives can be determined by dynamic ray tracing (DRT) along the ray. Two basic DRT systems have been broadly used in applications: the DRT system in Cartesian coordinates and the DRT system in ray-centred coordinates. In this paper, the transformation relations between the second-order derivatives of the travel-time field in Cartesian and ray-centred coordinates are derived. These transformation relations can be used both in isotropic and anisotropic media, including computations of complex-valued travel times necessary for the evaluation of Gaussian beams.     

Calculating complex‐valued P‐wave first‐arrival traveltimes in attenuative vertical transversely isotropic media with an irregular surface

The complex‐valued first‐arrival traveltime can be used to describe the properties of both velocity and attenuation as seismic waves propagate in attenuative elastic media. The real part of the complex‐valued traveltime corresponds to phase arrival and the imaginary part is associated with the amplitude decay due to energy absorption. The eikonal equation for attenuative vertical transversely isotropic media discretized with rectangular grids has been proven effective and precise to calculate the complex‐valued traveltime, but less accurate and efficient for irregular models. By using the perturbation method, the complex‐valued eikonal equation can be decomposed into two real‐valued equations, namely the zeroth‐ and first‐order traveltime governing equations. Here, we first present the topography‐dependent zeroth‐ and first‐order governing equations for attenuative VTI media, which are obtained by using the coordinate transformation from the Cartesian coordinates to the curvilinear coordinates. Then, we apply the Lax–Friedrichs sweeping method for solving the topography‐dependent traveltime governing equations in order to approximate the viscosity solutions, namely the real and imaginary parts of the complex‐valued traveltime. Several numerical tests demonstrate that the proposed scheme is efficient and accurate in calculating the complex‐valued P‐wave first‐arrival traveltime in attenuative VTI media with an irregular surface.     

Ray-centred coordinate systems in anisotropic media

Whereas the ray-centred coordinates for isotropic media by Popov and Pšenčík are uniquely defined by the selection of the basis vectors at one point along the ray, there is considerable freedom in selecting the ray-centred coordinates for anisotropic media. We describe the properties common to all ray-centred coordinate systems for anisotropic media and general conditions, which may be imposed on the basis vectors. We then discuss six different particular choices of ray-centred coordinates in an anisotropic medium. This overview may be useful in choosing the ray-centred coordinates best suited for a particular application. The equations are derived for a general homogeneous Hamiltonian of an arbitrary degree and are thus applicable both to the anisotropic-ray-theory rays and anisotropic common S-wave rays.     

Second-Order and Higher-Order Perturbations of Travel Time in Isotropic and Anisotropic Media

The partial derivatives of travel time with respect to model parameters are referred to as perturbations. Explicit equations for the second-order and higher-order perturbations of travel time in both isotropic and anisotropic media are derived. The perturbations of travel time and its spatial derivatives can be calculated by simple numerical quadratures along rays.     

Review of the Anisotropic Interface Ray Propagator: Symplecticity,Eigenvalues, Invariants and Applications

A review of the 6 × 6 anisotropic interface ray propagator matrix in Cartesian coordinates and within the framework of the Hamiltonian formalism shows that there is one unique propagator satisfying the symplectic property. This is essential, since the symplecticity furnishes an exact inverse, while an eigenvalue analysis indicates that the propagator may be arbitrarily ill-conditioned. As such, the symplectic interface propagator naturally connects to symplectic ray integration algorithms for smooth media, designed to maintain accuracy. Moreover, several ray invariants for smooth media remain invariant across interfaces. It is straightforward to derive expressions for the interface propagator, both explicit and implicit. Symplecticity is equivalent to the condition that the propagator preserves the eikonal constraint across the interface. The symplectic interface propagator complies with phase matching of the incident and reflected/transmitted ray field, and is therefore in accordance with the earlier derived 4 × 4 matrix in ray-centred coordinates. The symplectic property is related to the symmetry of the second derivative matrix of the reflected/transmitted traveltime field. Thanks to the analytic expression of the symplectic interface propagator, relating interface curvature directly to second derivatives of traveltimes observed at a datum level, numerous applications are available in the area of processing and inversion.     

Paraxial ray methods in inhomogeneous anisotropic media: Initial conditions

Paraxial ray methods have found broad applications in the seismic ray method and in numerical modelling and interpretation of high-frequency seismic wave fields propagating in inhomogeneous, isotropic or anisotropic structures. The basic procedure in paraxial ray methods consists in dynamic ray tracing. We derive the initial conditions for dynamic ray equations in Cartesian coordinates, for rays initiated at three types of initial manifolds given in a three-dimensional medium: 1) curved surfaces (surface source), 2) isolated points (point source), and 3) curved, planar and non-planar lines (line source). These initial conditions are very general, valid for homogeneous or inhomogeneous, isotropic or anisotropic media, and for both a constant and a variable initial travel time along the initial manifold. The results presented in the paper considerably extend the possible applications of the paraxial ray method.     

Time‐domain analysis of unbounded media using mixed‐variable formulations

Formulation of a matrix‐valued force–displacement relationship which can take radiation damping into account is of major importance when modelling unbounded domains. This can be done by means of fundamental solutions in space and time in connection with convolution integrals or by means of a frequency dependent boundary element representation, but for discrete frequencies Ω only. In this paper a method for interpolating discrete values of dynamic stiffness matrices by a continuous matrix valued rational function is proposed. The coupling between interface degrees of freedom is fully preserved. Another crucial point in soil–structure interaction analysis is how to implement an approximation in the spectral domain into a time‐domain analysis. Well‐known approaches for the scalar case are based on the partial‐fraction expansion of a scalar rational function. Here, a more general procedure, applicable to MDOF‐systems, for the transformation of spectral rational approximations into the time‐domain is introduced. Evaluation of the partial‐fraction expansion is avoided by using the so‐called mixed variables. Thus, unknowns in the time‐domain are displacements as well as forces. Copyright © 2001 John Wiley & Sons, Ltd.     

Dynamic analyses of multilayered poroelastic media using the generalized transfer matrix method

Since the generalized transfer matrix method overcomes the intrinsic instability of the Thomson–Haskell transfer matrix method for both high frequencies and/or thick layers, it can produce stable and accurate solutions for the dynamic analysis of viscoelastic media as well as anisotropic media. This paper extends the generalized transfer matrix method to the dynamic analysis of multilayered poroelastic media. Main improvements include the presentation of the concisely explicit general solutions for the dynamic analysis of multilayered poroelastic media and the derivation of an analytical inversion of 8×8 order layer matrix corresponding to the general solutions. The explicit solutions are valid for the dynamic analysis of one-dimensional, two-dimensional and three-dimensional poroelastic medium problems. In addition, an efficient recursive scheme is proposed for accurate determination of the equivalent interface sources for multilayered poroelastic media due to excitation by a source at an arbitrary depth. For the dynamic analysis of multilayered poroelastic media, the generalized transfer matrix method recursively transfers both the interface stiffness matrix and equivalent source starting from the bottom half space to the top layer, without resort to the numerical solution of the assembled global equations as the exact stiffness matrix method does. While keeping the simplicity of the Thomson–Haskell transfer matrix method, the generalized transfer matrix method is both accurate and stable. The related numerical examples have demonstrated that the generalized transfer matrix method is an alternative approach to conducting the dynamic analysis of multilayered poroelastic media.     

Errors due to the anisotropic-common-ray approximation of the coupling ray theory

The common-ray approximation eliminates problems with ray tracing through S-wave singularities and also considerably simplifies the numerical algorithm of the coupling ray theory for S waves, but may introduce errors in travel times due to the perturbation from the common reference ray. These travel-time errors can deteriorate the coupling-ray-theory solution at high frequencies. It is thus of principal importance for numerical applications to estimate the errors due to the common-ray approximation applied. The anisotropic-common-ray approximation of the coupling ray theory is more accurate than the isotropic-common-ray approximation. We derive the equations for estimating the travel-time errors due to the anisotropic-common-ray (and also isotropic-common-ray) approximation of the coupling ray theory. The errors of the common-ray approximations are calculated along the anisotropic common rays in smooth velocity models without interfaces. The derivation is based on the general equations for the second-order perturbations of travel time.     

Wave propagation in semi-infinite media with topographical irregularities using Decoupled Equations Method

In this paper, a novel semi-analytical method, called Decoupled Equations Method (DEM), is presented for modeling of elastic wave propagation in the semi-infinite two-dimensional (2D) media which are involved surface topography. In the DEM, only the boundaries of the problem are discretized by specific subparametric elements, in which special shape functions as well as higher-order Chebyshev mapping functions are implemented. For the shape functions, Kronecker Delta property is satisfied for displacement function. Moreover, the first derivatives of displacement function with respect to the tangential coordinates on the boundaries are assigned to zero at any given node. Employing the weighted residual method and using Clenshaw–Curtis numerical integration, coefficient matrices of the system of equations are transformed into diagonal ones, which leads to a set of decoupled partial differential equations. To evaluate the accuracy of the DEM in the solution of scattering problem of plane waves, cylindrical topographical features of arbitrary shapes are solved. The obtained results present excellent agreement with the analytical solutions and the results from other numerical methods.     

Vibration of surface foundations of arbitrary shapes

Presented is a systematic procedure for generating impedance (or compliance) matrices for foundations with arbitrary shapes, resting on an elastic half-space medium. A technique to decompose prescribed harmonic tractions on the half-space medium is employed to solve analytically the differential wave equations in cylindrical coordinates. However, the interaction stresses due to the vibration of a foundation with arbitrary shape are described in rectangular coordinates, and assumed to be piecewise constant in the region of the arbitrary shape. A coordinate transformation matrix is introduced for the piecewise constant tractions in order to use the solution of the differential wave equations in cylindrical coordinates. Finite element modelling is assumed in rectangular coordinates for the foundation itself. The impedance matrix is then obtained for the finite element model, using a variational principle and the reciprocal theorem. A simple example of a rigid square plate resting on a half-space medium and subjected to vertical excitation is used to demonstrate the efficiency and effectiveness of the procedure. Some numerical aspects are investigated and some possible extensions of the procedure are also discussed.     

Introduction to seismic travel time methods in anisotropic media

Summary Section 1 (and 11) develops the concepts of the front velocity, the front gradient, the travel time in space and on seismometric profiles, the profile velocity and the profile gradient in connection with the propagation of the fronts of elastic waves in solid isotropic and anisotropic media. The sectional velocity and the sectional gradient are defined in terms of the motion of the curve of intersection of a front with a fixed surface. Section 2 (and 12) relates the coefficients of elasticity of the medium, the front types, and their respective rays. In section 12, the theory of fronts of arbitrary shape and of the corresponding rays for any anisotropic, homogeneous or inhomogeneous solid medium is summarized. In section 3 (and 13), the law of reflection and refraction of fronts on surfaces of discontinuity of arbitrary shape is presented. Sections 4 to 6 (and 14 to 16) treat some elementary applications of seismic travel time methods to homogeneous, uniaxially anisotropic media (=transverse isotropy) in greater detail. In section 4 (and 14), the travel time of a direct front generated by a point source is considered and it is shown how the coefficients of elasticity of the medium can be found based on travel time measurements. The seismic prospection of a plane reflector and of a reflecting boundary of arbitrary shape and position are discussed in section 5 (and 15). In section 6 (and 16), the seismic refraction method is used to locate a plane boundary between a homogeneous, uniaxially anisotropic and a homogeneous isotropic medium, where the boundary is perpendicular or at an arbitrary angle to the direction of anisotropy.     

A note on two-point paraxial travel times

Recently, several expressions for the two-point paraxial travel time in laterally varying, isotropic or anisotropic layered media were derived. The two-point paraxial travel time gives the travel time from point S′ to point R′, both these points being situated close to a known reference ray Ω, along which the ray-propagator matrix was calculated by dynamic ray tracing. The reference ray and the position of points S′ and R′ are specified in Cartesian coordinates. Two such expressions for the two-point paraxial travel time play an important role. The first is based on the 4 × 4 ray propagator matrix, computed by dynamic ray tracing along the reference ray in ray-centred coordinates. The second requires the knowledge of the 6 × 6 ray propagator matrix computed by dynamic ray tracing along the reference ray in Cartesian coordinates. Both expressions were derived fully independently, using different methods, and are expressed in quite different forms. In this paper we prove that the two expressions are fully equivalent and can be transformed into each other.     

A discontinuous Galerkin method for seismic wave propagation in coupled elastic and poroelastic media

Numerical simulation in coupled elastic and poroelastic media is important in oil and gas exploration. However, the interface between elastic and poroelastic media is a challenge to handle. In order to deal with the coupled model, the first-order velocity–stress wave equations are used to unify the elastic and poroelastic wave equations. In addition, an arbitrary high-order discontinuous Galerkin method is used to simulate the wave propagation in coupled elastic–poroelastic media, which achieves same order accuracy in time and space domain simultaneously. The interfaces between the two media are explicitly tackled by the Godunov numerical flux. The proposed forms of numerical flux can be used efficiently and conveniently to simulate the wave propagation at the interfaces of the coupled model and handle the absorbing boundary conditions properly. Numerical results on coupled elastic–poroelastic media with straight and curved interfaces are compared with those from a software that is based on finite element method and the interfaces are handled by boundary conditions, demonstrating the feasibility of the proposed scheme in dealing with coupled elastic–poroelastic media. In addition, the proposed method is used to simulate a more complex coupled model. The numerical results show that the proposed method is feasible to simulate the wave propagation in such a media and is easy to implement.     

2D inversion of refraction traveltime curves using homogeneous functions

A method using simple inversion of refraction traveltimes for the determination of 2D velocity and interface structure is presented. The method is applicable to data obtained from engineering seismics and from deep seismic investigations. The advantage of simple inversion, as opposed to ray‐tracing methods, is that it enables direct calculation of a 2D velocity distribution, including information about interfaces, thus eliminating the calculation of seismic rays at every step of the iteration process. The inversion method is based on a local approximation of the real velocity cross‐section by homogeneous functions of two coordinates. Homogeneous functions are very useful for the approximation of real geological media. Homogeneous velocity functions can include straight‐line seismic boundaries. The contour lines of homogeneous functions are arbitrary curves that are similar to one another. The traveltime curves recorded at the surface of media with homogeneous velocity functions are also similar to one another. This is true for both refraction and reflection traveltime curves. For two reverse traveltime curves, non‐linear transformations exist which continuously convert the direct traveltime curve to the reverse one and vice versa. This fact has enabled us to develop an automatic procedure for the identification of waves refracted at different seismic boundaries using reverse traveltime curves. Homogeneous functions of two coordinates can describe media where the velocity depends significantly on two coordinates. However, the rays and the traveltime fields corresponding to these velocity functions can be transformed to those for media where the velocity depends on one coordinate. The 2D inverse kinematic problem, i.e. the computation of an approximate homogeneous velocity function using the data from two reverse traveltime curves of the refracted first arrival, is thus resolved. Since the solution algorithm is stable, in the case of complex shooting geometry, the common‐velocity cross‐section can be constructed by applying a local approximation. This method enables the reconstruction of practically any arbitrary velocity function of two coordinates. The computer program, known as godograf , which is based on this theory, is a universal program for the interpretation of any system of refraction traveltime curves for any refraction method for both shallow and deep seismic studies of crust and mantle. Examples using synthetic data demonstrate the accuracy of the algorithm and its sensitivity to realistic noise levels. Inversions of the refraction traveltimes from the Salair ore deposit, the Moscow region and the Kamchatka volcano seismic profiles illustrate the methodology, practical considerations and capability of seismic imaging with the inversion method.     

Wave front construction in smooth media for prestack depth migration

We implemented a wave front construction algorithm specifically designed for smooth media for application to prestack depth migration. The highest priority was given to maximum computational speed to allow an extension of the techniques to 3D media. A simple grid-based model representation in combination with fast bilinear interpolation is used. It is shown that this procedure has no distorting effects on the ray tracing results for smooth media. In our implementation, wave front construction (WFC) has proven to be as fast as some of the recently developed methods for travel time computations. WFC has advantages over these methods, since amplitudes and other ray theoretical quantities are available, and it is not restricted to the calculation of only first arrivals. Thus, it meets the requirements for migration in complex media. Furthermore, WFC allows for introduction of a perturbation scheme for computing travel times for slightly varying models simultaneously. This has applications for, e.g., prestack velocity estimation techniques. The importance of later arrivals for migration in complex media is demonstrated by prestack images of the Marmousi data set.     

Application of a perfectly matched layer in seismic wavefield simulation with an irregular free surface

Recently, an effective and powerful approach for simulating seismic wave propagation in elastic media with an irregular free surface was proposed. However, in previous studies, researchers used the periodic condition and/or sponge boundary condition to attenuate artificial reflections at boundaries of a computational domain. As demonstrated in many literatures, either the periodic condition or sponge boundary condition is simple but much less effective than the well‐known perfectly matched layer boundary condition. In view of this, we intend to introduce a perfectly matched layer to simulate seismic wavefields in unbounded models with an irregular free surface. We first incorporate a perfectly matched layer into wave equations formulated in a frequency domain in Cartesian coordinates. We then transform them back into a time domain through inverse Fourier transformation. Afterwards, we use a boundary‐conforming grid and map a rectangular grid onto a curved one, which allows us to transform the equations and free surface boundary conditions from Cartesian coordinates to curvilinear coordinates. As numerical examples show, if free surface boundary conditions are imposed at the top border of a model, then it should also be incorporated into the perfectly matched layer imposed at the top‐left and top‐ right corners of a 2D model where the free surface boundary conditions and perfectly matched layer encounter; otherwise, reflections will occur at the intersections of the free surface and the perfectly matched layer, which is confirmed in this paper. So, by replacing normal second derivatives in wave equations in curvilinear coordinates with free surface boundary conditions, we successfully implement the free surface boundary conditions into the perfectly matched layer at the top‐left and top‐right corners of a 2D model at the surface. A number of numerical examples show that the perfectly matched layer constructed in this study is effective in simulating wave propagation in unbounded media and the algorithm for implementation of the perfectly matched layer and free surface boundary conditions is stable for long‐time wavefield simulation on models with an irregular free surface.     

Fractal travel time estimates for dispersive contaminants

Alternative fractional models of contaminant transport lead to a new travel time formula for arbitrary concentration levels. For an evolving contaminant plume in a highly heterogeneous aquifer, the new formula predicts much earlier arrival at low concentrations. Travel times of contaminant fronts and plumes are often obtained from Darcy's law calculations using estimates of average pore velocities. These estimates only provide information about the travel time of the average concentration (or peak, for contaminant pulses). Recently, it has been shown that finding the travel times of arbitrary concentration levels is a straightforward process, and equations were developed for other portions of the breakthrough curve for a nonreactive contaminant. In this paper, we generalize those equations to include alternative fractional models of contaminant transport.