Bence ; Solymosi
email: bence.solymosi@gmail.com

Favretto

Cristini

Arntsen

Nathalie Favretto-Cristini

Vadim Monteiller

Dimitri Komatitsch

Paul Cristini

Børge Arntsen

Bjørn Ursin

How to adapt numerical simulation of wave propagation and ultrasonic laboratory experiments to be comparable? A case study for a complex topographic model

Keywords: ultrasonic, finite element, modeling, diffraction wave propagation, spectral-element modeling, laboratory experiments, curved interfaces, diffraction, transducer characteristics

Numerical methods are widely used in seismic exploration to simulate wave propagation, however the algorithms are based on various assumptions. The accuracy of numerical simulations is of particular interest in case of realistic geological setups. The direct comparison of numerical results can have limitations, and an alternative approach can be the comparison of synthetic results with experimental data, obtained for a small-scale physical model in laboratory conditions. Laboratory experiments are repeatable and provide high-quality data for a known configuration. The goal of this work is to provide a possible workflow to adapt the numerical simulations and the laboratory experiments to each other, such that the two can be easily compared with high accuracy. The model is immersed in a water tank and a conventional pulse-echo technique is used to collect reflection data in both zero-offset and offset configurations. We use a spectral-element method for the numerical modeling. The model geometry is implemented using a nonstructured mesh and the computational cost can be optimized using larger elements and higher-order basis functions. The real source transducer characteristics are implemented based on a new approach: laboratory characterization of the impulse response, followed by an inversion step to obtain a numerically-equivalent source. The comparison of the zero-offset synthetic and laboratory results reveals an excellent fit in terms of arrival time, phase and amplitude. Minor amplitude mismatches may be attributed to the noise recorded in the laboratory data, and to the possible inaccuracy of the proposed source implementation. Comparison of simulated and laboratory offset traces also exhibits a good fit in general, but with significantly less accuracy for some arrivals than in the zerooffset case. This can be mainly attributed to the inaccuracies of the transducer positions

INTRODUCTION

