Elastic Finite-difference Modeling Research Articles

Removing Time Dispersion from Elastic Wave Modeling with the pix2pix Algorithm Based on cGAN

The finite-difference (FD) method is one of the most commonly used numerical methods for elastic wave modeling. However, due to the difference approximation of the derivative, the time dispersion phenomenon cannot be avoided. This paper proposes the use of pix2pix algorithm based on a conditional generative adversarial network (cGAN) for removing time dispersion from elastic FD modeling. Firstly, we analyze the time dispersion of elastic wave FD modeling. Then, we discuss the pix2pix algorithm based on cGAN, improve the loss function of the pix2pix algorithm by introducing a Sobel operator, and analyze the parameter selection of the network model for the pix2pix algorithm. Finally, we verify the feasibility and effectiveness of the pix2pix algorithm in removing time dispersion from elastic wave FD modeling through testing some model simulation data.

Improving accuracy and efficiency in elastic finite difference modelling

Multicomponent technology has known to have a better imaging quality as compared to P-wave imaging. However, the cost of acquiring and processing the multicomponent data is more costly than P-wave imaging. Although P-wave imaging may have poorer imaging capability as compared to shear wave imaging, a good P-wave imaging can still be obtained if a good seismic modelling technique is used. Thus, this paper aims to enhance seismic modelling technique, in terms of accuracy, suppression of grid dispersion, and reduction in computation time. The finite difference (FD) technique is used as the modelling scheme, and the elastic modelling approach is used in this paper. The new method is applied on homogeneous model and synthetic data, and then compared to conventional FD modelling (4th order), higher order FD modelling (14th order), optimization of staggered grid coefficients such as least squares method and Wang & Liang’s method. The comparison will then be shown and discussed.

Low-frequency ambient distributed acoustic sensing (DAS): case study from Perth, Australia

SUMMARYAmbient wavefield data acquired on existing (so-called ‘dark fibre’) optical fibre networks using distributed acoustic sensing (DAS) interrogators allow users to conduct a wide range of subsurface imaging and inversion experiments. In particular, recorded low-frequency (<2 Hz) surface-wave information holds the promise of providing constraints on the shear-wave velocity (VS) to depths exceeding 0.5 km. However, surface-wave analysis can be made challenging by a number of acquisition factors that affect the amplitudes of measured DAS waveforms. To illustrate these sensitivity challenges, we present a low-frequency ambient wavefield investigation using a DAS data set acquired on a crooked-line optical fibre array deployed in suburban Perth, Western Australia. We record storm-induced microseism energy generated at the nearby Indian Ocean shelf break and/or coastline in a low-frequency band (0.04−1.80 Hz) and generate high-quality virtual shot gathers (VSGs) through cross-correlation and cross-coherence interferometric analyses. The resulting VSG volumes clearly exhibit surface wave energy, though with significant along-line amplitude variations that are due to the combined effects of ambient source directivity, crooked-line acquisition geometry and the applied gauge length, fibre coupling, among other factors. We transform the observed VSGs into dispersion images using two different methods: phase shift and high-resolution linear Radon transform. These dispersion images are then used to estimate 1-D near-surface VS models using multichannel analysis of surface waves (MASW), which involves picking and inverting the estimated Rayleigh-wave dispersion curves using the particle-swarm optimization global optimization algorithm. The MASW inversion results, combined with nearby deep borehole information and 2-D elastic finite-difference modeling, show that low-frequency ambient DAS data constrain the VS model, including a low-velocity channel, to at least 0.5 km depth. Thus, this case study illustrates the potential of using DAS technology as a tool for undertaking large-scale surface wave analysis in urban geophysical and geotechnical investigations to depths exceeding 0.5 km.

Anisotropic elastic finite-difference modeling of sources and receivers on Lebedev grids

The Lebedev grid finite-difference (FD) method allows modeling of anisotropic elastic-wave propagation. On Lebedev grids, erroneous point-source excitations can create spurious (nonphysical) waves. The only known remedy for such artifacts in the literature is the Lisitsa-Vishnevsky method. This method uses a distributed array to create point sources and point receivers on the FD grid. However, the Lisitsa-Vishnevsky method does not fully eliminate spurious artifacts. A novel approach is found in the FD-consistent point source, which suppresses spurious artifacts entirely. The method requires no array recording to create point receivers. The advantage of this method over the Lisitsa-Vishnevsky method is determined with two anisotropic modeling examples.

