Numerical modeling of zero-offset laboratory data in a strong topographic environment: results for a spectral-element method and a discretized Kirchhoff integral method

Abstract

Accurate simulation of seismic wave propagation in complex geological structures is of particular interest nowadays. However conventional methods may fail to simulate realistic wavefields in environments with great and rapid structural changes, due for instance to the presence of shadow zones, diffractions and/or edge effects. Different methods, developed to improve seismic modeling, are typically tested on synthetic configurations against analytical solutions for simple canonical problems or reference methods, or via direct comparison with real data acquired in situ. Such approaches have limitations, especially if the propagation occurs in a complex environment with strong-contrast reflectors and surface irregularities, as it can be difficult to determine the method which gives the best approximation of the “real” solution, or to interpret the results obtained without an a priori knowledge of the geologic environment. An alternative approach for seismics consists in comparing the synthetic data with high-quality data collected in laboratory experiments under controlled conditions for a known configuration. In contrast with numerical experiments, laboratory data possess many of the characteristics of field data, as real waves propagate through models with no numerical approximations. We thus present a comparison of laboratory-scaled measurements of 3D zero-offset wave reflection of broadband pulses from a strong topographic environment immersed in a water tank with numerical data simulated by means of a spectral-element method and a discretized Kirchhoff integral method. The results indicate a good quantitative fit in terms of time arrivals and acceptable fit in amplitudes for all datasets.

1 Introduction

Accurate simulation of seismic wave propagation in complex geological structures is widely used for environmental and industrial applications for subsurface structure evaluation and in seismic exploration as a core tool of seismic imaging and inversion methods (Robertsson et al. 2007; Virieux et al. 2011). In models with simple structures and slowly varying material properties, conventional methods (e.g., ray methods, finite-difference methods) are efficient. However, difficulties arise for complex geological structures with large and rapid structural changes, and conventional methods may then fail to simulate realistic wavefields due to the presence of shadow zones, diffractions and/or edge effects. Different methods have thus been developed to improve seismic modeling in complex geological environments with structural complexities like faults with steep dips or curved reflectors. In the literature, they are typically tested on synthetic configurations against analytical solutions for simple canonical problems or reference methods, and several projects focusing on verification and validation of numerical methods have been conducted in the last few years (Igel et al. 2000; Moczo et al. 2006; Chaljub et al. 2010; Fehler and Keliher 2011). Such an approach has limitations, especially if the propagation occurs in a complex environment with strong-contrast reflectors and surface irregularities, as it can be difficult to determine the method which gives the best approximation of the “real” solution given by a reference method. Another approach is to validate these methods via direct comparison with real data acquired in situ (Houbiers et al. 2012). Unfortunately, without a priori (good) knowledge of the geological environment, the interpretation of the obtained results may be a tedious task due to the existence of diffraction and sideswipe events.

An alternative approach for seismics consists in comparing the synthetic data with data obtained in the laboratory. The main assumption underlying this approach is that the scaled physical mechanisms are identical to the field physical mechanisms (Ebrom and McDonald 1994). Such an approach holds in our case because we will resort to a linear wave equation. In contrast with in situ experiments, high-quality data are collected under controlled conditions for a known configuration, which is crucial for comparisons with numerical models. Moreover, unlike synthetic data, laboratory data possess many of the characteristics of field data (random and signal-generated noise, multiples, mode conversions), as real waves propagate through models with no numerical approximations. Since the mid part of the twentieth century small-scaled physical modeling approaches have been extensively exploited in seismics for a better understanding of wave propagation phenomena (Oliver et al. 1954; Angona 1960), for validation of theoretical predictions (Howson and Sinha 1984; Pant et al. 1992), or for testing seismic processing and imaging techniques (French 1974; Macdonald et al. 1987; Lo et al. 1988). Nowadays, physical modeling is still used for complex 3D configurations whose response is difficult to model numerically, such as anisotropic (Stewart et al. 2011), fractured (Ekanem et al. 2013), or anelastic (Lines et al. 2012) media. Data from laboratory experiments are also employed to test new data processing algorithms (Campman et al. 2005), in time-lapse 3D studies (Sherlock et al. 2000), and to benchmark numerical model solutions (Bretaudeau et al. 2011).