Numerical simulation of seismic wave propagation is widely used for environmental and industrial applications for subsurface evaluation in seismic exploration (e.g., for survey design, data processing and interpretation) [START_REF] Robertsson | Introduction to the supplement on seismic modeling with applications to acquisition, processing, and interpretation[END_REF]. Simulation is also a core tool of seismic imaging and inversion [START_REF] Virieux | A review of the spectral, pseudospectral, finite-difference and finite-element modelling techniques for geophysical imaging[END_REF]. Conventional methods (e.g., ray-tracing, Kirchhoff integral and finite-difference methods), widely used in seismic exploration, are efficient to simulate realistic wavefields in environments with simple structures and slowly-varying material properties. However, difficulties arise for environments with large and rapid structural changes, due to shadow zones and (multiple) diffractions. Thus, different methods have been developed to improve seismic modeling in realistic geological environments, including steeply-dipping faults, curved interfaces, salt bodies, etc (e.g., [START_REF] Mittet | On the internal interfaces in finite-difference schemes[END_REF]. Before real applications, new methods are typically tested against other numerical methods using synthetic configurations. Several projects have focused on the comparison and validation of different numerical results (e.g., [START_REF] Igel | The COSY project: verification of global seismic modeling algorithms[END_REF][START_REF] Moczo | Comparison of numerical methods for seismic wave propagation and source dynamics-the SPICE code validation[END_REF]Fehler & Kehiler, 2011;[START_REF] Chaljub | 3-D numerical simulations of earthquake ground motion in sedimentary basins: testing accuracy through stringent models[END_REF].

Since each numerical algorithm is based on mathematical or physical assumptions, their direct comparison for realistic and complex models can have limitations, as it can be difficult to determine the one that gives the best approximation of a physically unknown solution. Therefore there is a strong interest in using physical datasets to benchmark synthetic results. However, since in real life the subsurface of the Earth is never accurately known, the synthetic results cannot be directly compared with any seismic or seismological dataset from real measurements.

An alternative approach to test and validate the performance of numerical methods in realistic cases can be the comparison of synthetic results with experimental data, obtained for a small-scale physical model in laboratory conditions. This implicitly assumes that the scaled physical mechanisms are identical to those at seismic scale (Ebrom & MacDonald, 1994), and this assumption is fulfilled in case of the linear wave equation. Considered as obsolete in the 1990's -essentially due to the drastic increase in computing capacities -, laboratory experiments have recently been re-introduced into the ideas-to-applications pipeline. Laboratory can be considered as a halfway house between numerical modeling and field observations. Indeed, laboratory experiments are repeatable, more controllable than real seismic surveys, versatile in terms of acquisition setup and provide high-quality data for a known configuration. Furthermore, similarly to real seismic acquisitions and unlike some numerical data, laboratory measurements contain random and signal-generated noise, multiples, mode conversions, and uncertainties due to position inaccuracies. As these sources of noise and uncertainties can be better assessed than in case of field datasets, laboratory experiments also provide a higher fidelity data than real seismic surveys. Therefore, laboratory experiments are a valuable tool to validate numerical simulations against real physical datasets.

Initially, small-scale physical modeling was extensively developed for a better understanding of wave propagation phenomena (e.g., Wapenaar & Berkhout, 1987;[START_REF] Pant | Physical and numerical model study of diffraction effects on seismic profiles over simple structures[END_REF] and for the validation of theoretical predictions [START_REF] Favretto-Anrès | Excitation of the Stoneley-Scholte wave at the boundary between an ideal fluid and a viscoelastic solid[END_REF]. Laboratory experiments are still used nowadays to investigate physics that is not sufficiently understood to be numerically modeled with confidence [START_REF] Cooper | The wedge model revisited: A physical modeling experiment[END_REF][START_REF] Stewart | Physical modeling of anisotropic domains: Ultrasonic imaging of laser-etched fractures in glass[END_REF][START_REF] Ekanem | P-wave attenuation anisotropy in fractured media: A seismic physical modeling study[END_REF]Xu et al., 2016;[START_REF] Chang | Azimuthal variation of converted-wave amplitude in a reservoir with vertically aligned fractures -a physical model study[END_REF]. Data from laboratory experiments are also used as input to inverse problems [START_REF] Pratt | Seismic waveform inversion in the frequency domain[END_REF][START_REF] Favretto-Anrès | Identification of shear wave parameters of viscoelastic solids by laboratory measurements of Stoneley-Scholte waves[END_REF][START_REF] Bretaudeau | 2D elastic full-waveform imaging of the near-surface: application to synthetic and physical modelling data sets[END_REF][START_REF] Chai | Reflectivity inversion for attenuated seismic data: Physical modeling and field data experiments[END_REF], to test new data processing algorithms [START_REF] Campman | Imaging and suppressing near-receiver scattered surface waves[END_REF], and in time-lapse 3D studies [START_REF] Sherlock | Time-lapse 3D seismic physical 53 modeling[END_REF]. More recently, small-scale modeling approaches have been developed as tools to validate numerical modeling and seismic imaging methods in the context of onshore and offshore seismics [START_REF] Bretaudeau | Small-scale modeling of onshore seismic experiment: A tool to validate numerical modeling and seismic imaging methods[END_REF]Tantsereva et al., 2014;[START_REF] Favretto-Cristini | Numerical modeling of zero-offset laboratory data is a strong topographic environment: results for a spectral-element method and a discretized Kirchhoff integral method[END_REF]. In particular, Tantsereva et al. (2014) evaluated the ability of a 3D discretized Kirchhoff integral method (DKIM) to accurately simulate complex diffractions using a zero-offset laboratory dataset, measured for a small-scale model with strong topography, immersed in a water tank. Comparisons of numerical and laboratory datasets showed that the DKIM could correctly reproduce the wavefield, except in the vicinity of secondary shadow boundaries created by the interaction of the edges of the topographic structures. As a follow up of the work of Tantsereva et al. (2014[START_REF] Favretto-Cristini | Analysis of wave scattering from a viscoelastic layer with complex shape: accepted[END_REF] quantitatively analyzed the effect of multiple scattering and surface curvature on the wavefield, in order to define the cases where these effects may be neglected in the numerical modeling without a significant loss of accuracy. These works clearly show the importance of laboratory datasets as part of the benchmarking options for numerical algorithms.

An efficient benchmarking procedure requires a carefully chosen and jointly-adapted approach to both the laboratory experiments and the numerical modeling. On the one hand, the numerical tools must be adapted to the experimental configuration (e.g., implementing the real source characteristics, material properties and acquisition conditions). On the other hand, the laboratory experiments have to be carried out in accordance with the capabilities of the numerical tools (e.g., choosing the acquisition geometry such that the future computational cost is the lowest possible). In addition, experimental uncertainties must be identified, in order to be subsequently minimized as The goal of this paper is to address the above challenges by providing a possible way to adapt both the numerical simulations and the laboratory experiments such that the two can be easily compared with high accuracy. To achieve this goal, we resorted to spectralelement modeling (SEM) as a full-wave method, which is not based on strong underlying assumptions (such as for example the DKIM). SEM also becomes more and more popular with time in the seismic community, since it is well-suited for high-performance computing (Komatitsch et al., 2003). As it combines the accuracy of a pseudo-spectral method with the flexibility of a finite-element method, SEM allows the handling of complex geometries by using a non-structured mesh and different element sizes in the computational domain [START_REF] Komatitsch | Spectral-element simulations of global seismic wave propagation-I. Validation[END_REF]Oliviera & Seriani, 2011). Since the necessary trade-off between accuracy and computational cost is usually sought-after, we are also interested in the ability of SEM to accurately simulate complex 3D wavefields including (multiple) diffractions at the lowest possible computational and man-hour cost.

Few papers are devoted to the comparison of SEM with other numerical methods (e.g., [START_REF] Capdeville | A new coupled spectral element and modal solution method for global seismology: A first application to the scattering induced by a plume-like anomaly[END_REF][START_REF] Capdeville | Coupling the spectral 48 element method with a modal solution for elastic wave propagation in global Earth models[END_REF][START_REF] Moczo | On accuracy of the finite-difference and finite-element schemes with respect to P-wave to S-wave speed ratio[END_REF][START_REF] Chaljub | 3-D numerical simulations of earthquake ground motion in sedimentary basins: testing accuracy through stringent models[END_REF][START_REF] De Basabe | A comparison of finite-difference and spectralelement methods for elastic wave propagation in media with a fluid-solid interface[END_REF]. [START_REF] Pageot | Improving the seismic small-scale modelling by comparison with numerical methods[END_REF] recently compared laboratory and synthetic datasets in an onshore configuration to investigate surface wave propagation and amplitude transformation between 2D and 3D. But to the best of our knowledge, no paper uses laboratory experiments as reference datasets for comparison in offshore configuration.

In this paper, we thus compare 3D spectral-element simulations with laboratory measurements in both 3D zero-offset and 3D offset offshore reflection configurations for the so-called Marseille-Benchie model. To our knowledge, it is the first time that such a comparison for the general offset case is presented. The Marseille-Benchie model is

F o r P e e r R e v i e w
based on French's model [START_REF] French | Two-dimensional and three-dimensional migration of modelexperiment reflection profiles[END_REF] and also includes additional structures with steep flanks, sharp edges, corners and curved interfaces (Tantsereva et al., 2014). It enhances multiple reflections and diffractions, as well as shadow zones and interactions between different structures. This complexity provides a significant challenge to any numerical method to reproduce the wavefield.

The paper is organized as the following. Section "Small-scale seismic modeling" briefly describes the laboratory experiments, i.e., the small-scale model and the experimental setup. The laboratory datasets obtained in zero-offset and offset configurations are also presented, together with their physical interpretations. Section "Numerical modeling" is devoted to a short description of SEM. More importantly, we also discuss in details the input data/parameters, namely the model geometry, the material properties, and the characteristics of the source and receiver transducers. In "Data comparisons", we compare the laboratory and numerical datasets, and then discuss the results in "Discussion".

Small-scale geological model

The Marseille-Benchie model contains various complex topographic features, such as a dome, a truncated dome, a truncated pyramid and two flat parts separated by a ramp (Figure 1). The model is made of polyvinyl chloride (PVC) and has a size of 600 mm x 400 mm, corresponding to 12 km x 8 km at seismic scale. Its thickness varies between 30 and 70 mm, depending on the geometry. The measured material properties of the PVC are Vp=2220±10 m/s, Vs=1050±10 m/s, ρ=1412±17 kg/m 3 . PVC is considered to be homogeneous and isotropic for the frequency range of interest of this work (250-650 kHz). For example the P-and S-wave velocities were found to be quasi-constant for this frequency range, suggesting a negligible dispersion. Attenuation was measured using the amplitude of transmitted monochromatic sine waves at different frequencies through different PVC samples of varying thickness. Following the measurements, the attenuation can be transformed into quality factors: 30 < Qp < 60 and 27 < Qs < 31 for P-and Swaves, respectively. As that is known, the measurement of attenuation is one of the most difficult laboratory tasks among all the material property characterizations. According to the measurements and accounting for the uncertainties of this parameter, we consider the PVC to have the same quality factor between 250-650 kHz which is in good agreement with the literature [START_REF] Favretto-Anrès | Excitation of the Stoneley-Scholte wave at the boundary between an ideal fluid and a viscoelastic solid[END_REF]. For further details on the model, see also Tantsereva et al. (2014).

Experimental setup

Since our work is concerned with offshore measurements, the model was immersed in A conventional pulse-echo technique was used to collect reflection data in both zerooffset and offset configurations (Figure 1). Zero-offset measurements were performed by 11 F o r P e e r R e v i e w using a custom-made Imasonic® transducer as both the source and the receiver. It has a diameter of 3 mm and was located 180±0.5 mm above the flat part of the model (corresponding to 3.6 km at seismic scale). The source signal generated by this transducer and its associated frequency spectrum for the frequency range of interest of this work are shown in Figure 2. The transducer has a dominant frequency of 500 kHz and a broad-beam radiation pattern, as the width of the main lobe is 35° at -3 dB. This radiation pattern allows for a large area to be illuminated and therefore more 3D effects to be captured, e.g., interaction of waves with multiple topographic features and multiple wave scattering [START_REF] Favretto-Cristini | Analysis of wave scattering from a viscoelastic layer with complex shape: accepted[END_REF].

Offset measurements were performed by using two transducers, namely the abovementioned transducer as the source and an omnidirectional Teledyn Reson® hydrophone as the receiver. Both transducers were located 150±0.5 mm above the flat part of the model. The hydrophone has an active diameter of 4 mm and its sensitivity is constant between 50-800 kHz. The source can be tilted with various angles to illuminate different parts of the model or to enhance some particular effects, such as shadow zones.

Laboratory datasets

We performed acquisitions along x-lines (Figure 1) with a sampling distance of 0.5 mm (corresponding to 10 m at seismic scale). The collected dataset thus consists of numerous parallel profiles, providing a set of reflection data for a dense grid (with a receiver spacing of 0.5 mm in both directions). In order to enhance the signal-to-noise ratio, 256 acquisitions were performed at each grid point and then the individual registrations were stacked. The speed of sound in water was 1485±0.1 m/s and 1484±0.1 m/s for zero-offset and offset measurements, respectively, due to the different water temperatures during the two measurements. The tilt angle of the source transducer was 39° for the offset measurements.

For the sake of brevity, we focus only on one line in this paper, located above the full dome, the truncated pyramid and the ramp (red line in Figure 3). This acquisition line provides complex diffraction effects, multiple reflections and arrivals corresponding to curved interfaces. Figure 4 shows the laboratory zero-offset dataset for the study line, overlaid by the interpretation of the recorded events. For the sake of clarity, only the main events are shown in Figure 4. Events a) and c) represent the arrivals related to the top surface of the PVC, corresponding to the pyramid and the dome, and to the flat parts, respectively. Event a) mainly consists of reflections (see positions between 90-115 mm for the pyramid and 95-410 mm for the dome). There are also some diffraction hyperbolas corresponding to the edges of the top surface of the pyramid, for positions below 90 mm and above 115 mm. Moreover due to the angle of the truncation, there are two smaller hyperbolas in the center of the pyramid, partly overlapping with the reflections from the top flat surface of the object. Similarly to event a), event c) mainly

