Interactive comment on “ Finite difference modelling to evaluate seismic P wave and shear wave field data ” by T .

This paper is well written and organized. Some points need to be clarified, however, before publishing in my opinion. 1) In the P wave imaging workflow (Table 1), what is the finite-difference migration? It seems it is a time migration, since there is time-todepth conversion after. If so, some difference between the field P wave data image and the first synthetic P wave image may be caused by the time migration. Have you try depth migration? Please clarify this point. 2) It is a big jump in Figure 11. How to interpret the measured field data and then modify the input model? Do you use travel time tomography to update the model, or simply draw the lines by geological knowledge? It is mysterious to me. 3) Please spell out SOFI in the abstract. 4) section 6.1 could be shorter. Some paragraphs are too general, and not closely related.


Introduction
Near-surface geophysical methods constitute a nondestructive means to investigate the shallow subsurface.Especially engineering tasks, for instance geohazard assessment or groundwater prospecting, profit from structural and parametrical methods (e.g.Miller, 2013;Kirsch, 2008).In some cases, results of geophysical prospecting are compiled into 3-D models and can act as input for, e.g.groundwater flow modelling.
Among other geophysical methods, high-resolution reflection seismics constitute a valuable tool to extract structural information.The reflection seismic method using compressional waves (P waves) is an established method and often produces accurate results for near-surface tasks (e.g.Steeples and Miller, 1990;Rumpel et al., 2009;Jørgensen et al., 2012).In many cases, however, a higher spatial resolution is desired than P waves offer.This is one important reason, why shear waves have become popular.Since shear waves have a lower velocity than P waves, they offer shorter wavelengths, i.e. higher spatial resolution, within the same frequency band as P waves.Also, shear wave velocity information yields an additional parameter than the P-wave velocity information alone and offers more detailed studies of elastic properties, like for instance Poisson's ratio.The combination of P-and S-wave velocities can help characterize the lithology or pore fluid (e.g.Sheriff and Geldart, 1995).For shallow application, recent developments have been successful (e.g.Inazaki, 2004;Pugin et al., 2009a, b;Polom et al., 2010Polom et al., , 2013;;Krawczyk et al., 2013).However, we experience that shear wave reflection seismics more strongly varies in quality compared with compressional wave seismics.
Seismic modelling is an important tool to evaluate seismic field data, yet not often applied for near-surface tasks.Most studies applying wave-field modelling concentrate on deep reservoirs and crustal structures (Robertsson et al., 2007).A comprehensive approach is reported by Bellefleur et al. (2012), which we use accordingly in this study.One current method to create synthetic seismograms of complex structures is finite-difference (FD) modelling (Alterman and Karal, 1968).The advantage of the FD method is the ability to choose arbitrary heterogeneous input models without fundamental restrictions.Drawbacks like high-computational requirements have become less restricting during the past years.The FD code by Bohlen (2002) has the ability to include intrinsic seismic wave absorption, an advantage with respect to the study of unconsolidated material.
Full waveform inversion with either synthetic or real field data is an evolving field in exploration seismology.It aims at the automatic subsurface model generation from seismic field data.Although there has been a lot of progress in its evolution, and it can be applied routinely for simple subsurface structures, especially in marine environments, its application for near-surface studies is still experimental (Romdhane et al., 2011;Groos et al., 2012;Plessix, 2012).This has to do with the importance of a good starting model for the inversion process but also with the large parameter space in the near surface.Also, many seismic surveys try to avoid low frequencies in order to not generate large surface waves.Yet, the low frequencies are an important prerequisite to find the global optimum in the optimization process of the inversion (e.g.Virieux and Operto, 2009;Fichtner et al., 2011).Therefore, we focus here on direct modelling.
In this paper we describe the acquisition and processing of P-wave and SH-wave reflection seismic field data along an example profile.We generated synthetic P-wave data of the same profile with a FD algorithm on the basis of the Pwave depth section interpreted with a priori geological and geophysical information.The synthetic data were compared with the field data.Main reflectors can be reproduced.To assess the lower signal-to-noise ratio of the recorded shearwave field data, we simulated SH-wave seismograms on the basis of this model to comprehend quality variations in the SH-wave reflection field data.

