First-arrival Traveltime Research Articles

Source-receiver Marchenko redatuming: Obtaining virtual receivers and virtual sources in the subsurface

By solving the Marchenko equations, one can retrieve the Green’s function (Marchenko Green’s function) between a virtual receiver in the subsurface and points at the surface (no physical receiver is required at the virtual location). We extend the idea behind these equations to retrieve the Green’s function between any two points in the subsurface, i.e., between a virtual source and a virtual receiver (no physical source or physical receiver is required at either of these locations). This Green’s function is called the virtual Green’s function, and it includes all primary, internal, and free-surface multiples. Similar to the Marchenko Green’s function, this virtual Green’s function requires the reflection response at the surface (single-sided illumination) and an estimate of the first-arrival traveltime from the virtual locations to the surface. These Green’s functions can be used to image the interfaces from above and below.

Slope tomography based on eikonal solvers and the adjoint-state method

Velocity macromodel building is a crucial step in the seismic imaging workflow as it provides the necessary background model for migration or full waveform inversion. In this study, we present a new formulation of stereotomography that can handle more efficiently long-offset acquisition, complex geological structures and large-scale data sets. Stereotomography is a slope tomographic method based upon a semi-automatic picking of local coherent events. Each local coherent event, characterized by its two-way traveltime and two slopes in common-shot and common-receiver gathers, is tied to a scatterer or a reflector segment in the subsurface. Ray tracing provides a natural forward engine to compute traveltime and slopes but can suffer from non-uniform ray sampling in presence of complex media and long-offset acquisitions. Moreover, most implementations of stereotomography explicitly build a sensitivity matrix, leading to the resolution of large systems of linear equations, which can be cumbersome when large-scale data sets are considered. Overcoming these issues comes with a new matrix-free formulation of stereotomography: a factored eikonal solver based on the fast sweeping method to compute first-arrival traveltimes and an adjoint-state formulation to compute the gradient of the misfit function. By solving eikonal equation from sources and receivers, we make the computational cost proportional to the number of sources and receivers while it is independent of picked events density in each shot and receiver gather. The model space involves the subsurface velocities and the scatterer coordinates, while the dips of the reflector segments are implicitly represented by the spatial support of the adjoint sources and are updated through the joint localization of nearby scatterers. We present an application on the complex Marmousi model for a towed-streamer acquisition and a realistic distribution of local events. We show that the estimated model, built without any prior knowledge of the velocities, provides a reliable initial model for frequency-domain FWI of long-offset data for a starting frequency of 4 Hz, although some artefacts at the reservoir level result from a deficit of illumination. This formulation of slope tomography provides a computationally efficient alternative to waveform inversion method such as reflection waveform inversion or differential-semblance optimization to build an initial model for pre-stack depth migration and conventional FWI.

Model regularization for seismic traveltime tomography with an edge-preserving smoothing operator

The solutions of the seismic first-arrival traveltime tomography are generally non-unique, and the Tikhonov model regularization for the inversion is commonly used to stabilize the inversion. However, the Tikhonov regularization for traveltime tomography often produces a low-resolution velocity model. To sharpen the velocity edges for the traveltime tomographic results and fit data at the same time, we should apply the edge-preserving concepts to regularize the inversion. In this study, we develop a new model regularization method by introducing an edge-preserving smoothing operator to detect the model edges in traveltime tomography. This edge-preserving smoothing operator is previously used in processing seismic images for enhancing resolution. We design three synthetic velocity models with sharp interfaces and with or without velocity gradients to study the performance of the regularization method with the edge-preserving smoothing operator. The new edge-preserving regularization not only sharpens the model edges but also maintains the smoothness of the velocity gradient in the layer.

Landstreamer seismics and physical property measurements in the Siilinjärvi open-pit apatite (phosphate) mine, central Finland