Anisotropic elastic finite-difference modeling of sources and receivers on Lebedev grids

Anisotropic elastic finite-difference modeling of sources and receivers on Lebedev grids

Cross-profile seismic data acquisition, imaging, and modeling of iron-oxide deposits: A case study from Blötberget, south-central Sweden

Two 2D reflection seismic profiles were acquired in Blötberget, south-central Sweden, for deep targeting and delineation of sheet-like iron-oxide deposits, known to dip toward the southeast and extend down to at least 0.8 km depth from core drilling observations. The two perpendicular profiles recorded shots at every receiver station along the main and cross profiles. To obtain more information on the lateral extent of the mineralized horizons, data from the two profiles, including the cross-profile records, were binned together in a 3D grid and further processed as a 3D data set. Processing results suggest that more information is retrieved when 3D processing is used and the mineralization lateral extent can be inferred for at least 0.3 km. The seismic response of the mineralization was further studied through forward reflection traveltime modeling, using a 3D ray-tracing approach; thus, the 3D geometry of several planar reflectors was validated. Additionally, 2D elastic finite-difference modeling work showed that the observed reflection pattern in the seismic data may originate from several mineralized horizons, suggesting potential resources in the footwall of the known deposits and large-scale geologic structures. The results encourage the use of seismic methods for direct delineation of mineral deposits even from 2D profiles and prompted a 3D survey in the area.

Acoustic and elastic finite-difference modeling by optimal variable-length spatial operators

Time-space domain finite-difference modeling has always had the problem of spatial and temporal dispersion. High-order finite-difference methods are commonly used to suppress spatial dispersion. Recently developed time-dispersion transforms can effectively eliminate temporal dispersion from seismograms produced by the conventional modeling of high-order spatial and second-order temporal finite differences. To improve the efficiency of the conventional modeling, I have developed optimal variable-length spatial finite differences to efficiently compute spatial derivatives involved in acoustic and elastic wave equations. First, considering that temporal dispersion can be removed, I prove that minimizing the relative error of the phase velocity can be approximately implemented by minimizing that of the spatial dispersion. Considering that the latter minimization depends on the wavelength that is dependent on the velocity, in this sense, this minimization is indirectly related to the velocity, and thus leads to variation of the spatial finite-difference operator with velocity for a heterogeneous model. Second, I use the Remez exchange algorithm to obtain finite-difference coefficients with the lowest spatial dispersion error over the largest possible wavenumber range. Then, dispersion analysis indicates the validity of the approximation and the algorithm. Finally, I use modeling examples to determine that the optimal variable-length spatial finite difference can greatly increase the modeling efficiency, compared to the conventional fixed-length one. Stability analysis and modeling experiments also indicate that the variable-length finite difference can adopt a larger time step to perform stable modeling than the fixed-length one for inhomogeneous models.

Two-grid full-waveform Rayleigh-wave inversion via a genetic algorithm — Part 2: Application to two actual data sets

We have applied our two-grid genetic-algorithm Rayleigh-wave full-waveform inversion (FWI) to two actual data sets acquired in Luni (Italy) and Grenoble (France), respectively. Because our technique used 2D elastic finite-difference modeling for solving the forward problem, the observed data were 3D to 2D corrected prior to the inversion. To limit the computing time, both inversions focused on predicting low-resolution, smooth models by using quite coarse inversion grids. The wavelets for FWI were estimated directly from the observed data by using the Wiener method. In the Luni case, due to the strong dispersion effects on the data, to strengthen the inversion, envelopes and waveforms were considered in the objective function and an offset-marching strategy was applied. Though no a priori information was exploited, the outcomes of the Luni and Grenoble data inversion were fair. The predicted Luni [Formula: see text] model indicates a strong velocity increase from approximately 3 to 6 m, and velocity inversions have been detected at approximately 2 and 9 m depths. Analyzing the dispersion spectra, it results that the predicted Luni data reasonably reproduced the waveforms related to the fundamental mode and, likely, a small part of those related to the first higher mode. Concerning the Grenoble example, the predicted [Formula: see text] model coincides reasonably well with the long-wavelength structures presented in the [Formula: see text] profiles obtained from nearby boreholes. The data reconstruction is generally satisfactory, and when mismatches occur between the predicted and observed traces, the phase differences are always within half-periods. The fair inversion outcomes suggest that the predicted Luni and Grenoble models would likely be adequate initial models for local FWI, which could further increase the resolution and the details of the estimated [Formula: see text] models.

