Thomas Bohlen

Paralel 3-D viscoelastic finite difference seismic modelling

Computational power has advanced to a state where we can begin to perform wavefield simulations for realistic (complex) 3-D earth models at frequencies of interest to both seismologists and engineers. On serial platforms, however, 3-D calculations are still limited to small grid sizes and short seismic wave traveltimes. To make use of the efficiency of network computers a parallel 3-D viscoelastic finite difference (FD) code is implemented which allows to distribute the work on several PCs or workstations connected via standard ethernet in an in-house network. By using the portable message passing interface standard (MPI) for the communication between processors, running times can be reduced and grid sizes can be increased significantly. Furthermore, the code shows good performance on massive parallel supercomputers which makes the computation of very large grids feasible. This implementation greatly expands the applicability of the 3-D elastic/viscoelastic finite-difference modelling technique by providing an efficient, portable and practical C-program.

Finite‐difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid

We describe the application of the rotated staggered-grid (RSG) finite-difference technique to the wave equations for anisotropic and viscoelastic media. The RSG uses rotated finite-difference operators, leading to a distribution of modeling parameters in an elementary cell where all components of one physical property are located only at one single position. This can be advantageous for modeling wave propagation in anisotropic media or complex media, including high-contrast discontinuities, because no averaging of elastic moduli is needed. The RSG can be applied both to displacement-stress and to velocity-stress finite-difference (FD) schemes, whereby the latter are commonly used to model viscoelastic wave propagation. With a von Neumann-style anlysis, we estimate the dispersion error of the RSG scheme in general anisotropic media. In three different simulation examples, all based on previously published problems, we demonstrate the application and the accuracy of the proposed numerical approach.

Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves

Heterogeneous finite-difference (FD) modeling assumes that the boundary conditions of the elastic wavefield between material discontinuities are implicitly fulfilled by the distribution of the elastic parameters on the numerical grid. It is widely applied to weak elastic contrasts between geologic formations inside the earth. We test the accuracy at the free surface of the earth. The accuracy for modeling Rayleigh waves using the conventional standard staggered-grid (SSG) and the rotated staggered grid (RSG) is investigated. The accuracy tests reveal that one cannot rely on conventional numerical dispersion discretization criteria. A higher sampling is necessary to obtain acceptable accuracy. In the case of planar free surfaces aligned with the grid, 15 to 30 grid points per minimum wavelength of the Rayleigh wave are required. The widely used explicit boundary condition, the so-called image method, produces similar accuracy and requires approximately half the sampling of the wavefield compared to heterogeneous free-surface modeling. For a free-surface not aligned with the grid (surface topography), the error increases significantly and varies with the dip angle of the interface. For an irregular interface, the RSG scheme is more accurate than the SSG scheme. The RSG scheme, however, requires 60 grid points per minimum wavelength to achieve good accuracy for all dip angles. The high computation requirements for 3D simulations on such fine grids limit the application of heterogenous modeling in the presence of complex surface topography.

1.5D inversion of lateral variation of Scholte-wave dispersion

Reliable models of in-situ shear-wave velocities of shallow-water marine sediments are important for geotechnical applications, lithological sediment characterization, and seismic exploration studies. We infer the 2D shear-wave velocity structure of shallow-water marine sediments from the lateral variation of Scholte-wave dispersion. Scholte waves are recorded in a common receiver gather generated by an air gun towed behind a ship away from a single stationary ocean-bottom seismometer. An offset window moves along the common receiver gather to pick up a local wavefield. A slant stack produces a slowness-frequency spectrum of the local wavefield, which contains all modes excited by the air gun. Amplitude maxima (dispersion curves) in the local spectrum are picked and inverted for the shear-wave velocity depth profile located at the center of the window. As the window continuously moves along the common receiver gather, a 2D shear-wave velocity section is generated. In a synthetic example the smooth lateral variation of surficial shear-wave velocity is well reconstructed. The method is applied to two orthogonal common receiver gathers acquired in the Baltic Sea (northern Germany). The inverted 2D models show a strong vertical gradient of shear-wave velocity at the sea floor. Along one profile significant lateral variation near the sea floor is observed.

Shear wave anisotropy of laminated lower crust beneath Urach (SW Germany): a comparison with xenoliths and with exposed lower crustal sections