We tested the applicability of a newly developed broadband (0–800 Hz) digital-based seismic landstreamer for open-pit mine planning in the apatite-bearing Siilinjärvi mine in central Finland. Four seismic profiles, in total approximately 2.5 km long (2–4 m source and landstreamer receiver spacing), two inside the pit and two on its margins, were acquired in combination with wireless recorders connected to 10 Hz geophones and fixed at every 10 m spacing along the seismic profiles while the streamer data were being acquired. Downhole logging and laboratory physical property measurements on core and rock samples were carried out to not only support the seismic interpretations but also to provide information about the possible geophysical signature of these unique types of deposits. In spite of a highly noisy mining environment, seismic data of high quality were acquired; however, reflection processing and interpretations were challenged by the geologic complexities of several generations of basic and carbonatite dikes. To complement the reflection data imaging and to account for the steep elevation changes and crookedness of some of the seismic profiles, 3D first-arrival traveltime tomography and 3D swath reflection imaging were also carried out. Clear refracted arrivals from the open-pit profiles suggest the possibility of low-velocity zones associated with either blasting or several shear zones intersecting the seismic profiles. In terms of reflectivity, reflections have a different appearance from short and flat to longer and steep ones. The downhole- and borehole logging data suggest that some of these reflections are associated with diabase dikes and some are likely from zones of weaknesses in the alkaline-carbonatite complex. We determine the potential of using seismic streamers for cost- and time-effective open-pit mine planning and encourage further testing in simpler geologic settings to be established.

Joint Traveltime and Waveform Envelope Inversion for Near-surface Imaging

In the linearized seismic full waveform inversion (FWI), it is well known that a good initial model is needed to avoid cycle-skipping issue, especially when the low-frequency components of data are lacking. To solve the problem, we develop a joint first-arrival traveltime and early arrival envelope inversion method (JTE) to build a good near-surface velocity model with low-wavenumber components. The JTE method combines the robustness of nonlinear traveltime inversion and the low-frequency information embedded in the waveform envelope. In two synthetic experiments, we confirm that traveltime inversion constrains the top near-surface velocity structures tightly, while the envelope inversion recovers the low-wavenumber structures with low-velocity objects or layers. Using the results of JTE as starting models, we find that FWI can produce accurate solutions for complex numerical models. In a real-data example, JTE followed by FWI resolves a near-surface velocity model to improve statics corrections for the subsurface stacking image.

The use of seismic modeling for the geologic interpretation of deep seismic reflection data with low signal-to-noise ratios

We have performed a deep seismic reflection study, DACIA-PLAN, based on the data recorded along a crooked line across the southeastern Romanian Carpathians. The signal-to-noise ratio (S/N) of these data varies along the seismic profile, and its variation is considered to be an effect of the rough topography, complex subsurface geology, and varying surface conditions encountered during seismic data acquisition. The migrated time section that covers the mountainous area is clear, without visible reflections, making the geologic interpretation very difficult. We used a seismic modeling technique to explain the poor S/N of the recorded data and to generate synthetic seismic sections that can be useful for the geologic interpretation of the field seismic section (migrated time section). We used ray-tracing modeling to obtain the expected seismic expression of horizons of interest. Subsurface illumination modeling indicates that the complex subsurface geology and irregularly deployed sources and receivers are responsible for the incomplete and/or uneven illumination of the subsurface and can lead to strong amplitude variations. We then used 2.5D acoustic finite-difference modeling to analyze the effect of a crooked line on seismic wave propagation. The synthetic shot gathers prove that crooked line arrival times for reflected and head waves contain static time shifts relative to a straight line regular sampling geometry. Some geologic interfaces of interest are not well-imaged on the synthetic seismic section, and this is considered to be an effect of poor positioning during seismic data acquisition. We used the velocity model from the tomographic inversion of first-arrival traveltimes and synthetic and field crooked line deep seismic reflection data to create a structural image for the southeastern Romanian Carpathians and the Focsani Basin, which tie well with the geologic model built for this area on the basis of geologic and well data only.