Finite-difference modelling
Seismic wave-field modelling reveals the propagation of seismic waves within a subsurface model.So far, no exact analytical solution exists for the calculation of such wave fields in arbitrary media.Therefore, a number of approximation methods have been introduced over the years to solve the specific wave-equation (Carcione et al., 2002;Fichtner et al., 2011).One method is FD modelling that works for arbitrary media.The medium is divided into a grid small enough to represent the natural phenomena of elastic waves.Changes of elastic parameters for each grid point are approximated for defined time steps and simulate seismic waves travelling through the gridded model.Forces at any chosen location within the model excite the respective grid points and start the wave propagation.These forces are independent in time and place and constitute seismic sources.Several elastic properties can be extracted at any grid point and any time.These properties represent, e.g. the seismic response to a receiver at the surface or in a borehole.FD modelling simulates the entire wave-field and thus contains surface waves naturally.FD modelling has only recently become feasible for routine application, because of enhanced computing capabilities.
The FD software we used in this work is described in Bohlen (2002).This code, Seismic mOdelling with FInite differences called SOFI, is based on earlier work by Virieux (1986) and Levander (1988) for the elastic wave simulation but also on work by Robertsson et al. (1994) for the viscoelastic case.It allows the user to also consider intrinsic wave absorption (viscoelasticity, Q), and it provides an alternative rotated grid representation of the subsurface, based on work by Gold et al. (1997) and Saenger et al. (2000) for a more exact simulation of surface waves.Absorption is a common phenomenon in unconsolidated near-surface rock units, and surface waves pose a typical problem in nearsurface seismic data processing.The consideration of these phenomena in the simulation offers a better separation of the different subsurface parameters responsible for seismic field data signatures.SOFI provides full wave-field simulations for 3-D media and 2-D modelling.The 2-D codes simulate either P and SV waves within the propagation plane or SH waves oscillating orthogonally to the propagation plane.For geometrical reasons SH waves will never convert into P or SV waves in this case.This kind of simulation is often unrealistic because the subsurface has a 3-D structure in which arbitrary wave conversion takes place.Nonetheless, it is a valuable aid for the interpretation of the SH-wave field surveys.SOFI requires v P , v S , and density (ρ) as input information and optionally accepts absorption models for P and S waves.The SH version does not require compressional wave parameters.
Another important feature of the code is its ability to run simulations in parallel threads on multiprocessor computers to save time.For such parallel simulation, the model can be split up into subgrids that the different CPUs (central processing units), called processing elements (PEs), calculate independently.For every time step each PE updates, the wave field in its subgrid and parameter exchange needs to be carried out with the grid points of neighbouring subgrids.This exchange can slow down the calculation so much that the use of all CPUs and a respective number of subgrids not necessarily delivers the fastest result.Here, the minimum computational time for one test model was achieved with 48 processors and a corresponding number of subgrids.

Test site Föhr
The test site is the North Sea island of Föhr (Fig. 1) that was investigated as a pilot area in the Interreg project CLIWAT (Harbo et al., 2011), which was co-financed by the European Union.The aim of CLIWAT was to analyse the influence of climate change to groundwater systems, and one major outcome was a hydrogeological model to forecast the groundwater evolution (Burschil et al., 2012a).

Geological evolution
The investigation area is located on the German North Sea island of Föhr, which is part of the UNESCO world heritage Wadden Sea (Fig. 1).In general, the deeper subsurface was sedimented as part of the North German Basin while the landscape of Föhr was formed during the glacial and post glacial epoch (Scheer, 2012).The older sediments were deposited under marine conditions since the Cretaceous age until the youngest Tertiary.For this reason, marine clay (mica clay) is found up to Miocene age.Sedimentation changed during the Pliocene, and sandy material (kaolin sand) was deposited until a climate change marked the beginning of the Pleistocene age.Glaciations from the Baltic-Scandinavian area reworked the shallow underground by alternating processes of glacial advance followed by erosion and sedimentation.In the region of Föhr, push moraines as well as a system of tunnel valleys were formed and refilled with glacial deposits (Scheer, 2012).The great outwash plains were increasingly flooded during the Holocene so that tidal mud deposits were accumulated and formed large marshland areas.Finally, heavy floods in historical times eroded large parts and formed the present North Frisian Islands and Wadden Sea.