In order to investigate the laminated lower crust beneath Urach (SW Germany), we analysed (1) the anisotropy and reflectivity of seismic shear waves, and (2) the composition of xenoliths samples. We conclude that the directional dependence of shear waves of the lower crust is effectively transversely isotropic with a maximum velocity difference of 3 to 6%. Both seismic signature and petrological information are compared with reference profiles deduced from the exposed lower crustal sections of Ivrea and Calabria (Italy). The comparison, based on synthetic seismograms, shows that the observed combination of shear wave anisotropy and reflectivity pattern corresponds to alternating layers containing a high amount of metapelites, such as found in the Calabria profile and in the Val Strona subsection of Ivrea. The main portion of the whole rock anisotropy is caused by the lattice-preferred orientation of anisotropic minerals (so-called intrinsic anisotropy). The effect of alternating rock layering (so-called layering anisotropy) contributes only a minor portion to the whole rock anisotropy. Layered mafic intrusions, such as found in the Val Sesia subsection of Ivrea, are nearly isotropic. Both seismic anisotropy and xenoliths indicate that mafic layers do not play a major role for the laminated lower crust at Urach. They may be restricted to a thin layer at the crust-mantle boundary.

Superdeep vertical seismic profiling at the KTB deep drill hole (Germany): Seismic close-up view of a major thrust zone down to 8.5 km depth

The lowermost section of the continental superdeep drill hole German Continental Deep Drilling Program (KTB) ( south Germany) has been investigated for the first time by vertical seismic profiling (VSP). The new VSP samples the still accessible range of 6 - 8.5 km depth. Between 7 and 8.5 km depth, the drill hole intersects a major cataclastic fault zone which can be traced back to the Earths surface where it forms a lineament of regional importance, the Franconian line. To determine the seismic properties of the crust in situ, in particular within and around this deep fault zone, was one of the major goals of the VSP. For the measurements a newly developed high-pressure/high-temperature borehole geophone was used that was capable of withstanding temperatures and pressures up to 260 degreesC and 140 MPa, respectively. The velocity-depth profiles and reflection images resulting from the VSP are of high spatial resolution due to a small geophone spacing of 12.5 m and a broad seismic signal spectrum. Compared to the upper part of the borehole, we found more than 10% decrease of the P wave velocity in the deep, fractured metamorphic rock formations. P wave velocity is similar to 5.5 km/s at 8.5 km depth compared to 6.0 - 6.5 km/s at more shallow levels above 7 km. In addition, seismic anisotropy was observed to increase significantly within the deep fracture zone showing more than 10% shear wave splitting and azimuthal variation of S wave polarization. In order to quantify the effect of fractures on the seismic velocity in situ we compared lithologically identical rock units at shallow and large depths: Combining seismic velocity and structural logs, we could determine the elastic tensors for three gneiss sections. The analysis of these tensors showed that we need fracture porosity in the percent range in order to explain seismic velocity and anisotropy observed within the fault zone. The opening of significant pore space around 8 km depth can only be maintained by differential tectonic stress combined with intense macroscopic fracturing. VSP reflection imaging based on PP and PS converted reflected waves showed that the major fault system at the KTB site is wider and more complex than previously known. The so-called SE1 reflection previously found in two- and three-dimensional surface seismic surveys corresponds to the top of an similar to1 km wide fault system. Its lower portion was not illuminated by surface seismic acquisition geometry. VSP imaging shows that the fault zone comprises two major and a number of smaller SE dipping fault planes and several conjugate fracture planes. The previously recognized upper fault plane is not associated with a strong velocity anomaly but indicates the depth below which the dramatic velocity decrease starts. Regarding the complexly faulted crustal section of the KTB site as a whole, we found that fluctuation spectra of rock composition and seismic velocity show similar patterns. We could verify that a significant amount of P wave energy is continuously converted into shear energy by forward scattering and that multipathing plays an important role in signal formation. The media behaves effectively smoothly only at wavelength larger than 150 m. It was shown by moving source profiling that the media is orthorhombic on are regional scale. The tilt of the symmetry axes of anisotropy varies with depth following the dip of the geological structure.

3D elastic full-waveform inversion of small-scale heterogeneities in transmission geometry