An efficient 3D traveltime calculation using coarse-grid mesh for shallow-depth source

3D Kirchhoff pre-stack depth migration requires an efficient algorithm to compute first-arrival traveltimes. In this paper, we exploited a wave-equation-based traveltime calculation algorithm, which is called the suppressed wave equation estimation of traveltime (SWEET), and the equivalent source distribution (ESD) algorithm. The motivation of using the SWEET algorithm is to solve the Laplace-domain wave equation using coarse grid spacing to calculate first-arrival traveltimes. However, if a real source is located at shallow-depth close to free surface, we cannot accurately calculate the wavefield using coarse grid spacing. So, we need an additional algorithm to correctly simulate the shallow source even for the coarse grid mesh. The ESD algorithm is a method to define a set of distributed nodal sources that approximate a point source at the inter-nodal location in a velocity model with large grid spacing. Thanks to the ESD algorithm, we can efficiently calculate the first-arrival traveltimes of waves emitted from shallow source point even when we solve the Laplace-domain wave equation using a coarse-grid mesh. The proposed algorithm is applied to the SEG/EAGE 3D salt model. From the result, we note that the combination of SWEET and ESD algorithms can be successfully used for the traveltime calculation under the condition of a shallow-depth source. We also confirmed that our algorithm using coarse-grid mesh requires less computational time than the conventional SWEET algorithm using relatively fine-grid mesh.

First-arrival traveltime tomography for anisotropic media using the adjoint-state method

Traveltime tomography using transmission data has been widely used for static corrections and for obtaining near-surface models for seismic depth imaging. More recently, it is also being used to build initial models for full-waveform inversion. The classic traveltime tomography approach based on ray tracing has difficulties in handling large data sets arising from current seismic acquisition surveys. Some of these difficulties can be addressed using the adjoint-state method, due to its low memory requirement and numerical efficiency. By coupling the gradient computation to nonlinear optimization, it avoids the need for explicit computation of the Fréchet derivative matrix. Furthermore, its cost is equivalent to twice the solution of the forward-modeling problem, irrespective of the size of the input data. The presence of anisotropy in the subsurface has been well established during the past few decades. The improved seismic images obtained by incorporating anisotropy into the seismic processing workflow justify the effort. However, previous literature on the adjoint-state method has only addressed the isotropic approximation of the subsurface. We have extended the adjoint-state technique for first-arrival traveltime tomography to vertical transversely isotropic (VTI) media. Because [Formula: see text] is weakly resolvable from surface seismic alone, we have developed the mathematical framework and procedure to invert for [Formula: see text] and [Formula: see text]. Our numerical tests on the VTI SEAM model demonstrate the ability of the algorithm to invert for near-surface model parameters and reveal the accuracy achievable by the algorithm.

Edge-preserving traveltime tomography with a sparse multiscale imaging constraint

Solving the near-surface statics problem is often the first step in land or shallow marine seismic data processing. Near-surface velocity structures can be very complex, with large velocity contrasts within a small depth range. First-arrival traveltime tomography is a common approach for near-surface imaging. However, first-arrival traveltime tomography generally produces smooth model solutions due to the Tikhonov regularization, which constrains the model for minimum structures. Failing to resolve high velocity contrasts may result in inaccurate static values for reflection imaging. In this study, we develop a sparse multiscale imaging constraint for traveltime tomography to address this issue. In this method, we assume that the velocity model is sparse under a known wavelet basis. According to the model sparse representation, we first obtain the low wavenumber velocity structures, followed by the finer features, by alternately solving two sets of inversion problems. The synthetic tests and two real data applications show that this method exhibits better performance in reconstructing near-surface models with high velocity contrasts.

A two-dimensional analysis of the sensitivity of a pulse first break to wave speed contrast on a scale below the resolution length of ray tomography.