Sub-basalt Imaging of Hydrocarbon-Bearing Mesozoic Sediments Using Ray-Trace Inversion of First-Arrival Seismic Data and Elastic Finite-Difference Full-Wave Modeling Along Sinor\u2013Valod Profile of Deccan Syneclise, India

Imaging below the basalt for hydrocarbon exploration is a global problem because of poor penetration and significant loss of seismic energy due to scattering, attenuation, absorption and mode-conversion when the seismic waves encounter a highly heterogeneous and rugose basalt layer. The conventional (short offset) seismic data acquisition, processing and modeling techniques adopted by the oil industry generally fails to image hydrocarbon-bearing sub-trappean Mesozoic sediments hidden below the basalt and is considered as a serious problem for hydrocarbon exploration in the world. To overcome this difficulty of sub-basalt imaging, we have generated dense synthetic seismic data with the help of elastic finite-difference full-wave modeling using staggered-grid scheme for the model derived from ray-trace inversion using sparse wide-angle seismic data acquired along Sinor–Valod profile in the Deccan Volcanic Province of India. The full-wave synthetic seismic data generated have been processed and imaged using conventional seismic data processing technique with Kirchhoff pre-stack time and depth migrations. The seismic image obtained correlates with all the structural features of the model obtained through ray-trace inversion of wide-angle seismic data, validating the effectiveness of robust elastic finite-difference full-wave modeling approach for imaging below thick basalts. Using the full-wave modeling also allows us to decipher small-scale heterogeneities imposed in the model as a measure of the rugose basalt interfaces, which could not be dealt with ray-trace inversion. Furthermore, we were able to accurately image thin low-velocity hydrocarbon-bearing Mesozoic sediments sandwiched between and hidden below two thick sequences of high-velocity basalt layers lying above the basement.

Full elastic finite-difference modeling and interpretation of karst system in a subsalt carbonate reservoir

We evaluated an advanced forward-modeling-based reservoir characterization technique that uses full elastic finite-difference simulation to investigate the limits of karst identification in stacked seismic data. Identification of karsts is important for field development in carbonated reservoirs because paleokarst features can result in a loss of circulation and/or sometimes lost drill bits. Our primary objective was to verify whether we can detect and interpret the location and size of karsts from seismic data, especially given a complex overburden. We constructed an elastic reservoir model consisting of compressional velocity ([Formula: see text]), shear velocity ([Formula: see text]), and density for the study area using interpreted horizons and well log information. Karsts with varying widths, thicknesses, dip angles, and porosities were inserted to generate multiple versions of the model. We also evaluated the imaging impact of overlying faults and salt on karst detection. Full elastic simulation was performed on the various reservoir models using a realistic acquisition geometry to generate gathers, which were then prestack time migrated to quantify the impact of different karst properties on the seismic images and study the effect of reservoir property changes on the seismic response. Finally, a wave-equation target-oriented analysis was presented to improve the understanding of subsalt amplitude and illumination. From the finite-difference modeling and analysis that we performed, we obtained an uncertainty range on karst property estimation from seismic images and gained insights into future survey design for subsalt interpretation and amplitude analysis. For our specific model, we found the limit of karst identification from seismic data is a 30-m-wide horizontal karst or a 500 m karst dipping at 60°. Also, the karst image width reflected its true width only when the actual karst width was larger than the P-wave wavelength (240 m in this case). With a dipping overburden above the reservoir, apparent positions of karsts were shifted in the updip direction by prestack time migration up to 50 m from their true position. This lateral uncertainty should be kept in mind in well planning to avoid karst features interpreted from a time-migrated seismic section.

Scaling behavior and the effects of heterogeneity on shallow seismic imaging of mineral deposits: A case study from Brunswick No. 6 mining area, Canada