Our aim is to study 3D complications in zero-offset reflection profiles acquired over a strong topographic environment in order to improve the understanding of the physical mechanisms involved in the interaction of the waves with irregular surfaces. As noted previously, in such a complex environment the numerical methods based on approximations may fail to simulate accurately the seismic wavefields and produce different results depending on their intrinsic hypotheses. The main purpose of this article is, therefore, to test the approach using laboratory data as reference data for benchmarking 2D and 3D numerical methods and techniques. Using the indoor tank facilities of the Laboratory of Mechanics and Acoustics (LMA, Marseille, France), we have performed laboratory-scaled measurements of zero-offset reflection of broadband pulses on a model containing topographic structures with several edges and corners and immersed in a water tank. The presence of these structures is expected to complicate the wavefields significantly. In what follows, we will present comparisons of these measurements with numerical data simulated by means of a spectral-element method (SEM) and a discretized Kirchhoff integral method (DKIM). Here, we will also only be concerned with propagation problems; discussion of migration problems is left for further articles.

The article is organized as follows. In the first section, we describe the experimental setup, and we present the laboratory data obtained in zero-offset configurations with narrow-beam and broad-beam sources. The second section is devoted to a short description of the numerical modeling methods. In the third section, we compare the laboratory and synthetic data and discuss the results obtained.

2 Small-scaled seismic experiments

We carried out laboratory experiments at LMA in Marseille, France. The model used in these experiments, called the “Marseille model,” is partly based on French’s model (French 1974), but contains additional topographies such as a truncated dome and a truncated pyramid (Fig. 1). The model of size 600 × 400 × 70 mm³ is made of PVC material, which is isotropic at ultrasonic frequencies and whose measured properties are in the same range as those of typical geological media. The thickness of the model varies from 30 to 70 mm, the difference between the two levels separated by a planar fault being 40 mm. The model is immersed in a water tank equipped with a computer-controlled system which allows for accurate positioning of the source and receiver. The measured properties of the materials are V_P = 1,476–1,493 m/s (depending on the water temperature) and ρ = 1,000 kg/m³ in the water layer, and V_P = 2,220 m/s, V_S = 1,050 m/s, and ρ = 1,412 kg/m³ in the PVC material. Attenuation in the PVC layer is described by quality factors 40 < Q_P < 60 and 27 < Q_S < 31 for P and S-waves, respectively. Attenuation in the water is negligible. As we consider zero-offset seismic configuration only for these experiments, we illuminate the model using a piezoelectric transducer which operates both as a source and a receiver. We use different kinds of transducers with a dominant frequency f_c = 500 kHz and a large frequency bandwidth: one transducer with diameter D = 25.4 mm and narrow-beam (NB) aperture (about 8.3° at −3 dB), and one transducer with diameter D = 3 mm and broad-beam (BB) aperture (about 45° at −3 dB), which allows us to obtain 2D and 3D zero-offset data.

Fig. 1

Left The Marseille model and specific acquisition lines. Right Area covered by acquisition with the narrow-beam transducer (light gray) and the broad-beam transducer (dark gray)

We use a conventional ultrasonic pulse-echo technique to obtain the reflection data. The distance from the transducer to the flat part of the model surface is either 105 mm (±1 mm), or 150 mm (±1 mm), the far-field condition being thus fulfilled for any position of the transducer (Zemanek 1971). Wave propagation occurs in small-scale conditions, i.e., for instance, if a scale ratio of 2 × 10⁴ is considered an experimental frequency of 500 kHz corresponds to a real frequency of 25 Hz, and an experimental distance of 10 mm corresponds to a real distance of 200 m, velocities as well as densities and attenuations remaining unchanged.

We perform acquisitions along Y-lines with a spatial sampling Δx = 2 mm (Fig. 1). The collected data thus consist of numerous parallel profiles composed of a collection of reflection data for dense grids of source–receiver locations. We process reflection data to produce seismograms corresponding to different cross-sections of the model. We pay more attention to a few specific profiles because they present a high interest as they cross the main structures of the model. The data collected along these profiles might thus contain reflections and diffractions from all the structures. These profiles are represented by heavy-dashed lines in Fig. 1. The seismograms are obtained after application of a low-pass filter to raw data in order to eliminate the harmonic resonances of the transducers. Additionally, for visualization purposes, we apply a clipping procedure with different clipping numbers x to all seismograms presented below, i.e., saturation of all the signals whose amplitude is greater than x % of the maximum amplitude in order to enlighten weaker signals. Small variations in the two-way travel times for reflections from the plane parts of the model can be observed. They are due to non-constant tilt of the Marseille model during acquisition, with position errors of less than 2 %. Varying water temperature in the tank during long-lasting acquisitions can also have some effect on the time arrivals (Bilaniuk and Wong 1993).