Mapping the speed of mechanical waves traveling inside a medium is a topic of great interest across many fields from geoscience to medical diagnostics. Much work has been done to characterize the fidelity with which the geometrical features of the medium can be reconstructed and multiple resolution criteria have been proposed depending on the wave-matter interaction model used to decode the wave speed map from scattering measurements. However, these criteria do not define the accuracy with which the wave speed values can be reconstructed. Using two-dimensional simulations, it is shown that the first-arrival traveltime predicted by ray theory can be an accurate representation of the arrival of a pulse first break even in the presence of diffraction and other phenomena that are not accounted for by ray theory. As a result, ray-based tomographic inversions can yield accurate wave speed estimations also when the size of a sound speed anomaly is smaller than the resolution length of the inversion method provided that traveltimes are estimated from the signal first break. This increased sensitivity however renders the inversion more susceptible to noise since the amplitude of the signal around the first break is typically low especially when three-dimensional anomalies are considered.

Tomogram-based comparison of geostatistical models: Application to the Macrodispersion Experiment (MADE) site

We propose a new framework to compare alternative geostatistical descriptions of a given site. Multiple realizations of each of the considered geostatistical models and their corresponding tomograms (based on inversion of noise-contaminated simulated data) are used as a multivariate training image. The training image is scanned with a direct sampling algorithm to obtain conditional realizations of hydraulic conductivity that are not only in agreement with the geostatistical model, but also honor the spatially varying resolution of the site-specific tomogram. Model comparison is based on the quality of the simulated geophysical data from the ensemble of conditional realizations. The tomogram in this study is obtained by inversion of cross-hole ground-penetrating radar (GPR) first-arrival travel time data acquired at the MAcro-Dispersion Experiment (MADE) site in Mississippi (USA). Various heterogeneity descriptions ranging from multi-Gaussian fields to fields with complex multiple-point statistics inferred from outcrops are considered. Under the assumption that the relationship between porosity and hydraulic conductivity inferred from local measurements is valid, we find that conditioned multi-Gaussian realizations and derivatives thereof can explain the crosshole geophysical data. A training image based on an aquifer analog from Germany was found to be in better agreement with the geophysical data than the one based on the local outcrop, which appears to under-represent high hydraulic conductivity zones. These findings are only based on the information content in a single resolution-limited tomogram and extending the analysis to tracer or higher resolution surface GPR data might lead to different conclusions (e.g., that discrete facies boundaries are necessary).

An iterative, fast-sweeping-based eikonal solver for 3D tilted anisotropic media

Computation of first-arrival traveltimes for quasi-P waves in the presence of anisotropy is important for high-end near-surface modeling, microseismic-source localization, and fractured-reservoir characterization — and it requires solving an anisotropic eikonal equation. Anisotropy deviating from elliptical anisotropy introduces higher order nonlinearity into the eikonal equation, which makes solving the eikonal equation a challenge. We addressed this challenge by iteratively solving a sequence of simpler tilted elliptically anisotropic eikonal equations. At each iteration, the source function was updated to capture the effects of the higher order nonlinear terms. We used Aitken’s extrapolation to speed up convergence rate of the iterative algorithm. The result is an algorithm for computing first-arrival traveltimes in tilted anisotropic media. We evaluated the applicability and usefulness of our method on tilted transversely isotropic media and tilted orthorhombic media. Our numerical tests determined that the proposed method matches the first arrivals obtained by wavefield extrapolation, even for strongly anisotropic and highly complex subsurface structures. Thus, for the cases where two-point ray tracing fails, our method can be a potential substitute for computing traveltimes. The approach presented here can be easily extended to compute first-arrival traveltimes for anisotropic media with lower symmetries, such as monoclinic or even the triclinic media.

Fast Huygens’ sweeping methods for multiarrival Green’s functions of Helmholtz equations in the high-frequency regime