Geophysical and geological framework
In the project CLIWAT we accomplished a multidisciplinary geophysical acquisition (Burschil et al., 2012a).Between 2009 and 2011, we acquired seven seismic reflection profiles with P waves (in total 8 km) and three profiles with SH waves (in total 2.4 km).Five vertical seismic profiles (VSPs) were recorded with a 3C borehole geophone and the small electrodynamical vibrator system ELVIS (Polom et al., 2011) with excitation in the vertical and horizontal-transversal directions relative to the borehole location (Fig. 1).Maximum depths of five VSPs were in the range of 39-102 m, depending on the borehole.Additionally, in 2008 the island was mapped within the airborne geophysics mapping project of LIAG (Leibniz Institute for Applied Geophysics; Wiederhold et al., 2010) with the airborne electromagnetic system SkyTEM (Sørensen and Auken, 2004).The result of the Pwave seismic survey was used as a priori information for the electromagnetic data inversion (Burschil et al., 2012b).The information transfer between the different geophysical methods improved the electric resistivity model of the island.Borehole logging data were evaluated statistically regarding electric resistivity as well as seismic velocity.This allowed for a petrophysical classification of sand, till, and clay (Burschil et al., 2012a).These lithological units form structures such as push-moraines and buried valleys, which are consistent with the geological evolution of the region.A simplified hydrogeological model was compiled from all geophysical and geological results for groundwater flow modelling (Burschil et al., 2012a).m 4 Seismic reflection field data

Seismic acquisition
Seismic equipment varied with surveys due to different wave types used for exploration.For two P-wave surveys we used the hydraulic vibrator systems of the LIAG, the MHV2.7 in 2009 and the new HVP-30 in 2010.As source signal we chose a linear sweep in the range of 30-240 Hz with 10 s duration.The source excited the seismic signal on a paved street in the western part of the profile and on grassland in the eastern part.The receivers were vertical geophones (20 Hz), planted in the green strip next to the line of the source locations (cross-line offset ca. 1 m).We chose a receiver spacing of 2 m to enhance the near-surface resolution and fold.A combination of split-spread/roll along geometry was used (Fig. 2a) with a source spacing of twice the receiver spacing (4 m).High-quality data had been acquired with this geometry before.We operated up to 268 active channels with Geometrics Geode seismographs with a maximum in-line source-receiver offset of 360 m.This resulted in a fold between 47 and 72 in the main parts of the profiles.
To avoid large data volumes, we recorded the seismic traces after vibro-seismic correlation with the Pelton Vib Pro input signal.
For the shear-wave survey in 2011, an established system comprising LIAG's hydraulic shear-wave vibrator MVP-4S in combination with a 120 m land-streamer system was used (Polom et al., 2010).Horizontal 10 Hz geophones (SH-wave polarization) are mounted every metre on the land-streamer.The whole system was towed by the recording vehicle.The acquisition geometry (Fig. 2b) was similar to the P-wave surveys, except for a shorter spread and that we vibrated within the middle third of the streamer and sequentially pulled thirds of its length along the profile.Shot-point spacing was 4 m.We used a linear sweep between 30 and 160 Hz for 10 s.In contrast to the P-wave surveys, we recorded the uncorrelated traces to allow for precorrelation processing.The shear-wave seismic profile only covers the western half of the P-wave   , 6, 33-47, 2015 www.solid-earth.net/6/33/2015/profile.Here, the surface is paved, which helps avoiding the generation of surface waves.All P-wave profiles show a good signal-to-noise ratio (Fig. 3).Seismic reflections can be detected down to 1 s two-way-travel time (TWT).In contrast, the shear waves offer a smaller signal-to-noise ratio (Fig. 4).Reflection hyperbola signals are faint and within a reverberating background.Chevron patterns, also called herring-bone pattern, appear irregularly among shot gathers as part of the reverberations.The reflection signal bandwidth decreases with time.Below 0.7 s the signal vanishes, and further analysis in combination with seismic FD modelling is suggested.