Abstract We have studied the scaling behavior of compressional-wave velocity and density logs from an exploration borehole that extends down to about 700 m depth in the Brunswick No. 6 mining area, Bathurst Mining Camp, Canada. Using statistical methods, vertical and horizontal scale lengths of heterogeneity were estimated. Vertical scale length estimates from the velocity, density and calculated acoustic impedance are 14 m, 33 m, and about 20 m, respectively. Although the estimated scale length for the acoustic impedance implies a weak scattering environment, elastic finite difference modeling of seismic wave propagation in 2D heterogeneous media demonstrates that even this weak scattering medium can mask seismic signals from small, but yet economically feasible, massive sulfide deposits. Further analysis of the synthetic seismic data suggests that in the presence of heterogeneity, lenticular-shaped targets may only exhibit incomplete diffraction signals whereby the down-dip tails of these diffractions are mainly visible on the stacked sections. Therefore, identification of orebody generated diffractions is much easier on the unmigrated stacked sections than on migrated stacked sections. The numerical seismic modeling in 2D heterogeneous media indicates that in the presence of large horizontal, but small vertical scale lengths (structural anisotropy), identification of massive sulfide deposits is possible, but their delineation at depth requires detailed velocity modeling and processing algorithms which can handle the anisotropy.

Elastic finite-difference modeling of volcanic-hosted massive sulfide deposits: A case study from Half Mile Lake, New Brunswick, Canada

We present elastic finite-difference modeling results over a geologically realistic 2D representation of the Half Mile Lake volcanic-hosted massive sulfide deposit, New Brunswick, Canada. The model is constrained by geologic information from surface mapping and boreholes, whereas petrophysical properties are provided by wireline logging data acquired in two boreholes intersecting different parts of the deposit. We analyzed the P-P, P-S, S-P, and S-S responses of the lower and deep mineralized zones and assessed some compositional effects by substituting massive sulfides with gabbro properties in the model. Finite-difference modeling results predict complex scattering signature associated with the lower and deep sulfide zones. Both zones scattered back P-P, P-S, S-P, and S-S waves generally having strongest amplitudes in the stratigraphy down-dip direction. The P-S, S-P, and S-S scattered waves, if properly recorded on multicomponent data, represent useful signal that could help the targeting of deep sulfide mineralization. Finite-difference simulations further reveal phase-reversals on P-P wavefields scattered at the lower and deep zones. The phase reversals are not observed for gabbro inclusions, suggesting that this signature could be used to discriminate gabbro units from sulfide mineralization. The finite-difference simulation successfully reproduces many events of the VSP data, in particular P-S and S-S events on the radial component and P-P and S-P events on the vertical component. Comparison with 3D data is rather poor and only shows weak correlation with P-P events from the lower and deep zones. Despite the poor correlation, a prestack time migrated S-P section displays an amplitude anomaly at the location of the deep zone, suggesting that S-P waves were recorded on the 3D data, although this survey was acquired with explosive sources and vertical geophones.

Rayleigh-to-shear wave conversion at the tunnel face — From 3D-FD modeling to ahead-of-drill exploration

For safe tunnel excavation, it is important to predict lithologic and structural heterogeneities ahead of construction. Conventional tunnel seismic prediction systems utilize body waves (P- and S-waves) that are directly generated at the tunnel walls or near the cutter head of the tunnel boring machine (TBM). We propose a new prediction strategy that has been discovered by 3D elastic finite-difference (FD) modeling: Rayleigh waves arriving at the front face of the tunnel are converted into high-amplitude S-waves propagating further ahead. Reflected or backscattered S-waves are converted back into Rayleigh waves which can be recorded along the sidewalls. We name these waves RSSR waves. In our approach, the front face acts as an S-wave transceiver. One technical advantage is that both the sources and the receivers may be placed behind the cutter head of the TBM. The modeling reveals that the RSSR waves exhibit significantly higher amplitudes than the directly reflected body waves. The excavation damage zone causes dispersion of the RSSR wave leading to multimodal reflection response. For the detection of geologic interfaces ahead, RSSR waves recorded along the sidewalls are corrected for dispersion and stacked. From the arrival times, the distance to the S-S reflection point can be estimated. A recurrent application, while the tunnel approaches the interface, allows one to quantify the orientation of the reflecting interfaces as well. Our approach has been verified successfully in a field experiment at the Piora adit of the Gotthard base tunnel. The distance to the Piora fault zone estimated from stacked RSSR events agrees well with the information obtained by geologic surveying and exploratory drilling.