Multiarrival Green’s functions are essential in seismic modeling, migration, and inversion. Huygens-Kirchhoff (HK) integrals provide a bridge to integrate locally valid first-arrival Green’s functions into a globally valid multiarrival Green’s function. We have designed robust and accurate finite-difference methods to compute first-arrival traveltimes and amplitudes, so that first-arrival Green’s functions can be constructed rapidly. We adapted a fast butterfly algorithm to evaluate discretized HK integrals. The resulting fast Huygens’ sweeping method has the following unique features: (1) it precomputes a set of local traveltime and amplitude tables, (2) it automatically takes care of caustics, (3) it constructs Green’s functions of the Helmholtz equation for arbitrary frequencies and for many point sources, and (4) for a fixed number of points per wavelength, it constructs each Green’s function in nearly optimal complexity [Formula: see text] in terms of the total number of mesh points [Formula: see text], where the prefactor of the complexity only depends on the specified accuracy, and is independent of the frequency. The 2D and 3D examples revealed the performance of the method.

A study of fuzzy c-means coupling for joint inversion, using seismic tomography and gravity data test scenarios

Joint inversion, the inversion of multiple geophysical data sets containing complementary information about the subsurface, has the potential to significantly improve inversion results by reducing the nonuniqueness of the inverse problem. One of the challenges of joint inversion is deciding how to couple the multiple physical property models. If a coupling approach is used that is inconsistent with the physical truth, then inversion artifacts can occur and may lead to incorrect interpretations. In this paper, we investigated the fuzzy c-means (FCM) clustering approach to provide a lithological coupling of the seismic velocity and density models in joint 2D inversions of first-arrival traveltimes and gravity data. Even though this coupling approach has been used in previous works, recommendations for its effective use have not yet been developed. We conducted a suite of joint inversion tests on synthetic data generated from a geologically realistic model based on magmatic massive sulfide deposits. There is a known relationship between seismic velocity and density for the silicate rocks and sulfide minerals involved; this lithological relationship was used to design a clustered coupling strategy in the joint inversions. The tests we conducted clearly exhibited the benefits of joint inversion using FCM coupling. Our work revealed the effects of including inaccurate a priori physical property information. We also evaluated approaches to assess whether such inaccurate information may have been used.

Topography-Dependent Eikonal Traveltime Tomography for Upper Crustal Structure Beneath an Irregular Surface

Seismic modeling of the crust with nonflat topography can be made by first-arrival traveltime tomography, which faces the challenge of an irregular free surface. A feasible way to deal with this problem consists of expanding the physical space by overlapping a low velocity layer above the irregular surface in order to have a flat topography, besides using the classical eikonal equation solver for traveltime computation. However, the undesirable consequences of this method include seismic ray deviations due to the transition from an irregular surface that is the free boundary to an inner discontinuity lying in the expanded computational space. An alternative solution, called irregular surface flattening, which involves the transformation between curvilinear and Cartesian coordinate systems, has been recently proposed through the formulation of the topography-dependent eikonal equation (TDEE) and a new solver for forward modeling of traveltimes. Based on the solution of this equation, we present topography-dependent eikonal traveltime tomography (hereafter TDETT) for seismic modeling of the upper crust. First-arrival traveltimes are calculated using the TDEE solver and the raypaths with the minimum traveltime that can be found by following the steepest traveltime gradient from the receiver to the source. By solving an algebraic equation system that connects the slowness perturbations with the already determined traveltimes, these variables can be obtained making use of the back-projection algorithm. This working scheme is evaluated through three numerical examples with different topographic complexities that are conducted from synthetic data and a fourth example with somewhat more complicated topography and real data acquired in northeastern Tibet. The comparison of the results obtained by both methods, i.e., physical space expansion above the irregular surface and irregular surface flattening, fully validates the tomography scheme that is proposed to construct seismic velocity models with nonflat topography.

First-arrival traveltime sound speed inversion with a priori information.