Seismic processing
P-wave and shear-wave data processing differ due to different signal locations within the wave field.For the P-wave data, we set up a processing scheme (Table 1) focusing on the enhancement of the reflections (e.g.Yilmaz, 2001).Processing was carried out with Landmark's ProMAX 2D.The most important processing steps turned out to be muting of the surface waves, spectral whitening, and dynamic corrections, including dip move-out corrections.The detailed velocity analysis left a certain tolerance within the semblance calculations.We therefore included geological a priori knowledge to reduce the uncertainty so that the resulting interval velocity distribution better corresponds to the reflectors (Burschil et al., 2012a).At the end of the processing sequence, we applied a time-migration to the CDP (common-depth-point)stacked section.The migrated section was time-to-depth converted, using a single velocity function from one borehole, yielding a seismic depth section (Fig. 5).
To test the robustness of the result, a post-stack depth migration of the non-migrated CDP-stacked section was also carried out (not shown here).It provided a similar seismic image as time-migration and time-to-depth conversion.Only small differences in the location of some reflectors of up to 3 m can be identified, but the reflector continuity remains.Therefore, we continued with the time-migrated, depth-converted section.
The shear-wave data contain surface wave interferences related to the specific shear-wave reflection move-out.This is the reason why we cannot simply mute the surface wave noise and purely focus on reflection signals, as we could for P waves.To enhance the reflection signals, we applied several techniques, e.g.fk filtering, spectral whitening, and deconvolution.Finally, automatic gain control (AGC, 300 ms) and fk-filtering with low-cutoff velocities below 350 m s −1 provided the best results (Table 1).This filter also removes a wide range of the chevron patterns that are present in the lower part of the seismograms around the shot location, depicted in Fig. 4. In contrast to Polom (1997), who investigated a chevron pattern that originated in ghost sweeps, we cannot identify a comparable increase of the frequency with time in the data.Here, the chevron pattern is rather monofrequent and no ghost sweeps can be detected within the pattern.
Velocity analysis was very difficult because reflections can only be identified sporadically and can hardly be traced through to neighbouring shots.This restrains velocity analysis in CMP (common-mid-point)-sorted data.After normal move-out correction and common midpoint stacking, we converted the section to depth with a velocity function derived from VSP measurements.The result is rather monofrequent, but the till layer as well as deeper reflectors can be identified (Fig. 6).However, the quality of this shear-wave survey is less compared to the P-wave seismic results.

Synthetic data from finite-difference modelling
To understand data quality differences between shear-wave measurements and compressional wave measurements on the island of Föhr, seismic wave-field modelling is introduced here for further data analysis.We chose the 2-D P/SV-version of SOFI for P-wave simulations and the 2-D SH version of SOFI for SH-wave modelling.
To create the input velocity and density models, we interpreted the major structures of the P-wave seismic field data depth section (Fig. 5) and assigned P-wave velocities and densities according to geological and geophysical a priori information (Fig. 7; Table 2).To supplement detailed shearwave velocity information, we generated a cross-plot of velocity data (v P and v S ) from the VSP surveys on the island of Föhr.This allowed us to calculate mean and median values and the linear regression for different lithological classes.For statistical confidence, the number of classes was limited to those for which more than 20 velocity samples were available.We then picked shear-wave velocities from the linear regression line of this cross-plot and thereby defined 14 lithological units with known v P and v S values that constitute the input model (Fig. 8; Burschil et al., 2012a).
For comparison of the synthetic data with the field data, we calculated a number of single shots with geometries only slightly different from the field geometries (Fig. 2).Differences are related to modelling requirements.Although receivers were simulated 1 m apart and shots 4 m apart (Fig. 2c), only every second receiver was used for later Pwave processing.The maximum offset of shear wave channels was restricted in later S-wave processing to respect the field geometry.Additionally, the efficient use of computational capacities required the splitting of the whole model into a number of model segments of which the westernmost is depicted in Fig. 9. Our FD-modelling approach required the use of absorbing frames at the left, right, and lower boundaries of each model segment.Therefore, the model segments were 600 m long and contained an absorbing frame width of 45 m (Table 3).Instead of a free surface, we also implemented an absorbing frame at the surface (Fig. 9).For this purpose, we expanded the model top layer by 50 m that www.solid-earth.net/6/33/2015/Solid Earth, 6, 33-47, 2015  hosts the absorbing frame of 45 m.We did this to suppress surface multiples, which is the same approach that was used by Jones (2013).The effect of this absorbing frame is comparable to a weathering layer with high parameter gradients.Additionally, we introduced a low-velocity layer with a free surface to test the simulation of a weathering layer.We used a zero-phase Ricker wavelet for the simulations instead of a vibro-seismic correlated sweep signal, i.e. a Klauder wavelet.This is a practical compromise which takes into account that the field data signals are absorbed to a certain degree and do not show a white frequency spectrum.In contrast to the field data source signals, we used the same central frequency of the Ricker wavelet for P-wave simulation as well as for SHwave simulation.We used the same input models (Fig. 7; Table 2) and settings (Table 3) for P-wave modelling and SH-wave modelling, and only varied the modelled trace length.Due to lower shear-wave velocities, we had to consider longer travel times for the shear-wave signals and simulated data up to 2.5 s TWT.Table 2. Seismic velocities and densities for refined input model (cf.Fig. 7b).Shear-wave velocity was calculated according to the cross-plot relation (Fig. 8).