Origin of upper-mantle seismic scattering - evidence from Russian peaceful nuclear explosion data

SUMMARY Detailed interpretation of upper-mantle scattered waves recorded in the 800–1400 km offset range of four reversed ‘peaceful nuclear explosion’ seismic sections along the 3500 km long profile Kraton in Russia has revealed an ∼85 km thick seismically inhomogeneous low-velocity zone below the 8° discontinuity at ∼100 km depth. The scattered waves make up a wave train of high amplitude that decreases in duration from ∼8 s at 800 km offset to ∼4 s at 1400 km offset. 2-D elastic finite-difference modelling of the seismic wavefield has previously shown that the inhomogeneous zone from ∼100 to ∼185 km depth may be described by continuous, random velocity fluctuations represented by an exponential medium with a horizontal correlation length of 5–10 km, a vertical correlation length of <5 km and a standard deviation of ∼2 per cent of the average background velocity value. Here, we demonstrate that other proposed velocity fluctuations in the lower crust and in the uppermost mantle from the Moho to ∼100 km depth cannot explain the observed scattered phases. Our study shows that upper-mantle velocity fluctuations below the 8° discontinuity at 100 km depth are needed in order to explain the observed seismic wavefield, and that they lead to scattering that is significantly different from the scattering caused by the fluctuations at shallower levels.

Imaging of groundwater resources in glacial deposits using high-resolution reflection seismics, Sweden

Two high-resolution reflection seismic profiles were acquired in the Heby area of eastern Sweden over glacial deposits for the purpose of mapping groundwater resources. The majority of shot points were located in clay resulting in good quality data along most of the profiles. On stacked and migrated sections, the uppermost clay is about 20 m thick and is characterized by its subhorizontal reflectivity. Sand/gravel deposits below it contain more dipping interfaces and have a chaotic reflectivity pattern. Depth to bedrock is interpreted to be 90 and 65 m on the respective profiles and occurs in about a 100-m-wide trough on both profiles. Reflections from the tops of sandy gravel zones generally have higher amplitudes. Clear reflections from a thin silt layer (20 cm thick) at about 10-m depth are observed on one of the profiles. Elastic finite difference modeling and the observation of this reflection in shot gathers show that the reflection is not an artifact of the acquisition nor the processing. The modeling also shows that there is no marked low-velocity waveguide in the near surface, but that an effective low Q zone may be present. Comparison with refraction profiling on the other profile shows that there is better agreement between the reflection seismic results and penetration tests than the refraction results with these tests. Both profiles allow the thickness of the overlying clay layers to be determined, as well as the thickness of the underlying sand/gravel deposits. This is important for estimating the amount of groundwater resources in an area.

3D Elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing

AbstractThis article provides a technique to model seismic motions in 3D elastic media using fourth-order staggered-grid finite-difference (FD) operators implemented on a mesh with nonuniform grid spacing. The accuracy of the proposed technique has been tested through comparisons with analytical solutions, conventional 3D staggered-grid FD with uniform grid spacing, and reflectivity methods for a variety of velocity models. Numerical tests with nonuniform grids suggest that the method allows sufficiently accurate modeling when the grid sampling rate is at least 6 grid points per shortest shear wavelength. The applicability for a finite fault with non-uniform distribution of point sources is also confirmed. The use of nonuniform spacing improves the efficiency of the FD methods when applied to large-scale structures by partially avoiding the spatial oversampling introduced by the uniform spacing in zones with high velocity. The significant reduction in computer memory that can be obtained by the new technique improves the efficiency of the 3D-FD method at handling shorter wavelengths, larger areas, or more realistic 3D velocity structures.

Elastic finite-difference modeling of strong motion in Ashigara Valley for the 1990 Odawara, Japan, earthquake

Abstract Finite-difference seismograms are calculated for the 1990 Odawara earthquake to model transverse and radial components of strong-motion data. Using a basin structure derived from a geotechnical survey and double-couple solution determined from P first-motion data and SH/SV amplitude ratio, two-dimensional SH and P/SV finite-difference calculation, along with correction for soft surface layers, reproduces the observed spatial variation in amplitude and waveform of the transverse and radial records at the Ashigara Valley blind prediction sites for frequencies up to 2 Hz. Comparison of the synthetic and the observed seismograms suggests that the observed spatial variation in the SH amplitude is attributed to the source radiation and amplification due to the soft surface layers, while the SV amplitude variation is controlled by the amplification due to the sedimentary basin as well as the source and surface-layer effects.