Figure 2 shows seismograms corresponding to line Y150 obtained with the NB and BB transducers. We can note that the recorded signals are of high quality due to a high signal-to-noise ratio. Events with a time arrival smaller than 210 μs correspond to primary reflections from the top surface of the model, whereas events with greater time arrival correspond to either multiples or reflections from the bottom surface. Noticeable discrepancies between data obtained with the two kinds of transducers can be identified in Fig. 2. More specifically, diffractions at the edges of the topographic structures can be clearly observed only on data obtained with the BB transducer. The steep slope of the fault, however, remains invisible whatever the transducer used.

Fig. 2

Laboratory zero-offset data (seismograms) obtained for line Y150 with the narrow-beam transducer (left) and the broad-beam transducer (right). The transducer is located at 150 mm above the model. Clipping with x = 15 (left) and x = 40 (right) is applied

3 Numerical methods

We use two kinds of numerical methods for synthetic modeling of the zero-offset experiments: a SEM and a DKIM.

The SEM is based upon a high-order piecewise polynomial approximation of the weak formulation of the wave equation. It combines the accuracy of the pseudospectral method with the flexibility of the finite-element method (Tromp et al. 2008). In this method, the wavefield is represented in terms of high-degree Lagrange interpolants, and integrals are computed based upon Gauss–Lobatto–Legendre quadrature. This combination leading to a perfectly diagonal mass matrix leads in turn to a fully explicit time scheme which lends itself very well to numerical simulations on parallel computers. It is particularly well suited to handling complex geometries and interface matching conditions (Cristini and Komatitsch 2012). The typical element size that is required to generate an accurate mesh is of the order of λ, λ being the smallest wavelength of waves traveling in the model. Very distorted mesh elements can be accurately handled (Oliveira and Seriani 2011). The SEM may be computationally expensive depending on the size of the domain, especially for 3D domains and high-frequency simulations. We thus mesh the model in 2D with quadrangles using the open-source software package Gmsh (Geuzaine and Remacle 2009). We implement directional directivity of standard ultrasonic transducers using a set of equidistant omnidirectional sources (like a horizontal array) whose amplitude is weighted by a Hamming window. The radiation of the simulated source is directed along the vertical. It is obtained using 51 point sources distributed over a line length of 2.54 cm, which corresponds to the diameter of the NB transducer.

The DKIM is a method based on numerical evaluation of the Kirchhoff–Helmholtz surface integral, which is a powerful tool to model the scattered wavefield from a piecewise smooth interface (Tygel and Ursin 1999). Both the field and its normal derivative at the interface, appearing in the integral, are commonly computed using the Kirchhoff approximation, i.e., by multiplication of the incident wavefield and its normal derivative with the plane-wave reflection coefficient (Bleistein 1984). The discretization is performed exploiting the formulas by Ayzenberg et al. (2007). The detailed physical explanation of the method based on Huygens’ principle is presented in Ayzenberg et al. (2009).

4 Comparison of numerical simulations with experimental data

Figure 3 focuses on the wave reflection from, and transmission through, the flat and curved parts of the Marseille model, together with the diffraction effects induced by the truncated part of the small dome (line Y150), both illuminated by the NB transducer and modeled using the SEM. We can notice the cuspidal form of the interaction between the wave reflections from the (flat and curved) part of the Marseille model and the wave diffractions by the truncated part or the edge of the dome. The result of the modeling of the primary reflection from the top of the PVC material along line Y150 is shown in Fig. 4, which provides detailed comparison of two traces, corresponding to two different positions of the source above the truncated dome (i.e., just above the top of the dome and just above the truncated part of the dome). All datasets, including that illustrated in Fig. 4, are normalized by dividing the signal amplitude by the maximum amplitude of the whole dataset acquired for the same transducer, and we apply a static time-shift to all datasets. We also apply additional trace-dependent time-shifts to synthetic seismograms due to the above-mentioned data inconsistencies in travel times, i.e., we adjust the reflections from the plane parts of the model based on a normalized cross-correlation coefficient. Visual observation of the traces shows a good fit in terms of arrival times and an acceptable fit in shape and amplitude between 2D modeled and laboratory data. Nevertheless, as the simulations using the SEM are only 2D, the influence of the 3D diffraction effect produced by the truncated part of the small dome is largely overestimated (Fig. 4, right).