Computational facility
The computer we are currently using for the simulations is a DELL Poweredge 910.It accommodates four CPUs Intel Xeon E7-4870 with 2.4 GHz, 10 cores each, DDR3 memory, and hyper threading.That results in 80 usable single CPUs with SuSE 11.4 Linux.Our computer hosts 512 GB of RAM, which is by far larger than needed for our 2-D models but will be necessary for larger 3-D simulations.The computational time for modelling strongly depends on the number of modelled time steps and the size of the model.For P-wave modelling, we simulated the propagation for 1 s and 20 000 time steps (Table 3).The computational time for one single shot with the settings described was ca. 30 min on 2 CPUs using the 2-D P-SV version of SOFI.2-D shear-wave modelling with the SH version of SOFI was much faster per computational step.Simulating 2.5 s, i.e.50 000 time steps, of wave propagation took ca.60 min per shot.

P-wave modelling
In the simulated single shot data (Fig. 10), we detect direct P and SV waves as well as reflections of P and SV waves.The surface wave pattern of the field data is not present in the synthetic data, (cf.noise cone in Fig. 3).P-wave reflections can be tracked in the synthetic data, even in those regions where they are usually covered by surface waves.Also, the simulated data have an apparent higher frequency content, not showing the typical subsurface low-pass filter effect.
We applied a simple processing scheme to all 300 modelled single shots comprising amplitude control, frequency filtering, normal move-out correction, common midpoint stacking, time-migration, and time-to-depth conversion.Stacking velocities were picked via an interactive velocity analysis from CMP gathers.
In the migrated depth section (Fig. 11a), the uppermost reflector at 30 m is faint.Below, two strong reflectors mark the upper and lower boundaries of the till unit.The upper reflector can be traced through the entire section with varying amplitude.Also, a large number of reflectors with different amplitudes is imaged.In the central part of the seismogram at 150 m depth, two close reflectors mark the lower end of the complex geological units.Below, two nearly horizontal reflectors at 250 and 380 m depths are present.At the western and eastern edges of the section the migration process generated minor artefacts.Because this P-wave depth section was basically able to explain the P-wave field data, we continued with the simulation of the shear-wave field data instead of a further study of additional features present in the P-wave field data, like surface-wave ground roll.

Shear-wave modelling
Shear-wave propagation was modelled using the same input model and settings.We restricted the modelled data to a maximum offset of 80 m that corresponds to the maximum offset in the field data.Resulting shot gathers show the direct S-wave as well as several shear-wave reflections (Fig. 12).The sharp reflection signal does not change its shape with travel time.Because we used the same frequency band for the P-wave and S-wave source simulations, the stacked section (Fig. 13a) shows an even more detailed image of the subsurface due to the shorter wavelengths of the shear waves compared with the P waves.All structures of the input model are imaged with excellent resolution.

Discussion
In reflection seismic surveys normally only one type of wave is utilized, for instance P waves, shear waves, or converted waves.The use of more than one wave type or even the integration of multicomponent technology (Hardage et al., 2011) is rare but growing.This is even more pronounced for near-surface applications.Recently, a few studies reported the combination of shallow P-wave and shear-wave seismics, comparable to the field work presented in this article (e.g.Pugin et al., 2004Pugin et al., , 2009b;;Malehmir et al., 2013;Sauvin et al., 2014).Even 9C reflection seismics was tested in near-surface exploration (Pugin et al., 2009a).All of these studies successfully recorded high-resolution data of mainly good quality.
The signal-to-noise ratio of the P-wave data is similar to the P-wave data presented in this study, except in Pugin et al. (2009a).In contrast, the shear-wave data of these studies are of good quality and constitute prominent, coherent reflections.The equipment and acquisition geometry slightly vary among the reported surveys, but the combination of vibrator The P-wave field depth section is superimposed.Numbers 1-14 mark the units listed in Table 2. and land streamer is the favourite choice for shear-wave seismics.General data processing reported in some of the studies is similar to the processing we finally applied (Table 1).Sauvin et al. (2014) reported the application of elevation statics for only one of their shear-wave profiles.None of the authors reported refraction statics for shear waves.In some cases, differences in shear-wave processing are related to deconvolution and spectral whitening, which was applied by Pugin et al. (2009b) and Sauvin et al. (2014).Here, deconvolution and spectral whitening did not provide success.So far, we have not been able to reproduce the shear-wave field data with seismic modelling.A similar observation is reported by Bellefleur et al. (2012).They calculated synthetic P-wave data of four different input models from the surface to a mineralized zone and compared these data with field data.Their data excellently correlate with data from a VSP, but they consider the correlation with a surface 3-D reflection seismics data set as rather poor.They explain the poor surface data quality with the appearance of surface waves and a higher amount of scattering at inhomogeneities.Surface  longer for surface data than for VSP data.If this explanation is true, we will have to take small-scale inhomogeneities and a weathering layer as origin of surface waves into account for future studies.