3D eight-order elastic finite-difference modeling of refraction and strong-motion data from the Coyote Lake region, California

Abstract 3D elastic modeling of (vertical component, explosive source) refraction data is performed to constrain the shallow part of a velocity model, followed by 3D modeling of three-component strong-motion data from a (spatially and temporally extended) earthquake source in that model. Modeling is by eighth-order finite differencing on a staggered grid, which automatically generates all body and surface waves, including multiples and conversions. Application is to the 1981 USGS refraction survey and the 6 August 1979 Coyote Lake event in central California. 3D modeling was able to reproduce the main observed features in both the refraction and earthquake data. For the refraction data, the synthetics explain local variations in first-arrival travel times, reverberations in the near-surface sediments, and surface waves. For the strong-motion data, the initial arrival times and amplitudes are determined by the position and orientation of the initial rupture; the shapes of the trailing edges of the pulses are modified by the actual source distribution. Changes of source position of 1 to 2 km are easily resolvable through changes in the position of nodal planes. Discontinuities in material properties (e.g., layering) produce a coda of multiples and converted waves. Lateral variations in velocities, especially in the upper few kilometers, change the energy distribution among the trace components, by 3D refraction, and so shift the apparent source directivity and position, as well as contributing to the coda. Remaining differences between the observed and synthetic data are attributed to local seismometer coupling effects, very local site structures, and phenomena such as anelasticity and anisotropy, which are not explicitly modeled.

Nonlinear elastic waveform inversion of land seismic reflection data

Data from 53 shot records collected in a conventional, land, seismic reflection survey are inverted to obtain subsurface maps of short‐wavelength variations in elastic parameters (impedance). The inversion method employs a robust, nonlinear, iterative, gradient algorithm which inverts waveform data; thus there is no picking of travel times or amplitudes. The algorithm is also computationnally expensive as it relies heavily upon elastic finite difference modeling. In addition, because a gradient search is employed to minimize the misfit between observed and modeled data, we attempt to solve only the weakly nonlinear problem associated with the short wavelenghts of imedance (information on the short‐wavelengths of impedance is mainly derived from amplitude variations of the data). Resolving the long wavelengths of velocity (background) is thought to be a highly nonlinear problem, at least for waveform inversion, and is not attempted herein. For the inversion of this data set, long wavelenghts have been derived from well logs located along the seismic reflection profile. Although the data set is dominated by unconverted P wave arrivals, we obtain maps of the short wavelenghts of both P impedance and S impedance. Three well logs in the survey area confirm the validity of the P impedance map; however, no S velocity well log are available from this area. There are channel sands of economic interest which have a thickness of approximately 10 m in one well and are very thin or absent in the other two wells. The 10 m occurrence of channel sands is quite visible in both impedance maps. We note that because inversion in three dimensions of a realistically sized data set is simply not possible on present‐day computers, all modeling was performed in two dimensions, and a rather crude conversion from spherical to cylindrical energy spreading was applied to the original field data before inversion. Other seismic reflection data in the survey area suggest that there is little lateral variation in the subsurface except for the channel sands; thus we feel that a two‐dimensional inversion with energy compensation is justified.

Elastic finite difference modelling of lower crustal reflections

It is a common feature to observe on deep reflection profiles, short (less than the Fresnel area diameter) and strongly discontinuous lower crustal reflections. Such wave patterns can be attributed to variations in the crustal properties on the scale of a few wavelengths. In order to explain these characteristics, we have computed finite-difference synthetic seismograms for various types of deep heterogeneities. Focusing, defocusing and interference effects were especially studied. The results can only be interpreted in a qualitative manner because our models exhibit departures from the real cases (lower frequency content, and only two-dimensional structures). It is found that simple undulations of the deep interfaces can hardly produce the discontinuous pattern observed, but high velocity intrusions generate complex interferences matching the reflection character of real data. Furthermore, the reflection patterns strongly depend both on length and shape of the inclusions.