F o r P e e r R e v i e w
consists of reflections, as well as some diffraction hyperbolas related to the junctions of the the pyramid and the dome with the flat part. Events b) and d) correspond to reflections from the bottom of the PVC, below the pyramid and the dome, and below the flat part, respectively. Since Figure 4 is a time section, a classical velocity pull-up effect can be seen in the different arrival times of events b) and d). Indeed, depending on the overburden, reflections from the same horizontal bottom surface of the PVC arrive at different times. Event e) represents reflections from the small truncated dome (Figure 1). This out-of-plane arrival, which is due to the broad-beam radiation pattern of the source transducer, was not considered during the simulations. Event f) shows reflections from the ramp, which was partly accounted for during the simulations. Note that the diffractions on the right side of the section after 200 μs and for positions between 300-420 mm, are related to the side of the model.

NUMERICAL MODELING

We used a spectral-element method (SEM) for the numerical modeling of the experiments. Following a brief recall of the method, we discuss the strategy used to mesh the model geometry, the numerical implementation of the real transducer characteristics, as well as the preliminary numerical calibration of the material properties.

Recall of the spectral-element method

Here we focus only on some of the most important features of the method, and we refer the reader to [START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Fichtner | Full Seismic Waveform Modelling and Inversion[END_REF] or [START_REF] Peter | Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes[END_REF] for more details.

The strong form of the wave equation for viscoelastic materials is described as

ρ∂ t
2 u=∇⋅σ + f (1) where ρ denotes the mass density, u=u(x,t) the displacement field, x the spatial coordinates, t the time, σ the stress tensor and f the source term. The source term for viscoelastic materials can represent for example a point source:

f =-D⋅∇ δ(x-x s ) S(t) (2)
where D denotes the moment tensor, xs the source position, δ the Dirac delta distribution and S(t) the source wavelet. Under the assumption of small perturbations, the stress tensor σ is linearly related to the displacement field through the constitutive relationship (Hooke's law):

σ=c : ∇ u (3)
where c denotes the stiffness tensor, describing the elastic properties of the material.

Since the Marseille-Benchie model is made of a viscoelastic material, equation 3 has to be modified such that the stress is determined by the entire strain history:

σ (t )= ∫ -∞ t ∂ t c(t-t '): ∇ u (t ' )dt ' (4)
The technique of [START_REF] Liu | Velocity dispersion due to anelasticity: implications for seismology and mantle composition[END_REF] is used to approximate the absorption based on a set of three Zener standard linear solids. This standard approach usually assumes that the quality factor Q does not depend on the frequency. This assumption is confirmed by our laboratory measurements, which show the same Q factor for the PVC used in the model when measured between 250-650 kHz (see section "Small-scale geological model"). In practice, the implementation of the Zener model requires fitting Q in the frequency range of interest by using a set of relaxation mechanisms (relaxation times represented by points in the frequency range, and their associated weights). A linear approach can be used by setting the optimization points at given frequencies in the frequency range of interest and then optimizing the fit for the weights [START_REF] Emmerich | Incorporation of attenuation into timedomain computations of seismic wave fields[END_REF]. However, this approach does not ensure the positivity of the weights and thus the decay of the total energy.

An alternative approach is introduced by [START_REF] Blanc | Highly accurate stability-preserving optimization of the Zener viscoelastic model, with application to wave propagation in the presence of strong attenuation[END_REF] to ensure that the weights are always positive and their non-linear optimization is also more accurate in terms of fitting the constant Q approximation in the frequency range of interest. The non-linearity means that both the points and the weights are optimized to obtain a better fit of Q in the frequency range of interest. We use this approach for the simulations presented in this work.

SEM uses the weak form of the equations, which can be obtained by multiplying equation 1 with an arbitrary test function w and then integrating by parts over the total volume of the model Ω:

19

F o r P e e r R e v i e w ∫ Ω ρ w⋅∂ t 2 u d 3 x= ∫ ∂Ω n⋅σ⋅w d 2 x-∫ Ω ∇ w :σ d 3 x+ M : ∇ w (x s ) S(t ) (5)
Equation 5 can be reformulated as:

M ü(t)+ K u(t)=f (t) (6)
where M denotes the mass matrix and K the stiffness matrix. In SEM, the computational cost is optimized by combining high-degree Lagrange interpolants to represent the wavefield and Gauss-Lobatto-Legendre (GLL) quadrature to compute the integrals involved [START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF]. This combination leads to a perfectly diagonal mass matrix in equation 6, which then enables the use of an explicit time scheme that can be efficiently parallelized (Komatitsch et al., 2003;[START_REF] Carrington | High-frequency simulations of global seismic wave propagation using SPECFEM3D_GLOBE on 62 thousand processor cores[END_REF][START_REF] Vos | From h to p efficiently: implementing finite and spectral/hp element methods to achieve optimal performance for low-and high-order discretisations[END_REF]. On the one hand, SEM is efficient in handling complex geometries and fluid-solid coupling. On the other hand, the standard GLL quadrature requires a hexahedral mesh in 3D, which can be challenging to obtain in case of a non-structured mesh for realistic geological setups.

Considering only a part of the full physical domain results in the need for artificial boundaries in the simulations. In order to avoid spurious reflections from these boundaries, we use a perfectly matched layer (PML) technique (Komatitsch & Tromp, 2003;[START_REF] Festa | Interaction between surface waves and absorbing boundaries for wave propagation in geological basins: 2D numerical simulations[END_REF][START_REF] Kristek | A brief summary of some PML formulations and discretizations for the velocity-stress equation of seismic motion[END_REF]. As SEM uses the weak form of the wave equation, the PML equations have singularities, which need to be explicitly removed [START_REF] Blanc | Highly accurate stability-preserving optimization of the Zener viscoelastic model, with application to wave propagation in the presence of strong attenuation[END_REF].

In our work, we used Specfem3D-Cartesian, an open-source spectral-element software package [START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF], with explicit second-order Newmark time stepping [START_REF] Hughes | The finite element method, linear static and dynamic finite element analysis[END_REF].

20

F o r P e e r R e v i e w

Meshing

Creating a non-structured hexahedral mesh for a complex geometry is a challenging and lengthy task [START_REF] Shepherd | Hexahedral mesh generation constraints[END_REF][START_REF] Staten | Unconstrained plastering-hexahedral mesh generation via advancing-front geometry decomposition[END_REF]. We used Cubit/Trelis ® [START_REF] Blacker | CUBIT Mesh Generation Environment Users Manual[END_REF] to mesh the model. In case of a non-structured mesh, we need to keep in mind the future computational cost and the accuracy by considering three points. Firstly, the element size must be small enough to accurately model the highest frequencies to be considered. Secondly, the size of the different elements in one material should be as equal as possible (depending on the geometry) to avoid too small elements. This is important because smaller elements require smaller time steps (i.e., higher computational cost) and more memory. Finally, one needs to avoid creating too distorted/elongated elements, which could result in meshes of poor quality or even containing elements with negative Jacobians, making the simulation unstable. Since the model has a complex overall geometry, fully automatic hexahedral meshing algorithms could not be used. Our solution was to cut the domain into several subdomains, which were easier to handle for the meshing algorithms (Figure 3). The challenge of the task was to find the order in which the individual subdomains had to be meshed, such that the entire computational domain could be meshed at the end. Due to these difficulties -and also to reduce the computational cost -, we considered only a part of the full model for the numerical simulations (Figure 3). As a result, the truncated dome was excluded from the simulations, as this object proved to be too complex to be meshed, mainly due to its small dimensions combined with sharp edges and narrow corners. With our decomposition strategy, we first obtained a mesh containing about 15.6 million elements, including the water layer above the PVC. Figure 6 the model with the applied decomposition and meshing strategy. The maximum edge length is about four times bigger than the smallest one, being about 1.6 and 0.4 mm, respectively. It is important to note that the accuracy of the spectral-element simulations is not directly constrained by the element size, but rather by the number of GLL points per wavelength. The initial meshing strategy was designed for 4 th order polynomial basis functions, requiring 5 GLL points per wavelength, which is approximately 5 GLL points per element [START_REF] Mulder | Spurious modes in finite-element discretization of the wave equation may not be all that bad[END_REF]. Considering the minimum velocity of the model -namely, 1050 m/s for the S-waves in PVC -and the maximum target frequency -namely, 650 kHz -, the goal was to have all the edge lengths below 1.6 mm. We present an optimized meshing strategy in the Appendix, which allowed us to reduce the number of elements to about 1.4 million and the computational cost by a factor of 4 to 6.

Numerical implementation of the real transducer characteristics

One of the most important task in our work was to accurately implement the real characteristics of the transducers in the numerical simulations. Since it is omnidirectional and has a frequency-independent sensitivity, the Teledyn Reson® hydrophone (used as the receiver for offset measurements) was implemented as a point receiver. However, the

Imasonic® source transducer has a unique radiation pattern that cannot be described by the classical analytical formulation of Zemanek (1971) or by an approximate radiation pattern valid only for the dominant frequency (Tantsereva et al., 2014). The latter does not provide a good solution, as other frequencies also contribute to the radiation pattern.

Therefore, we propose a new, two-step approach to implement the real transducer characteristics in numerical simulations, accounting for all frequencies. The procedure consists of the laboratory characterization of the source followed by an inversion step in order to obtain a numerically-equivalent source.

The characterization of the radiation pattern of the source transducer was performed in a water tank. The source transducer was connected to a pulse generator and kept fixed, while the hydrophone was moved around the source to record the impulse response of the source at every 0.2°, covering an angle range of 200° at a constant distance of 259 mm (Figure 7). The recorded dataset was the input for the subsequent inversion step. Here, we need to distinguish the procedures to obtain the numerically-equivalent source in zero-offset and offset configurations, as they require a slightly different approach (Figure 8). We first introduce the procedure for the offset case, and the differences for the zerooffset case are discussed afterwards. As most of the numerical tools can use point sources, we built an equivalent disk source of several point sources. This approach is based on the theory of wave superposition [START_REF] Koopmann | A method for computing acoustic 51 fields based on the principle of wave superposition[END_REF]). The disk is described with three parameters: 7) where Φ (k) denotes the cost function after the k th iteration, and mij and cij are the recorded and the calculated impulse responses at the i th time sample and j th hydrophone position, respectively. During the inversion step, some parameters, namely the radius of the equivalent disk, the number of point sources distributed on the disk, the number of layers in the disk, and the number of iterations were tested to find the best fit with the measured data. We found that a radius of 3 mm with 253 point sources distributed on only one disk layer gives the best fit between the measured and the inverted radiation patterns (Figure 9). Using more than 10 000 iterations did not result in any significant change in the results. The comparison of simulated traces using the inverted equivalent source with the measured traces shows an angle dependency in the goodness of fit (Figure 9). Indeed, the main lobe of the radiation pattern (corresponding to ±30° around the center of the transducer) is accurately recovered. For example, the correlation coefficient between the measured and simulated traces at 0° and 20° are 0.99 and 0.98, respectively (Figure 10).

Φ (k) = ∑ i ∑ j (m ij -c ij ) 2 (
However, outside of the range of ±30°, the goodness of fit drastically decreases with increasing angle from the center of the transducer (Figure 9). For example the correlation The procedure to obtain the numerically-equivalent source in the zero-offset configuration is quite similar to the previous one, but requires a deconvolution process before the inversion step (Figure 8). This is due to the fact that for zero-offset measurements the same source transducer is used twice (first as the source, and then as the receiver). Therefore an extra measurement is needed, using only the source transducer to measure the reflected wavefield from the water-air interface at normal incidence (at a distance of half of 259 mm). The deconvolution process is done in the frequency domain after a Fourier transform of all the recorded signals. Let us denote the laboratory trace recorded by the hydrophone opposite to the source transducer as trace #1 (Figure 7), and the trace recorded by the source transducer corresponding to the reflected wavefield from the water-air interface at normal incidence as trace #2 (Figure 7). The process can be The dataset derived using this deconvolution process is the input to the inversion process, which is the same as described above for the general offset case. According to our tests, the best fit between the measured and the inverted radiation patterns can be obtained with a disk radius of 6 mm, using 253 point sources and only one disk layer.

Similarly to the offset case, using more than 10 000 iterations did not provide better results. The angle dependency in the goodness of fit is also valid for the zero-offset case.

However, the range of the accurate fit is broader (about ±35° around the center of the transducer). This is due to the fact that (relatively) more energy is focused in the central beam in this configuration, thus the outer region (with lower signal-to-noise ratio) has less influence on the inversion of the main lobe.

Numerical calibration of material properties

Since the characterization of material samples yields a range of values for each measured property, an initial calibration was necessary to find the values to be used for the final simulations. In order to have the least geometrical effect on the calibration, a test point above the flat part of the model was chosen (black cross in Figure 11), far from the sides of the model to avoid edge effects as much as possible. For this point the laboratory trace was considered to be the reference and the following properties were tested during the calibration simulations: velocity and Q factors for both P-and S-waves for the PVC.

The densities of PVC and water, and the velocity of P-waves in water were considered to be known.

We note here that the investigated parameters have a combined effect on the amplitude, phase and arrival times, so the result of this calibration provides one possible solution in the parameter-space. First the velocities of P-and S-waves for the PVC were calibrated, then all the remaining differences between the reference laboratory trace and

Comparison of zero-offset data

Here, we consider more specifically three traces of the laboratory zero-offset section,

shown in Figure 4 and corresponding to the red line in Figure 3. These traces (labeled as A, B, and C in Figure 4) are of particular interest, as they contain diffracted waves generated by the feature edges, as well as reflections from the flat and curved surfaces of the model.

Trace A is located above the lower edge of the pyramid, where its flank meets the flat part (Figure 12). It is composed of several groups of reflections from the side of the pyramid (denoted as event A1), from the flat part (A3 and A4), from the bottom of the model below the flat part (A8 and A9) and below the dome (A10). Several groups of diffracted waves are also present, corresponding to diffractions from the upper edge (A2)

and the lower edge (A7) of the pyramid, and from the edge of the dome (A6). Event A5 is related to the truncated dome, which is not included in the simulations.

Trace B is located halfway between the pyramid and the dome (Figure 13). It is composed of several groups of events associated with reflections from the side of the dome (events B1 and B2) and from the side of the pyramid (B3), from the flat part (B4), and from the bottom of the model below the flat part (B6). Event B5 is the superposition of diffractions from the upper edge of the pyramid and from the edge of the dome. Event B7 corresponds to the superposition of a reflection from the bottom of the model below the dome and a diffraction from the lower edge of the pyramid.

Trace C is located on the other side of the dome, compared to traces A and B (Figure andC7). Event C5 is a multiple of a diffraction from the lower edge of the dome, while event C8 is a multiple reflection.

It is important to note that the laboratory data are real data, as they contain unknown noise recorded by the acquisition system. Furthermore, as can be seen in Figure 2, the source wavelet is a long signal, composed of a main event followed by a few tens of μs long low-energy tail with small amplitudes (ringing effect). It may explain the fact that some events are composed of two parts (for instance, events A8 and A9 in Figure 12; B1 and B2 in Figure 13; C1 and C2, C3 and C4, C6 and C7 in Figure 14).

Qualitative comparisons between laboratory and synthetic traces show a quite good fit in amplitude, phase and travel time. Regardless of the trace location, synthetic traces show an almost perfect fit with the experimental data in time, phase and amplitude for the reflected events from the top and bottom of the flat part of the model. The early and late parts of these arrivals sometimes reveal minor amplitude misfits, most probably due to the low-energy second part of the source signal, which may not be perfectly reconstructed by the source inversion, as well as the uncertainties in the attenuation parameters chosen for the simulations. Events associated with the side reflections and/or diffractions from the topographic features are also well restored by the simulations (see events A1 in Figure 12, B5 and B6 in Figure 13, C1 and C8 in Figure 14). The correlation coefficients between measured and simulated traces are 0.91, 0.95 and 0.91 for traces A, B and C, respectively. These good results are due to the fact that SEM has the ability to honor the model geometry, even for tilted and curved interfaces by using a non-structured mesh.

Comparison of offset data

Here, we consider more specifically one trace of the laboratory offset section, labeled as D in Figure 5. The source location (represented by the asterisk in Figure 15) is above the flank of the dome, and the receiver (represented by the triangle in Figure 15) is located close to the ramp. Trace D contains several groups of events, mainly associated with the direct source-receiver travel path (labeled as D1 in Figure 15), reflection from the top of the truncated pyramid (D2) and reflection from the curved surface of the dome (D3).

A qualitative comparison between the simulated results and laboratory data shows a good fit in arrival time, phase and amplitude for all the events pointed out, except for the reflection from the dome (D3). Some important remarks must be highlighted here. First, the offset configurations are more sensitive to noise than zero-offset ones. Due to a more complex acquisition setup, using more electric devices and transducers, the signal-tonoise ratio of the laboratory data in offset configuration is generally lower than in the zero-offset case. Furthermore, the offset configurations are more sensitive also to uncertainties in the acquisition geometry, especially in case of a model with strongly tilted and curved interfaces. Indeed, it can be analytically shown that a small uncertainty in the source position and/or in the incidence angle of the incident wave may have a significant impact on the illumination of the curved surface of the dome. As a result, it can have a huge effect on the reflected and diffracted wavefields (see section "Discussion"). It explains the fact that event D3 in Figure 15 is quite well modeled by the simulation, but with both significant time and amplitude misfits. The correlation coefficient between the measured and simulated traces is 0.76 for trace D. Finally, as

F o r P e e r R e v i e w
mentioned above, the implementation of the numerically-equivalent source has higher uncertainties for both the low-energy later part of the source wavelet and the secondary lobes (due to the low signal-to-noise ratio of these events) than for the main lobe of the directivity pattern. Even if they carry less energy than the main lobe, the role of these higher-order lobes in the illumination of the model should not be neglected in case of offset configurations.

DISCUSSION

As shown above, SEM can reproduce laboratory zero-offset data in terms of arrival time, phase and amplitude with excellent accuracy. However, the fit between synthetic results and experimental offset data is significantly less accurate. We discuss here the possible explanations for this observation, including experimental uncertainties and numerical implementation of the directivity pattern of the source transducer, as well as the computational cost.

Laboratory data always contains noise, and offset datasets have a lower signal-tonoise ratio than zero-offset ones. Moreover, uncertainties in the transducer positions also have an effect on the data. It is more significant in the offset case, since there is a nonlinear combination of uncertainties in the source location, the source tilt angle and the receiver position. This non-linear combination makes it difficult to evaluate the role of one effect over the others. This effect is even more pronounced in the presence of strongly curved interfaces, causing wave defocusing. Let us analyze here this source of misfit and quantify the order of magnitude of the possible resulting error. As the problem is the most pronounced in case of strongly curved interfaces, we show an example for the Marseille-Benchie model using the dome.

The geometry of the problem is shown in Figure 16, and for now we suppose that the source is located exactly above the top of the dome. We are interested in the change in the offset of the arrival of a given ray due to an error in the tilt angle of the source ( γ

). Based on geometrical considerations, we can obtain the following relations: In ongoing work, we now focus on the design of an acquisition system providing much higher accuracy of transducer positions to significantly decrease these uncertainties in the recorded data. We have proposed a strategy to numerically implement the directivity pattern of the real source transducer. This strategy permits to recover the main lobe of the emitted beam, where most of the energy is concentrated. Excellent fit between synthetic results

δ x= R+ H cos 2 (θ) δθ (8)
and laboratory zero-offset data shows that the strategy is efficient in zero-offset configurations, as in that case mostly the main lobe illuminates the model and contributes to the recorded data. Nevertheless, for offset data recorded in strong topographic environments, the proposed strategy is less accurate, since both the low-energy later part of the source wavelet (see Figure 2) and the secondary lobes play an important role in the illumination of the model, even if they carry far less energy than the main event of the source wavelet in the main lobe. Indeed, in our case they may interact with the dome and the pyramid (depending on the source location), influencing the amplitude and phase of the wavefield. Future work will focus on a new approach to numerically implement real transducer properties, which can accurately account for both the main lobe and the lower- Computational cost of numerical methods is often an important point, particularly in the operational context of seismic exploration. Here, we differentiate between the manhour cost spent on the preparation of the simulations and the computational cost (i.e., the number of processors used for simulations times the actual running time of one simulation). In general, SEM has large man-hour cost, due to the non-structured hexahedral meshing. To illustrate this problem, it is interesting to note that we spent three weeks to find the optimal meshing strategy for the Marseille-Benchie model, without even finding a satisfying strategy for the truncated dome. Contrary to the common finitedifference gridding, there is no quick or cheap solution to non-structured hexahedral meshing in terms of man-hour cost.

Regarding the computational cost of the method, we used 2100 Intel® Xeon® Sandy Bridge EP (E5-2680) processor cores for the simulations. The initial meshing strategy resulted in a mesh that needed 10150 core-hours to simulate 350 μs of wave propagation (corresponding to 7 s at seismic scale). Using the optimized meshing strategy of the Appendix, this cost was reduced to 1611 core-hours with a relative element size of 2.5, or to 2538 core-hours with a relative element size of 2.25.

CONCLUSIONS AND FUTURE WORK

The goal of our work is to provide a possible workflow to adapt both the numerical simulations and the small-scale laboratory experiments to each other, such that the two can be easily compared with high accuracy. We are also interested in the ability of SEM to accurately simulate complex 3D wavefields including (multiple) diffractions at the lowest possible computational and man-hour cost. We have thus compared 3D SEM simulations with laboratory measurements in both 3D zero-offset and 3D offset offshore reflection configurations for a small-scale physical model. The model includes structures with steep flanks, sharp edges, corners and curved interfaces. This complexity provides a challenge to any numerical method to reproduce the wavefield.

Prior to the simulations, some works have focused on the input data/parameters, namely the material properties, the model geometry and the characteristics of the source and receiver transducers. The viscoelastic behavior of the material used in the model has been approximated with a set of Zener standard linear solids, while the other parameters are known from previous laboratory characterization of material samples. The real source transducer characteristics have been implemented based on a new approach: laboratory characterization of the impulse response, followed by an inversion step to obtain a numerically-equivalent source for the numerical simulations. The zero-offset measurement requires an additional deconvolution step before the inversion, as in that case only one transducer is used as both the source and the receiver.

We have suggested an optimization of the computational cost, by using larger elements in the non-structured mesh and higher-order polynomial basis functions. This technique helps to significantly reduce the computational cost, while obtaining a similar 40 Future work will focus on the development of a more accurate acquisition system to reduce inaccuracies in transducer positions during the laboratory measurements.

Furthermore, the noise level of the new acquisition system will have to be assessed (e.g. [START_REF] Valensi | Multicomponent reduced scale seismic modelling: upgrade of the MUSC laboratory with application to polarization observations[END_REF]. The proposed numerical implementation of the real transducer should also be investigated further, to account more for the low-energy secondary lobes, and the low-energy later arrivals of the source signal as well, since they may also contribute to the wavefield. The reduction of the man-hour cost of SEM, due to the lengthy meshing step, will be essential in the future. A more robust and automatic meshing tool, suitable for non-structured hexahedral meshing of arbitrary 3D geometries will be required. The laboratory datasets and the obtained numerical results should be also compared with other numerical methods, such as the commonly used finitedifference method.

Appendix

As mentioned in "Meshing", finding the appropriate meshing strategy for a complex geometry using a non-structured hexahedral mesh is a time-consuming procedure. Thus, our goal here is to reduce the computational cost without modifying the meshing strategy.

In other words, we use the same subdivision of the whole model into subdomains and the same meshing order of the different subdomains.

Our optimization takes advantage of the fact that the accuracy of the spectral-element simulations is not directly constrained by the element size, but rather by the number of GLL points per wavelength. It means that we can increase the element size and keep a similar level of accuracy of the simulations by increasing the order of the polynomial basis functions. This technique is widely discussed in the literature as h-, p-, or h-p convergence, where h stands for the element size and p denotes the polynomial order of the basis functions (e.g., [START_REF] Hughes | The finite element method, linear static and dynamic finite element analysis[END_REF]Maday & Rønquist, 1990;[START_REF] Seriani | A spectral element method for acoustic wave simulation in heterogeneous media[END_REF][START_REF] Vos | From h to p efficiently: implementing finite and spectral/hp element methods to achieve optimal performance for low-and high-order discretisations[END_REF]Oliveria & Seriani, 2011). For the initial meshing of the model we used basis functions of order N = 4, meaning that N + 1 = 5 GLL points are used in each element.

We consider the mesh presented in Figure 3 as the reference for this Appendix. As shown in section "Comparison of zero-offset data", the reference mesh yields accurate synthetic results compared with the zero-offset laboratory data. When creating a nonstructured mesh, the element size has a distribution as shown in Figure 6. This is due to the geometrical constraints on the meshing algorithm, resulting in various element sizes.

Cubit/Trelis ® uses a target element size, as the meshing algorithm aims to mesh the model such that the average of the element sizes is close to this target value. For the reference mesh, a target edge length of 1.1 mm is necessary to have the largest elements below the required threshold of 1.6 mm. For the sake of brevity, we use relative target element sizes hereafter, by considering the target element size of the reference mesh as being equal to 1. We note that small changes in the target element size result in the same distribution of the element sizes, only the size of each element is multiplied by the factor of the change. According to our experience, it remains true even for such a complex geometry as the Marseille-Benchie model if the change is at most a factor of 0.2-5. Since our meshing strategy involves the subdivision of the model into several subdomains, there is an upper limit for increasing the element size. With our decomposition strategy, using about 5 times larger elements than the reference value is the upper limit. It comes from the fact that above this value, the element size becomes too high compared with the dimensions of some of the subdomains. We emphasize that using larger elements does not result in a less accurate representation of even the curved interfaces, as curved hexahedral elements are used. It means that the numerical tool can account for curved edges and surfaces of each element, instead of only straight lines and planar surfaces [START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Fichtner | Full Seismic Waveform Modelling and Inversion[END_REF].

We tested two different approaches. The first one consists of fixing the polynomial order of the basis functions and changing the element size. The second one keeps the element size fixed and changes the polynomial order. To evaluate the results, we ran the same zero-offset simulations using the different meshes and compared the resulting synthetic traces to traces obtained with the reference mesh. To evaluate the tests we use zero-offset traces B and C (Figures 13 and14). To obtain a quantitative comparison, we computed the root mean square (RMS) difference between each synthetic trace and the reference trace. Figure 17 shows the relative computational cost and the RMS misfit for different element sizes, using 6 th -order polynomial basis functions. The results suggest that the optimum is around relative element sizes of 2.25-2.5. Using smaller elements increases the computational cost, while using larger elements increases the RMS misfit without any significant gain in the relative computational cost. Considering that the optimal element size is about 2.5 (based on Figure 17), the effect of the polynomial order has to be examined as well. Figure 18 shows the relative computational cost and the RMS misfit for different polynomial orders, using a relative element size of 2.5. The results show that the polynomial order of 8 yields the smallest RMS misfit (around 0.008). Even though the RMS misfit is somewhat higher for order 6 (around 0.05), its computational cost is more than three times lower than that of order 8. We note that the RMS misfit does not show a monotonous trend for any of the traces, neither in Figure 17, nor in Figure 18. This is probably due to the non-structured mesh. Maday & Rønquist (1990) mathematically proved that a monotonous trend in the misfit curves can be expected only for a structured mesh, but not necessarily for a nonstructured mesh. Figure 19 shows the comparison of the traces using the optimized mesh (relative element size of 2.5) and 6 th order basis functions with the reference traces. The optimized mesh provides an excellent fit with the reference solutions in general. Due to the larger elements, some minor oscillations can be seen, probably due to mesh dispersion. However their amplitude and difference from the reference solutions are negligible.

Considering that the relative computational cost is only 15.7 % of the reference simulation, we suggest to use the optimized meshing strategy. Moreover, the effect of mesh dispersion can be reduced by using a relative element size of 2.25 instead of 2.5. In that case the relative computational cost is still 25 % of the reference, in return for a a bit higher accuracy than with a relative element size of 2.5.

  measurements combined with the strong topography of the model.

  , and numerical simulations may also contribute to it.

  We carried out several laboratory experiments in a water tank, where wave propagation occurs in small-scale conditions. A scaling factor of 1:20 000 is used to scale down the real-life dimensions of typical seismic setups to the laboratory scale. Hence, an experimental frequency of 500 kHz corresponds to a seismic frequency of 25 Hz and an experimental distance of 1 mm corresponds to 20 m at seismic scale. Material properties, such as velocity, density and attenuation are not affected directly by the scaling.

  equipped with a computer-controlled acquisition system that allows for accurate positioning of the source and receiver transducers. The water temperature was continuously monitored during the measurements, providing accurate knowledge of the speed of sound waves in water (usually around 1480 m/s). Water was considered to have a density of 1000 kg/m 3 and negligible attenuation.

Figure 1 .

 1 Figure 1. Illustration of the small-scale seismic experiments performed in the water tank. The Marseille-Benchie model contains a dome (a), a truncated pyramid (b), a truncated dome (c), a ramp (d) and an elevated plateau (e). The model is illuminated by a piezoelectric transducer (on the left). In offset configuration, the wavefield is recorded by a hydrophone (on the right).

  Figure 2. Top: source signal (filtered between 250-650 kHz) generated by the source transducer and recorded by the hydrophone in opposite position, bottom: associated frequency amplitude spectrum.

Figure 5 Figure 3 .

 53 Figure5shows the laboratory offset dataset (i.e. a common shot gather) for the study line together with the interpretation of the main events. Event a) shows the direct arrival from the source and event b) the reflections from the dome. All the interpretations marked as event c) correspond to reflections/diffractions from the pyramid. Event d) illustrates the reflections from the flat part of the PVC. Event e) is related to reflections from the small truncated dome and event f) to reflections from the ramp. Event g) is not related to the model, but to spurious arrivals from the acquisition system (not considered during the simulations).