Influence on data quality in land seismic surveys
In general, a number of factors can influence the quality of land seismics data, listed and illustrated in Sheriff (1975).Typical factors, one would assume to be the most likely in our case, are source strength, inappropriate coupling of sources and receivers, superimposed surface waves, multiples, scattering, and intrinsic absorption (Q).
The vertical hydraulic vibrator sources MHV2.7 and HVP-30 have proven their ability to reach reflectors at least 2 km deep (Buness et al., 2009).The LIAG shear-wave vibrator source MHV-4S has been able to generate clear reflections at least as deep as 200 m in fluviatile and marine deposits (Polom et al., 2010).Sauvin et al. (2014) used a wheelbarrow-mounted microvibrator source for the analysis of quick clays.Their data show clear reflections from 40 m below ground level within fluviatile and marine sediments.In the underlying bedrock, reflections can be traced down to 120 m depth.Malehmir et al. (2013) report clear reflections in a similar environment from at least 40 m below ground level with the same source.This means that even small sources are able to create reflections from at least 40 m below the source.Pugin et al. (2009b) used an IVI Minivib on a minibuggy carrier, which is comparable to the MHV-4S.They show SH-wave reflections down to about 50 m in glacial deposits before reaching the bedrock.In the light of these studies, the MHV-4S source can be expected to be strong enough to image the glacially overworked deposits down to 150 m.In fact, the SH-wave image shows faint reflections even from ca. 270 m depth (Fig. 6).
The ground coupling of the sources is mainly affected by the vehicle mass and the driving peak force of the system.For P-wave sources this mass acts as a gravitational hold-down force that compensates the vertical force of the vibrator unit (Sheriff and Geldart, 1995).For S-wave sources, this mass increases horizontal friction of the base plate and thus prevents sliding of the shear-wave vibrator unit on the ground.We use the rule of thumb that the gravitational hold-down force of the vibrator mass on loose ground should at least be twice the peak shear force.Experience shows that under favourable conditions, the peak force of a vibrator with a rubber base plate on paved roads can reach 70 % of the hold-down force without sliding.The P-wave hydraulic vibrator systems MHV2.7/HVP-30 have masses of 2.6 t/4 t and maximum peak forces of 27.6 kN/30 kN.Therefore, these vibrator systems should not be affected by bouncing.The shear-wave vibrator system MHV-4S has a mass of 4 t and a maximum peak force of 30 kN.During the survey, we used a peak force of 17 kN on a paved road.Therefore, the gravitational hold-down force of nearly 40 kN should have been sufficient.
Receiver coupling with the ground is another factor.Carefully planted geophones typically offer proper receiverground coupling for P-wave surveys (Krohn, 1984).The Pwave field data we acquired support this expectation.For the shear-wave survey, we used the LIAG land streamer on which geophones are directly fixed to aluminium plates that have a gravitational three-point contact to the ground.This system has proven to receive good shear-wave signals from the subsurface before (Polom et al., 2010;Krawczyk et al., 2013;Malehmir et al., 2013).On one SH-wave profile on gravel, Sauvin et al. (2014) used the same streamer but on grassland they planted the geophones in the ground.They do not report any difference in data quality.Pugin et al. (2004Pugin et al. ( , 2009a, b) , b) use a different land streamer system that works successfully as well.In the field, we took care that the geophone coupling was good.Every time a new streamer position was reached, a noise test was carried out and noisy geophones were coupled to the ground by hand.We therefore expect coupling to be sufficient.
Surface waves are an important degradation factor for shear-wave reflection surveys.Here, the reflection hyperbola often interferes with surface wave (Love wave) signals; cf.surface waves (Fig. 3) and reflection hyperbolas (Fig. 4).This factor has complicated the application of shear waves in the past and still limits the application of this method.However, the observation that a high-velocity layer at the surface suppresses the generation of surface waves in shear-wave exploration (e.g.Inazaki, 2004) was an important step to overcome this limitation in many urban applications and even in rural environments if paved or consolidated roads are present.Even though surface waves are often excited during an SHwave surveys on unconsolidated surfaces, SH-wave reflection seismics surveys have been successful in these conditions (e.g.Malehmir et al., 2013;Sauvin et al., 2014).Polom et al. (2013) identify partly dispersive Love waves that show a similar signature as the chevron pattern depicted in Fig. 4. In their case and in our case, these waves seem to be linked to the shot-point location.Even neighbouring shot points show strong variations in this respect.For instance, in Fig. 4 we show shot gathers with shot-point locations of FFID's 1126 and 1127 that are just 4 m apart.To cancel surface waves with a linear move out or mild dispersive character fk filtering is often successfully used (e.g.Polom et al., 2013;Sauvin et al., 2014).We successfully applied an fk filter as well (Table 1), but this did not improve the coherency of reflectors in the final depth section (Fig. 6).
Further issues for data quality of land seismics data are intrinsic damping (e.g.Kang and McMechan, 1994) and scattering at complex structures (e.g.Cheraghi et al., 2013;Malehmir et al., 2013).Both aspects are crucial to simulate the influence of a near-surface weathering layer.