Fig. 3

Numerical simulations of the wave reflections, transmissions, and diffractions in the vicinity of the truncated dome using the 2D SEM. The source is located 150 mm above the flat part of the Marseille model (topleft), the truncated part of the dome (topright), the top of the dome (bottomleft), the edge of the dome (bottomright). Label R1 (respectively, R2) stands for the wave reflection from the flat (respectively, curved) surface of the Marseille model, whereas label T1 (respectively, T2) stands for the wave transmission through the flat (respectively, curved) surface. Label D stands for the wave diffraction by the truncated part or the edge of the dome

Fig. 4

Comparison of the synthetic traces (blue line) obtained with 2D simulations using the SEM and laboratory traces (green line) along line Y150, corresponding to wave reflection from the top of the dome (left) and from the truncated part of the dome (right). See caption of Fig. 3 for the label R, T, D convention

In Fig. 5 we show the results of modeling of the primary reflection from the top of the PVC material using the DKIM, and more specifically the single-scattering seismograms for synthetic data along line Y200, together with the total seismograms from laboratory data obtained with the NB and BB transducers. Direct observation of the laboratory and single-scattering seismograms shows a good fit in the modeling of reflection and diffraction events. For the case of the NB transducer, its focused beam cannot illuminate the out-of-plane structures which, therefore, are not observed on laboratory and synthetic data. On the contrary, for the case of the BB transducer, diffractions at the edges of the topographic structures and reflections from the out-of-plane structures and the fault can be clearly observed in the synthetic data, in accordance with the experimental data, even if the acquisition line does not cross the top of the structures. Though the data are zero-offset, they exhibit interesting 3D effects. Nevertheless, the steep slopes of the truncated pyramid remain invisible.

Fig. 5

Laboratory seismograms (left) and single-scattering synthetic seismograms (right) obtained with the narrow- (up) and the broad-beam (down) transducer along line Y200

In Fig. 6 we provide a more detailed qualitative comparison of the laboratory and synthetic traces, corresponding to three chosen source positions marked with arrows in Fig. 5. The source positions are the same for both NB and BB transducers. For the BB transducer, the left trace corresponds to the reflection from the slope of the out-of-plane full dome, the flat part of the model and the slope of the out-of-plane truncated dome (Fig. 6). The middle one corresponds to the reflections from the slope of the out-of-plane full dome and truncated dome, the flat part of the model and diffraction effects produced by the truncated part of the dome. The right one corresponds to diffraction effects produced by the fault, the reflection from the flat part of the model and diffraction effects produced by tips of the truncated pyramid. One can note a good fit of the traces in terms of the shape and the phase of the signal, but discrepancies in terms of amplitude, in particular, for traces obtained with the BB transducer, which contain diffractions and reflections from the out-of-plane structures.

Fig. 6

Comparison between laboratory (red) and numerical (blue) traces obtained with the narrow- (up) and broad-beam (down) transducer along line Y200, and corresponding to the source positions marked with arrows in Fig. 5. See caption of Fig. 3 for the label R, T, D convention. Label R3 stands for wave reflection from the full dome, whereas label DF (respectively DP) stands for wave diffractions produced by the fault (respectively, the pyramid)

Let us provide a quantitative analysis by performing a numerical comparison of these three traces selected. We use several error norms for quantitative analysis of the misfit between the synthetic and laboratory data, namely the single-valued normalized cross-correlation coefficient (Wilks 2011), the cross-correlation coefficient together with the root mean square (rms) misfit (Fehler and Keliher 2011), and misfit criteria EM and PM based on the time–frequency representation using the continuous wavelet transform (Kristeková et al. 2009).

Let us consider two seismograms s(t) and s_ref(t), where s(t) is the numerical signal, s_ref(t) the laboratory (reference) signal, and t time. The single-valued normalized cross-correlation coefficient cc is defined as (Wilks 2011)

$$\text{cc}=\frac{{\sum }_t\left[s(t)-\overline{s}\right]·\left[s_{\text{ref}}(t)-\overline{s_{\text{ref}}}\right]}{\sqrt{{\sum }_t{\left[s(t)-\overline{s}\right]}^2}\sqrt{{\sum }_t{\left[s_{\text{ref}}(t)-\overline{s_{\text{ref}}}\right]}^2}}$$

(1)