A first-arrival travel-time sound speed algorithm presented by Tarantola [Inverse Problem Theory and Methods for Model Parameter Estimation (SIAM, Philadelphia, PA, 2005)] is adapted to the medical ultrasonics setting. Through specification of a covariance matrix for the object model, the algorithm allows for natural inclusion of physical a priori information of the object. The algorithm's ability to accurately and robustly reconstruct a complex sound speed distribution is demonstrated on simulation and experimental data using a limited aperture. The algorithm is first demonstrated generally in simulation with a numerical breast phantom imaged in different geometries. As this work is motivated by the authors' limited aperture dual sided ultrasound breast imaging system, experimental data are acquired with a Verasonics system with dual, 128 element, linear L7-4 arrays. The transducers are automatically calibrated for usage in the eikonal forward model.A priori information such as knowledge of correlated regions within the object is obtained via segmentation of B-mode images generated from synthetic aperture imaging. As one illustration of the algorithm's facility for inclusion ofa priori information, physically grounded regularization is demonstrated in simulation. The algorithm's practicality is then demonstrated through experimental realization in limited aperture cases. Reconstructions of sound speed distributions of various complexity are improved through inclusion of a priori information. The sound speed maps are generally reconstructed with accuracy within a few m/s. This paper demonstrates the ability to form sound speed images using two opposed commercial linear arrays to mimic ultrasound image acquisition in the compressed mammographic geometry. The ability to create reasonably good speed of sound images in the compressed mammographic geometry allows images to be readily coregistered to tomosynthesis image volumes for breast cancer detection and characterization studies.

Seismic travel-time and attenuation tomography to characterize the excavation damaged zone and the surrounding rock mass of a newly excavated ramp and chamber

Seismic travel-time and attenuation tomography were applied to characterize the excavation damaged zone and the adjacent rock mass in the GFZ-Underground-Lab within the research and education mine Reiche Zeche of the Technical University Bergakademie Freiberg (Germany). The lab is situated in gneiss rocks at 150m depth and comprises three galleries which enclose an area of approximately 50m×100m. Along these galleries two seismic surveys were performed before and after the excavation of a new ramp and chamber. For both measurements, travel-time and attenuation tomographies for P-waves were performed with the ray-based inversion algorithm SIMULPS14. The seismic velocities were calculated from first-arrival travel times whereas a logarithmic-spectral-ratio approach was used to calculate the corresponding quality factors (Q) for attenuation tomography. A comparison of the tomograms reveals a decrease of average P-wave velocity values from 5.64km/s to 5.54km/s and of average Q-values from 29.8 to 26.5 in the whole area after the excavation of the new cavities. The maximum changes are located at already weakened zones either at the conjunction of two galleries or of a major fracture zone with a gallery. The attenuation tomography shows a higher sensitivity to rock mass changes than the travel-time tomography. However, the calculation of the Q-values demands a higher signal quality than the determination of the seismic travel times.

A new model of the upper mantle structure beneath the western rim of the East European Craton

Abstract. We present a new 1-D P wave seismic velocity model (called MP1-SUW) of the upper mantle structure beneath the western rim of the East European Craton (EEC) based on the analysis of the earthquakes recorded at the Suwałki (SUW) seismic station located in NE Poland which belongs to the Polish Seismological Network (PLSN). Motivation for this study arises from the observation of a group of reflected waves after expected P410P at epicentral distances 2300–2800 km from the SUW station. Although the existing global models represent the first-arrival traveltimes, they do not represent the full wavefield with all reflected waves because they do not take into account the structural features occurring regionally such as 300 km discontinuity. We perform P wave traveltime analysis using 1-D and 2-D forward ray-tracing modelling for the distances of up to 3000 km. We analysed 249 natural seismic events from four azimuthal spans with epicentres in the western Mediterranean Sea region (WMSR), the Greece and Turkey region (GTR), the Caucasus region (CR) and the part of the northern Mid-Atlantic Ridge near the Jan Mayen Island (JMR). For all chosen regions, except the JMR group for which 2-D modelling was performed, we estimate a 1-D average velocity model which will characterize the main seismic discontinuities. It appears that a single 1-D model (MP1-SUW model) explains well the observed traveltimes for the analysed groups of events. Differences resulting from the different azimuth range of earthquakes are close to the assumed picking uncertainty. The MP1-SUW model documents the bottom of the asthenospheric low-velocity zone (LVZ) at the depth of 220 km, 335 km discontinuity and the zone with the reduction of P wave velocity atop 410 km discontinuity which is depressed to 440 km depth. The nature of the regionally occurring 300 km boundary is explained here by tracing the ancient subduction regime related to the closure of the Iapetus Ocean, the Rheic Ocean and the Tornquist Sea.