Comparison of P-wave and S-wave simulations
We calculated synthetic seismograms for P-wave propagation through two different input models as well as shear-wave propagation through one of these models.In the following, we will compare these modelled data with the corresponding field data and discuss the observations.

P-wave shot gathers
The simulated P-wave and the corresponding field data of single shot seismograms (Figs. 3, 10) contain reflection events that show a similar basic waveform.However, the field data (Fig. 3) are much noisier, in particular before the first break as well as inside the surface-wave cone.No reflections are detectable inside this field data cone whereas in the synthetic data (Fig. 10) reflections can be traced through all parts of the seismogram.This can be explained to a large degree by the lack of a weathering layer in the models.The absence of this layer prevents the build-up of simulated surface waves and thus, a surface wave cone is missing in simulated shot records.Therefore, the refracted wave in the field data corresponds to the direct wave in the synthetic data.
Another noticeable observation is the signature of the reflections.In general, the reflection signals of the synthetic data seem to be more focused whereas many of the reflections within the field data seem to be made up of a number of oscillations/reverberations instead of single reflections (compare reflection signal at 0.3 s in Figs. 3 and 10).We expect an additional fine-layering within the 14 units of the input model to be able to reproduce this observation (Fig. 10).Perhaps, the simulation of a multilayer surface unit, related to the weathering zone might also help explaining the reverberation observation.We also simulated a low-velocity top layer (thickness = 10 m, v P = 500 m s −1 , v S = 200 m s −1 , and ρ = 2000 kg m −3 ) with a free surface, which corresponds to the unsaturated zone.In this case, the synthetic seismograms contain a large number of multiples between the first impedance contrast and the free surface.The multiples superimpose all reflection signals and cannot be identified as such in the field data.It turned out that it needs additional systematic tests to figure out a surface layer configuration that may be able to simulate this aspect of the field data.For this reason, we used an absorbing frame at the top of the model.Even though we cannot reproduce an exact copy of the reflections of the field data, the major reflections occur in good quality and quantity to use this model for additional shearwave simulations.

SH-wave shot gathers
Shear-wave shot simulations show larger differences compared with their corresponding field data (Figs. 4 and 12).The clear and continuous reflections in the synthetic data are not present in the field data.Some of the field records show the mentioned chevron pattern parallel to the first break with varying amplitude (Fig. 4) that does not show up in the synthetic data.In the model, we did not consider very shallow structures that could create Love waves.
The spectrum of the synthetic data does not change with time in the seismogram (Fig. 12).The Ricker wavelet has the same spectral shape for shallow reflections as for deep reflections.This is no surprise since up to this point we have not included any kind of signal damping in the model like intrinsic damping or sources for attenuation through scattering or interference.In contrast, the field data signals for longer travel times lose some of the higher-frequency components (Fig. 4).This indicates some kind of intrinsic or scattering attenuation.However, this cannot explain the lateral appearance variation among shot gathers of the shear-wave data.