Three-dimensional elastic full-waveform inversion aims to reconstruct elastic material properties of 3D structures in the subsurface with high resolution. Here we present an implementation of 3D elastic full-waveform inversion based on the adjoint-state method. The code is optimized regarding runtime and storage costs by using a time-frequency approach. The gradient is computed from monochromatic frequency-domain particle-velocity wavefields calculated with a time-domain velocity-stress finite-difference scheme. The 3D full-waveform inversion was applied to data of a complex random medium model, which resembles a realistic crystalline rock environment. We show synthetic inversion results of P-wave and S-wave velocities for two transmission geometries: (1) a 3D acquisition geometry with planes of sources and receivers and (2) a 2D geometry with two lines of sources and receivers, resembling a realistic two-borehole geometry. The 3D inversion of data acquired with 3D source-receiver geometry is capable to reconstruct differently sized 3D structures of shear and compressional velocities with resolution of about a wavelength. The 3D random medium data recorded with 2D acquisition geometry were inverted using 3D inversion and 2D full-waveform inversion for comparison. The 2D inversion suffers from strong artefacts that are caused by 3D scattering. The multiparameter 3D inversion, by contrast, is capable to invert the 3D scattered waves and to reconstruct 3D structures up to about 1–2 wavelengths adjacent to the plane between sources and receivers. The resolution is lower compared to the 3D acquisition geometry result. Still, a 3D inversion of cross-hole data can be beneficial compared to a 2D inversion in the presence of complex 3D small-scale heterogeneities, as it is capable to resolve 3D structures next to the source-receiver plane.

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.

Reverse Time Migration: A Seismic Imaging Technique Applied to Synthetic Ultrasonic Data

Ultrasonic echo testing is a more and more frequently used technique in civil engineering to investigate concrete building elements, to measure thickness as well as to locate and characterise built-in components or inhomogeneities. Currently the Synthetic Aperture Focusing Technique (SAFT), which is closely related to Kirchhoff migration, is used in most cases for imaging. However, this method is known to have difficulties to image steeply dipping interfaces as well as lower boundaries of tubes, voids or similar objects. We have transferred a processing technique from geophysics, the Reverse Time Migration (RTM) method, to improve the imaging of complicated geometries. By using the information from wide angle reflections as well as from multiple events there are fewer limitations compared to SAFT. As a drawback the required computing power is significantly higher compared to the techniques currently used. Synthetic experiments have been performed on polyamide and concrete specimens to show the improvements compared to SAFT. We have been able to image vertical interfaces of step-like structures as well as the lower boundaries of circular objects. It has been shown that RTM is a step forward for ultrasonic testing in civil engineering.

Variability of Scholte-wave Dispersion in Shallow-water Marine Sediments

Models of in situ shear-wave velocities of shallow-water marine sediments are of importance for geotechnical applications, sediment characterization, and seismic exploration studies. Here pseudo-2D shear-wave velocity models are inferred from the lateral variation of Scholte-wave dispersion at five different geological sites in the Baltic Sea (northern Germany). To explore Scholte-wave dispersion and the lateral variability of shear-wave velocities, Scholte waves were excited by air gun shots in the water layer and recorded by stationary ocean-bottom-seismometers or buried geophones. We analyze the recorded seismograms in a common-receiver-gather using offset-windowed, multichannel dispersion analysis. The observed local slowness-frequency spectra for the different study sites vary significantly with respect to excitation amplitudes and phase slownesses of different modes, as well as the excited frequency range. The excitation amplitudes are influenced by the local shear-wave velocity structure, absorption, length of Scholte-wave travel path, and the elevation of the source above the sea floor. The inverted shear-wave velocities range from 50 m/s to 600 m/s. Directly at the sea bottom, shear-wave velocities of 50 m/s for fine muddy sand and 300 m/s for glacial till were inferred. The maximum vertical gradient was +/- 80% (mean 250 m/s) within a depth range of 40 m, and horizontally +/- 33% (mean 350 m/s) within 300 m distance. The layer boundaries in the inverted shear-wave velocity models are in good agreement with high-frequency, zero-offset compressional-wave reflections. However, it was not possible to acquire the fundamental Scholte mode above very soft, unconsolidated sediment with shear-wave velocities smaller 50 m/s. The analysis of synthetic data shows that this is due to the elevation of the source and the receiver response function.