Figure 4 .Figure 5 .

 45 Figure 4. Cross-section of the laboratory zero-offset dataset, corresponding to the red line in Figure 3. The data was filtered between 250-650 kHz. Annotated events: (a) & (b) top & bottom of the pyramid and the dome, (c) & (d) top & bottom of the flat part, (e) truncated dome, (f) ramp. The vertical lines denote the zero-offset traces chosen for comparison.

Figure 6 .

 6 Figure 6. Distribution of the element size of the mesh for the part shown in Figure 3.

Figure 7 .

 7 Figure 7. Schematic diagrams of laboratory characterization of the source transducer using the source and receiver transducers (left) and the source transducer illuminating the air-water interface at normal incidence (right).

Figure 8 .

 8 Figure 8. Workflow of the procedures used to implement the numerical equivalent of the real source transducer in offset and zero-offset configurations.

  , and the number of point sources distributed on the surface. The different point sources are independent, i.e., they can have any arbitrary (smooth) source signal and are constrained only by the measured dataset and the inversion process. The goal of the inversion was to determine the source signal of each point source, such that the resulting overall source signal of the whole equivalent source is the same as the one measured at each angle. The initial guess for each point source was a constant zero pressure, and the cost function was computed using the L 2 -norm :

  and simulated traces at 60° is 0.64 (Figure10). Therefore, the secondary lobes of the real radiation pattern are less accurately recovered. This can be explained by the fact that the recorded signals corresponding to the secondary lobes have far less energy than those of the main lobe (less than -20 dB). Hence they can be overshadowed by the noise recorded in the laboratory data.