P-wave depth section
The post-stack time-migrated P-wave depth section of synthetic data shows less noise and less amplitude variability while the corresponding field data show natural levels of different signals and thus contain more information about smallscale and internal structures (Fig. 11).The modelled depth section reproduces the major features of the field data.The till top reflector between 50 and 80 m depth appears continuous and two deep reflectors show as well.Some of the dipping reflectors in this modelled section add details that similarly appear in the field data.
The field data section (Fig. 11b) can be divided into two parts: excitation on paved street and on grassland.Within the seismic field data section, we detected a lack of resolution in the grassland part (about 900-1400 m, marked with blue arrows in Fig. 11b) that is not present in the modelled data.However, we have not implemented the structural complexity of unconsolidated grassland as a near-surface weathering layer, i.e. large velocity contrasts, inhomogeneities, and intrinsic damping, in the model so far.
In general, the reflection signatures in the field data spatially vary more strongly than in the synthetic data (Fig. 11).Sources for these observations in the field data can be intrinsic damping, scattering attenuation including fine layering, and inhomogeneous lithological units.All of these features were not included in the input model (Fig. 7).In a natural environment, complex structures or fine layers can be sources of multiples and wave conversion that subsequently can pose similar challenges to data processing like random noise signals do.Since the model consisted of comparatively large units, this kind of noise was not simulated.

SH-wave stacked sections
The SH-wave stacked sections of synthetic data and field data (Fig. 13) differ more than the P-wave sections.While no clear interpretation can be carried out for the stacked field data, the stacked section from the modelled data set well reproduces the input model.For example reflector segments occur at 90 m depth, which correspond to the interface of the subhorizontal layer at 90 m (top till layer) and at 250 m depth corresponding to the first deep reflector in the P-wave seismic section.Nonetheless, the reflection signatures in the field data are much noisier than in the synthetic data.

Summary and outlook
Shear-wave field data recorded on the island of Föhr showed less quality compared with their compressional wave counterparts.To comprehend the reasons for this quality difference, we used seismic wave-field modelling within simple models of the subsurface, using the seismic field geometry.We chose finite-difference modelling to try to reproduce the field data because of its ability to simulate the entire wave field and to allow for arbitrary input models.We come to the following conclusions.
1. We were able to simulate P waves that show clear firstorder similarities compared with the P-wave field data.
2. Simplified FD modelling does not explain the small signal-to-noise ratio of the shear-wave field data.
For future analyses, we therefore suggest to consider additional complexity in the subsurface model that will presumably be able to explain the different quality of compressional and shear-wave field data.The most important additional factors are intrinsic damping, thin layers within the modelled units, a complex near-surface weathering layer structure, and heterogeneous material within the layers.While 2-D calculations gain faster results and allow testing the effect of different features, the full complexity of field acquisition may be understood using 3-D simulations in the future.

Figure 1 .
Figure 1.Overview maps with location of the island of Föhr (a, b).Detail map of the seismic locations on the island of Föhr (c).The profile discussed in this paper is labelled.

Figure 3 .
Figure 3. Seismic recordings of five single P-wave shots at different locations along the profile.The amplitude is displayed with an automatic gain control (AGC) of 150 ms.

Figure 4 .
Figure 4. Seismic recordings of five single shear wave shots with spatial divergence correction and AGC of 300 ms applied.

Figure 5 .
Figure 5. Final time migrated depth section of the P-wave seismic survey with AGC of 200 ms applied before time-to-depth conversion.The blue box marks the position of the SH-wave section (Fig. 6).

Figure 7 .
Figure7.Input model of P-wave velocity.The P-wave field depth section is superimposed.Numbers 1-14 mark the units listed in Table2.

Figure 8 .
Figure 8. v P /v S cross-plot from VSP data colour-coded for different lithologies.Additionally, median values, mean values and the linear regression (LinReg) are indicated.The shear-wave velocity was calculated for each P-wave velocity with the relation resulting from the linear regression.

Figure 9 .
Figure 9. Westernmost model segment of the input model.Individual single shots were modelled and the yellow stars mark the locations of the shot points (every 4 m).Black triangles mark the receiver position (each 1 m).The region of the absorbing frame is faded out.

Figure 10 .
Figure 10.Five different synthetic P-wave single shots.Shot gathers are displayed with 150 ms AGC.

Figure 11 .Figure 12 .
Figure 11.FD modelling of the P-wave section (a) and section from field data (b).Both sections are displayed with AGC of 100 ms.

Figure 13 .
Figure 13.Shear wave FD-modelling section (a) and section from field data (b).Both sections are displayed with AGC of 300 ms applied.

Table 1 .
Overview of processing of P-wave and SH-wave seismic field data including normal move-out (NMO) correction, dip move-out (DMO) correction, and common midpoint (CMP) stacking (rms -root mean square).