Accounting for imperfect forward modeling in geophysical inverse problems — Exemplified for crosshole tomography

Inversion of geophysical data relies on knowledge about how to solve the forward problem, that is, computing data from a given set of model parameters. In many applications of inverse problems, the solution to the forward problem is assumed to be known perfectly, without any error. In reality, solving the forward model (forward-modeling process) will almost always be prone to errors, which we referred to as modeling errors. For a specific forward problem, computation of crosshole tomographic first-arrival traveltimes, we evaluated how the modeling error, given several different approximate forward models, can be more than an order of magnitude larger than the measurement uncertainty. We also found that the modeling error is strongly linked to the spatial variability of the assumed velocity field, i.e., the a priori velocity model. We discovered some general tools by which the modeling error can be quantified and cast into a consistent formulation as an additive Gaussian observation error. We tested a method for generating a sample of the modeling error due to using a simple and approximate forward model, as opposed to a more complex and correct forward model. Then, a probabilistic model of the modeling error was inferred in the form of a correlated Gaussian probability distribution. The key to the method was the ability to generate many realizations from a statistical description of the source of the modeling error, which in this case is the a priori model. The methodology was tested for two synthetic ground-penetrating radar crosshole tomographic inverse problems. Ignoring the modeling error can lead to severe artifacts, which erroneously appear to be well resolved in the solution of the inverse problem. Accounting for the modeling error leads to a solution of the inverse problem consistent with the actual model. Further, using an approximate forward modeling may lead to a dramatic decrease in the computational demands for solving inverse problems.

Calculating Ray Paths for First-Arrival Travel Times Using a Topography-Dependent Eikonal Equation Solver

Abstract Irregular surfaces cause a number of problems for seismic processing and interpretation. The problems lie in the proper treatment of topography in the first‐arrival travel‐time calculation and ray‐path tracing, both of which are subject to preconditions in ray‐based seismogram synthesis, seismic tomography, and seismic migration calculations. Two treatment schemes for irregular surfaces have been used previously: (1) model expansion with the irregular surface treated as an inner discontinuity, and (2) flattening of the irregular surface using a transformation between curvilinear and Cartesian coordinates, while maintaining it as a free surface. In the first approach, first‐arrival travel times can be calculated using an eikonal equation solver, and rays are traced backward from the receiver to the source along the direction of the gradient of the travel‐time field. In the second scheme, a topography‐dependent eikonal equation is used to calculate irregular‐surface first‐arrival travel times. We present a ray‐path tracing scheme for irregular surfaces, which applies a travel‐time field calculated using a topography‐dependent eikonal equation. The scheme is realized using travel‐time gradients in curvilinear coordinates. The validity of the scheme for tracing ray paths in the presence of irregular topography is illustrated by five models containing differing degrees of topographical complexity. Comparison of ray paths and first‐arrival travel times between the two irregular‐surface treatment schemes suggests the topography‐flattened scheme avoids the difficulties of both discretizing the irregular surface and the velocity selection of the infill medium in the model expansion scheme. For the treatment of the inner discontinuity by model expansion, there is a possibility that ray paths will be traced outside the real physical model. Using the topography‐dependent eikonal equation solver with our ray‐path tracing scheme should provide an efficient means of dealing with irregular surfaces, which could be applied in the fields of seismic tomography, seismic migration, and tomographic static correction.