Figure 9 .

 9 Figure 9. Comparison of the measured (blue) and inverted (red) radiation patterns of the transducer in offset configuration. The amplitude is maximal opposite to the transducer (0°).

Figure 10 .

 10 Figure 10. Comparison of the measured traces (blue) with the simulated traces (red), using the inverted equivalent source, at different angles. The correlation coefficients for the corresponding traces are shown in the top right corner of each panel.

  note that trace #2 is the square of the source transfer function. Therefore the square root of trace #2 gives the source transfer function at 0°, i.e., in the opposite position to the source, -Step II: since the hydrophone is omnidirectional, the receiver transfer function is the same at each angle α. Therefore divide trace #1 by the source transfer function (Step I) to determine the receiver transfer function ( R(ω) ), -Step III: divide all the traces recorded by the hydrophone with the receiver transfer function (Step II) to get S α ( ω) for each angle α (deconvolution in time domain), -Step IV: take the square of S α (ω) for each α to obtain the zero-offset transfer function of the source transducer at each angle, -Step V: inverse Fourier transform all the obtained traces to the time domain.

Figure 11 .

 11 Figure 11. Comparison of zero-offset laboratory trace with synthetic results for the test point to calibrate the material properties before the final simulations. The position of the test point is shown with the black cross.

  is composed of several groups of reflections from the side of the dome (C1 and C2), from the flat part (C3 and C4), and from the bottom of the model below the flat part (C6

Figure 12 .

 12 Figure 12. Comparison of zero-offset laboratory trace A with synthetic results. Physical interpretation of the annotated events is provided in the text.

Figure 13 .

 13 Figure 13. Comparison of zero-offset laboratory trace B with synthetic results. Physical interpretation of the annotated events is provided in the text.

Figure 14 .

 14 Figure 14. Comparison of zero-offset laboratory trace C with synthetic results. Physical interpretation of the annotated events is provided in the text.

Figure 15 .

 15 Figure 15. Comparison of offset laboratory trace D with synthetic results. Physical interpretation of the annotated events is provided in the text. The red line shows the acquisition line, the yellow asterisk the source position and the yellow triangle the receiver corresponding to trace D.

  δ α cos (θ+ α)+ sin(θ+ α)(δ θ+ δ α)sin (α) cos 2 (θ+ α)(9)where the offset change is the sum of δ x and δ y , while R denotes the radius of the dome, H the distance between the top of the dome and the source transducer. To quantify the order of magnitude of the error in the offset, we consider H = 140 mm, R = 51.25 mm, and γ = 5°. Supposing an error of 0.5° in the source position, the point of illumination on the surface of the dome has an error of 1.3 mm. Using equations 8 and 9 above, we get δ x = 5.2 mm and δ y = 1.5 mm. It means that an uncertainty of 0.5° in the source tilt angle leads to a shift of 6.7 mm (134 m at seismic scale) in the arrival offset of the beam (i.e. in the receiver location), depending on the point of the illumination of the surface of the dome. This uncertainty may thus have a significant impact on both the arrival time and the amplitude of the reflected and diffracted events.

Figure 16 .

 16 Figure 16. The effect of the uncertainties in the source position on the wavefield.

  the zero-offset synthetic and laboratory results has revealed an excellent fit in terms of arrival time, phase and amplitude. Minor amplitude mismatches may be attributed to the noise recorded in the laboratory data, as well as to the inaccuracy of the proposed source implementation to restore the low-energy secondary lobes of the source transducer, and the uncertainties in the attenuation parameters chosen for the simulations. The comparison of simulated and laboratory offset traces has exhibited a good fit in terms of amplitude, arrival time and phase, but with significantly less accuracy for some arrivals than in the zero-offset case. This can be mainly attributed to the inaccuracies of the transducer positions during the laboratory measurements combined with the strong topography of the model, as well as to the smaller signal-to-noise ratio of the offset configurations.

Figure 17 .

 17 Figure 17. Relative computational cost (blue) and RMS misfit (green: trace B, red: trace C) for varying element size, using 6 th -order polynomial basis functions.

Figure 18 .

 18 Figure 18. Relative computational cost (blue) and RMS misfit (green: trace B, red: trace C) for varying order of the polynomial basis functions, using a relative element size of 2.5.

Figure 19 .

 19 Figure 19. Comparison of zero-offset synthetic traces using the reference mesh (blue) and the optimized mesh with a relative element size of 2.5 and 6 th -order basis functions (red): trace B (top), trace C (bottom).

F

  Wapenaar C. P.A., and A. J. Berkhout, 1987, 3D target-oriented pre-stack migration: First Break, 5(6), 217-227. Xie, Z., R. Matzen, P. Cristini, D. Komatitsch, and R. Martin, 2016, A perfectly matched layer for fluid-solid problems: Application to ocean-acoustics simulations with solid ocean bottoms, Journal of the Acoustical Society Of America, 140(1), 165-175, doi: 10.1121/1.4954736. Xu, C., D. Bangrang, and W. Jianxin, 2016, A physical modeling study of seismic features of karst cave reservoirs in the Tarim Basin, China: Geophysics, 81(1), B31-B41, doi: 10.1190/geo2014-0548.1.Zemanek, J., 1971, Beam behavior within the near-field of a vibrating piston: Journal of the Acoustical Society Of America, 49,181-191, doi: 10.1121/1.1912316.

	dispersion of spectral-element methods, Communications in Computational Physics,
	o r F o r F o r F o r F o r F o r F o r F o r
	P e e P e e P e e P e e P e e P e e P e e P e e
	r r r r r r r r
	R R R R R R R R
	e v i e w e v i e w e v i e w e v i e w e v i e w e v i e w e v i e w e v i e w
	55

Page 26 of 55 GEOPHYSICS

Page 29 of 55 GEOPHYSICS

Acknowledgements

The authors thank Zhinan Xie for his help with the implementation of stable PMLs for the 3D case, Emanuele Casarotti for his advices on non-structured meshing, and Régine Guillermin, Patrick Sanchez, Guy Rabau and Vincent Long for their help with the experiments. This work received funding from the European Union Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant, agreement no. 641943. Børge Arntsen and Bjørn Ursin received financial support from the Norwegian Research Council through the ROSE project. We also thank CNRS for its financial support through the PICS BENCHIE project. This work was granted access to the French HPC resources of TGCC under allocations no. 2016-gen7165 and 2017-gen7165 made by GENCI.

We also acknowledge the Associate Editor and the four anonymus reviewers for their careful reading, helpful suggestions and questions which all improved the paper.