where $\overline{\text{s}}$ and $\overline{s_{\text{ref}}}$ are the mean values of the corresponding series. That cross-correlation coefficient estimates the degree to which two signals are correlated in terms of phase and satisfies −1 ≤ cc ≤ 1. The equality cc = 1 is for perfect correlation, cc = 0 for uncorrelated series, and cc = −1 for negative correlation.

The rms misfit defined as

$$\text{rms}=\sqrt{\frac{{\sum }_t{\left|s(t)-s_{\text{ref}}(t)\right|}^2}{{\sum }_t{\left|s_{\text{ref}}(t)\right|}^2}}$$

(2)

gives the quantitative measure of the difference between the signals without indicating the reason of the difference. To define the phase and amplitude misfits separately, the cross-correlation coefficient cc can be used to determine the phase misfit (PM) and the rms misfit to determine the amplitude misfit after shifting the seismogram by the PM (Fehler and Keliher 2011).

In addition, we consider the misfit criteria based on the time–frequency representation using the continuous wavelet transform introduced by Kristeková et al. (2009). They define a time–frequency envelope difference

$$\mathrm{\Delta }E(t,f)=|W(t,f)|-|W_{\text{ref}}(t,f)|$$

(3)

and a time–frequency phase difference

$$\mathrm{\Delta }P(t,f)=\left|W(t,f)\right|\frac{arg\left[W(t,f)\right]-arg\left[W_{\text{ref}}(t,f)\right]}{\mathit{\pi }}$$

(4)

where W(t, f) and W_ref(t, f) are the complex functions of the continuous wavelet transforms of the numerical and reference signals, respectively. We also use single-valued measures of fit based on the error norms defined above, i.e., a single-valued envelope misfit

$$\text{EM}=\sqrt{\frac{{\sum }_f{\sum }_t{\left|\mathrm{\Delta }E(t,f)\right|}^2}{{\sum }_f{\sum }_t{\left|W_{\text{ref}}(f,t)\right|}^2}}$$

(5)

and a single-valued PM

$$\text{PM}=\sqrt{\frac{{\sum }_f{\sum }_t{\left|\mathrm{\Delta }P(t,f)\right|}^2}{{\sum }_f{\sum }_t{\left|W_{\text{ref}}(f,t)\right|}^2}}$$

(6)

We provide that comparison of the misfits between the three selected laboratory and numerical traces obtained with the NB and the BB transducer in Table 1. The qualitative results shown above are confirmed quantitatively by the high values of the normalized cross-correlation coefficient and the low values of PM on the one hand, indicating that the phase fit is good, and by relatively high values of the rms misfit and EM on the other hand. However, the results are slightly worse for the BB transducer than for the NB transducer, indicating that there are phase shifts though the amplitude fit is good. The time shifts observed for the reflections from the slope of the fault are due to a wrong tilt used for modeling.

Table 1 Comparison of the misfits between laboratory and numerical traces obtained with the narrow- (up) and broad-beam (down) transducer along line Y200, and corresponding to the source positions marked with arrows in Fig. 5

5 Conclusions and future work

We have tested an alternative approach for benchmarking numerical methods for 2D and 3D wave propagation. This approach consists in comparing synthetic data obtained using numerical modeling to laboratory data obtained for a known configuration. We have obtained the laboratory data by laboratory scale measurements of zero-offset reflection of broadband pulses, generated by narrow-beam and broad-beam sources, from a scaled representation of a geological model with strong 3D topographies immersed in water. The diffraction effects produced on the wavefields by the complicated features of the model, together with the existence of shadow zones, make the laboratory experiments under controlled conditions of interest. We have computed synthetic data by means of a spectral-element method and a discretized Kirchhoff integral method. The comparisons between synthetic and laboratory data exhibit a good qualitative and quantitative fit in terms of time arrivals and acceptable fit in amplitudes for both narrow-beam and broad-beam datasets. Nevertheless, they also reveal the necessity to switch to 3D in future work and use our numerical tools to model 3D zero-offset wave reflection from strong topographic structures and to better describe the reflection phenomena with more complex reflection coefficients.

The data sets obtained during this measurement campaign seem promising for future use as a real-data benchmark for 2D and 3D model comparisons. The next step will consist in a multi-offset experiment with broad-beam transducer. In future work, we plan to perform numerical cross-validation of 2D and 3D numerical methods and techniques on the obtained data in order to analyze the respective limitations of each method and to choose the right strategy for further development of the methods. This is the final goal of our project called BENCHIE (http://www.benchie.cnrs-mrs.fr/).