Characterization of a complex near-surface structure using well logging and passive seismic measurements

We combine geophysical well logging and passive seismic measurements to characterize the near-surface geology of an area located in Hontomin, Burgos (Spain). This area has some near-surface challenges for a geophysical study. The irregular topography is characterized by limestone outcrops and unconsolidated sediments areas. Additionally, the near-surface geology includes an upper layer of pure limestones overlying marly limestones and marls (Upper Cretaceous). These materials lie on top of Low Cretaceous siliciclastic sediments (sandstones, clays, gravels). In any case, a layer with reduced velocity is expected. The geophysical data sets used in this study include sonic and gamma-ray logs at two boreholes and passive seismic measurements: three arrays and 224 seismic stations for applying the horizontal-to-vertical amplitude spectra ratio method (H/V). Well-logging data define two significant changes in the P -wave-velocity log within the Upper Cretaceous layer and one more at the Upper to Lower Cretaceous contact. This technique has also been used for refining the geological interpretation. The passive seismic measurements provide a map of sediment thickness with a maximum of around 40 m and shear-wave velocity profiles from the array technique. A comparison between seismic velocity coming from well logging and array measurements defines the resolution limits of the passive seismic techniques and helps it to be interpreted. This study shows how these low-cost techniques can provide useful information about near-surface complexity that could be used for designing a geophysical field survey or for seismic processing steps such as statics or imaging.


Passive seismic (Array and H/V techniques)
The analysis of seismic noise is a valuable tool for subsurface characterization as shown by many authors such as Aki (1957) or Okada (2003). These techniques are founded on using the energetically dominant part of seismic noise which consists mainly of surface waves (Bonnefoy-Claudet et al., 2006). Two different types of methods exploit surface wave energy from seismic noise. 5 First, array techniques focus on the surface wave dispersion (frequency-dependent velocity) in order to obtain a shear-wave velocity profile by inversion process. These techniques require acquiring seismic noise simultaneously on a group of seismometers. The dispersion character of surface waves can be extracted using different methods. In this work, we focus on Frequency-Wavenumber (FK) method and Spatial Autocorrelation (SPAC) method. 10

Array data analysis with FK method (Horike (1985))
This method is based on the assumption that waves arrive in a plane across the receiver setup. The first step in the analysis of the signals according to the FK method is to obtain the so-called beam power plot. This is calculated after shifting the signal at each receiver according to a specific wavenumber (kx, ky), velocity and frequency. The stacking of the shifted signals is done in the frequency domain. Repeating this process for a wavenumber, velocity and frequency values within a defined 15 range, a complete beam power plot is retrieved. Maxima search methods of this plot provide an estimation of the travel velocity and direction of the waves.

Array data analysis with SPAC method (Aki (1957))
Another method to analyze array signals is the spatial autocorrelation method (SPAC) that assumes a random distribution of 20 seismic sources both in space and time. Aki (1957) shown that the autocorrelation ratios between two receivers is dependent on the phase velocity and the array geometry. The application of this method uses the mean of autocorrelation ratios ̅ obtained at each pair of receivers located at a distance r and considering all azimuths. Where J o is the zero-order Bessel function, c(ω) is the phase velocity for a certain frequency. Aki (1957) proposed circled 25 array configuration with different radii to obtain c(ω). Bettig et al. (2001) introduced a modification of the SPAC method that allows applying the method for different array configurations. With this modified SPAC method, the autocorrelation coefficients are obtained from station pairs separated a distance r within a certain range (r 1 , r 2 ) or rings instead of a fixed distance. Since this method is suitable for different configuration, the same geometry can be used to apply both FK and SPAC methods . 30 Solid Earth Discuss., doi:10.5194/se-2016-19, 2016 Manuscript under review for journal Solid Earth Published: 4 February 2016 c Author(s) 2016. CC-BY 3.0 License.
Dispersion or autocorrelation curves are subsequently inverted to obtain shear-wave velocity profile. Since we are dealing with 1D method, this information is assigned to the center of the array setup.

H/V technique
In addition to the array techniques, a second way of analyzing seismic noise is the horizontal-to-vertical spectral (H/V) ratio 5 method. The H/V method computes the ratio between the Fourier amplitude spectra of the horizontal and vertical components of seismic noise measurements at a single station. This technique is based on the idea that frequency-dependent ellipticity of surface wave motion can explain the H/V spectral ratio shape. In areas characterized by sediment over hardrock, it is widely accepted the association between the frequency corresponding to the H/V spectral ratio peak and the soil resonance frequency (Nogoshi and Igarashi (1970)). This technique was revised by Nakamura (1989) and has been proposed 10 as a quick, reliable and low-cost technique for site-effect characterization in earthquake seismology (Lermo and Chávez-García 1993;Bard 1999). Since the 90's, several authors have introduced the H/V method as suitable for exploration studies (e.g. Ibs-von Seht and Wohlenberg 1999;Delgado et al., 2000;Benjumea et al. 2011). These studies benefit from the relationship between the soil resonance frequency ( 0 ) and bedrock depth (H) to delineate bedrock geometry on basins (Gabàs et al., 2014). A relationship between these two quantities ( 0 and H) includes the average shear-wave velocity of the 15 sediments (V s ): 4 Data acquisition and processing

Geophysical well logging 20
During February 2012 geophysical logs were acquired in two boreholes in the study area: GW-1 up to 400 m depth and GW-3 that reached 150 m depth (Fig. 1). Two probes were used: dual induction with natural gamma sensor and three-receiver sonic probe. The logging equipment is from Robertson Geologging Ltd and includes a 500 m-winch. In this work we focus on total natural gamma values and P-wave and S-wave velocities obtained from the sonic dataset. The sonic probe has three receivers spaced 20 cm to record full-waveform data generated by a monopole source. Measurements for the two sondes 25 were made in the upward direction.
Data processing for both records includes depth matching and 11 point median filter to remove spikes. Full-wave sonic dataset was processed to measure the formation compressional and shear-wave velocities. To obtain P-wave velocities we Solid Earth Discuss., doi:10.5194/se-2016-19, 2016 Manuscript under review for journal Solid Earth Published: 4 February 2016 c Author(s) 2016. CC-BY 3.0 License. use a combination of manual first-arrival identification and semblance analysis (Kimball and Marzetta, 1984). The steps of this combined processing are: filtering of the signals for DC removing, obtaining semblance plot with the three filtered signals, identification of the first maximum of semblance map corresponding to the P-wave arrival, manual identification of the first arrival at the first receiver, obtaining theoretical first-arrival time for receiver two and three using the velocity calculated from semblance analysis, quality control of these theoretical arrivals: in case of difference between theoretical and 5 experimental arrival adjust first arrival of first receiver or adjust maximum selection on semblance plot. This procedure helps to select a maximum on the semblance analysis in case of multiple maxima and also helps to identify first arrivals with low signal-to-noise ratio.
Regarding shear-wave velocity, traditional acoustic log measurements use semblance analysis to retrieve this value from the 10 refracted shear wave. However, this is only possible in fast formations where shear-wave velocity is higher than P-wave velocity of the borehole fluid. If fluid velocity is higher than S-wave velocity no shear-wave arrival is detected (no critical refracted wave is generated). Some authors use Stoneley waves analysis to retrieve Vs information (e.g. Stevens and Day, 1986). In this work, the lack of low frequencies precludes the use of Stoneley waves for Vs estimation. Only some sectors of shear-wave velocity have been obtained in both boreholes corresponding to fast formations. 15

Array method
Previous geological information obtained from lithological description of oil well H2 (Alcalde et al., 2014) and P-and Swave velocity extracted from well logging in this work helps to define a seismic velocity model of the study area. This 20 model has been input to a dispersion curve modeling to plan the array geometry.
Three 2D arrays were deployed at the test site (Fig. 1). The procedure for each array consisted in recording simultaneous seismic noise at six stations forming two equilateral triangles with two different radii to a common center where a seventh sensor was located. One triangle was rotated 60 °respect to the other one in order to obtain good azimuthal coverage (Fig. 2). 25 After the first recording, the inner triangle was moved to a corresponding outer triangle keeping the same centre and a second time window was acquired. This procedure was repeated increasing the distance from the triangle vertex to the centre. Array 1 and 3 used radii of 25, 55, 100, 250, 400 and 1000 m. For Array 2, the radii were 10, 25, 55, 100, 400 and 750 m. The minimum radius for this array is smaller than for Arrays 1 and 3 since a significant thickness for unconsolidated sediments was expected. Seven Sara SL06 digitizers were connected to seven triaxial Lennartz LE-3D/5s sensors. The 30 sensors were covered by a plastic box to reduce wind noise. The sampling rate was fixed to 200 samples per second. All stations were equipped with GPS timing. To constitute our input dataset for retrieving shear-wave velocity, we can either use measured dispersion curves for FK method or extract autocorrelation curves using SPAC method. Both techniques have been tested and yield similar results.
Data processing has been carried out using Geopsy (http://geopsy.org/). 5 Figure 3 shows an example of FK histograms for each triangle combination corresponding to array 1-MT 12. Dotted and continuous black lines mark the limits for resolution and aliasing obtained from maximum wavenumber and minimum wavenumber respectively. These limits are calculated from the array response for the applied geometry (Wathelet et al., 2008). Continuous black line with error bars mark the maximum of the histogram related to the dispersion of surface waves energy of seismic noise. The combined dispersion curve corresponding to all the radii is shown in Fig. 4. Slowness slightly 10 decreases with frequency for frequencies higher than 4 Hz approximately. This relationship indicates the presence of high velocity layers overlying lower velocity material.
Regarding SPAC processing, we have used the modification proposed by Bettig et al. (2001) calculating the spatial autocorrelation curves on four rings of different radius ranges. An example of these curves for Array 1 is shown in Fig. 5. 15 An unique dispersion curve (DC) built from the DC of each configuration has been inverted using the neighbourhood algorithm that carries out an stochastic search through the model space (P-wave velocity, S-wave velocity, density and thickness layer). The algorithm was adapted and implemented by Wathalet (2008) in the DINVER software. The same software has been used for the inversion of the autocorrelation curves obtained from the SPAC method. 20 The parameterised model consists of four uniform layers over half-space at the base where velocity and thickness are constrained by certain ranges. Knowing the presence of velocity inversion helps parametrization for both DC and SPAC curves inversion. We also define a coupling between P and S-wave for each layer. The obtained dispersion and autocorrelation curves are inverted with the neighborhood algorithm generating a total of 12550 models for each technique. 25 Only the models with a misfit less than 0.9 will be considered. The misfit function is defined according to Wathelet et al. (2008). addition to 57 datasets coming from array stations have been analysed. In this case, the station included a SARA SL06 digitizer and the same seismometer than the first survey. Sampling frequency was set up to 200 Hz and record length ranges between 20 to 120 minutes.

H/V method
H/V spectral ratios were computed by dividing the noise recordings into 120 to 300 s-long windows. The resulting 224 H/V 5 spectral ratios were analysed considering the recommendations proposed by SESAME (Bard et al., 2004). observed ones in Fig. 6. Seismic properties combined with gamma ray measurements will be used in order to characterize Upper Cretaceous rocks following the next steps:

Geophysical logging and seismic velocities
 Zonation process based on two physical parameters (Vp and natural gamma). This process was made using WellCad software and identifies depth intervals with similar characteristics (Davis, 2002).
 Mean values calculation at each interval and for each parameter. 5  Crossplots of the mean values and visual cluster recognition. Vp and natural gamma coming from GW-1 and GW-3 are merged together in a single plot (Fig. 7a). The criterion for clustering has been identifying lowest gamma values and highest P-wave velocity values. The grey ellipsoid in Fig. 7a delineates this cluster interpreted as pure or massive limestone. Fast velocities (>4000 m/s) are reached if the clay content is below 5% (Eberli et al., 2003). The rest of the points are related to marly limestones and marls. 10  Projection of these two clusters to the corresponding depth interval in the well log plot (pure limestones in blue and marly limestones and marls in purple).
On the other hand, lower Cretaceous sediments are characterized by alternating gravel/sand and clays. Hence, seismic velocity changes are not suitable for distinguishing between layers since velocity ranges overlap for these lithological types. 15 This fact has been confirmed by the uniformity of P-wave velocity values within this sector (Fig. 6). However, natural gamma log can help to discriminate between gravel/sands and clays (Fig. 7b). Having established three ranges of gamma values based on the histogram, the projection to the depth column of these groups has been made (Fig. 6). The ranges of gamma values have been established as: 30-55 cps (gravel/sands-yellow colour), 55-85 cps (clays/gravels mix-orange colour) and 90-120 cps (clays-red colour). 20 Using the lithological description of both boreholes and based on the aforementioned geophysical parameters, we have delineated the passage from the Upper Cretaceous materials rock to the Lower Cretaceous sediments at 180 m depth at GW1 location. (iv) The shear-wave velocity of the last layer shows different values depending on the applied array technique. It ranges from 2300 m/s to 3000 m/s. The decrease of resolution of the method makes higher the uncertainty of shear-wave velocity at this depth. However, we can expect fast formation at depths higher than 600-700 m. 15 Velocity inversion stands out from these models, in particular between layer 2 and layer 3.
The first layer could be interpreted as Quaternary/Tertiary sediments. The second one could be related with the Upper Cretaceous over the third one interpreted as Utrillas and probably top of Escucha sediments. Finally the last layer can be interpreted as the bottom of Escucha Formation or Weld with Purbeck facies. The shallow part of the velocity model 20 including the velocity inversion will be compared with the logging results in the Discussion section.

H/V-frequencies and sediment thickness map
We have obtained the H/V spectral ratio for the 224 stations. The shape of these ratios can be classified in four different types.
(i) H/V spectral ratios show clear peaks that can be related to the soil resonance frequency at that station (Fig. 9a). 25 (ii) Multiple peaks can be identified. This H/V morphology is associated to multiples impedance contrasts and/or 2D/3D bedrock structure. In both cases, it is difficult to assign a soil resonance frequency although adjacent stations or other geophysical results can help to define that frequency (Fig. 9b).
(iii) Flat H/V spectral ratio indicates that the station is located on bedrock (Figure 9c).
(iv) H/V spectral ratio is characterized by a constant amplification over a wide frequency range, mainly at the high-30 frequency range (> 1 Hz), forming a step pattern. This constant amplification is not related to subsoil structure but probably to anthropic seismic noise and not to soil response (Figure 9d). For this reason, we do not use the results corresponding to this type of H/V shape. H/V frequency at each station corresponds to the frequency at the H/V spectral ratio maximum. Figure 10 shows the frequency values color coded on a map view from group (i). It also includes the location of the sites identified as rock sites (flat H/V spectral ratios-group ii) and the sites corresponding to an H/V step pattern. 5 We distinguish three zones with clear H/V peaks: two sectors in the NE and SE quadrants with fundamental frequencies ranging between 2 and 8 Hz and a third sector in the South part with frequencies values between 4 and 10 Hz. In all the sectors, the fundamental frequency gradually decreases with distance to the zone center. These sectors are related to soft soils with a significant thickness and impedance contrast with bedrock. Shallow bedrock or rock outcrops are located in the NW quadrant, SE end and central part of the study area. 10 In this study, GW3 location has been used as ground-truthing for H/V peak interpretation, i.e., to find which lithological can cause this fundamental frequency in the study area. One H/V station was located next to GW3 with the result of a fundamental frequency of 4.3 Hz. According to the GW3 lithological interpretation, a first layer of 20 m thickness composed by quaternary and Miocene sediments (clays and silts) overlies Upper Cretaceous limestone. This contact causes a 15 significant impedance contrast that produces the observed fundamental frequency of the soil. On the other hand, we can use this value to estimate a shear-wave velocity for the unconsolidated sediments that helps to produce sediment thickness maps.
Using equation (1), a shear-wave velocity value of 344 m/s is obtained. As presented in the array results, the Vs range of unconsolidated sediments is 290 m/s -500 m/s. The value obtained at GW3 applying H/V technique is included within this range. Hence, Vs varies along the study area as expected from the presence of Quaternary and Tertiary sediments of different 20 lithology. We have used the limits of the Vs range from arrays measurements to estimate a range for sediment thickness (Fig.   11). For the production of this map, the natural neighbour gridding method has been used. This method is suitable for datasets with high density in some areas and paucity in others. Maximum thickness for 500 m/s shear-wave velocity varies between 40 and 45 m ( Figure 11a) and for 290 m/s varies from 20 to 25 m (Fig. 11b).

Discussion 25
One of the limits of vintage well logging data is the dead zone in the sonic logs corresponding to the near surface. In the study area, three oil exploration boreholes lack the first 400 m of sonic data (Alcalde et al., 2014). In this framework, slim hole well logging is key to get information of the shallow part. Both Vp and Vs has been calculated for the shallow Upper Cretaceous limestone (first 90 m depth) and a complete P-wave velocity log is retrieved from the sonic log for the first 400 30 m. The velocity information increases understanding of physical properties of Upper Cretaceous materials and Utrillas (Lower Cretaceous). It will also help to interpret velocity information obtained with the array technique. On the other hand, Solid Earth Discuss., doi:10.5194/se-2016-19, 2016 Manuscript under review for journal Solid Earth Published: 4 February 2016 c Author(s) 2016. CC-BY 3.0 License. the combination with gamma ray log has enabled us to refine the geological interpretation of this shallow part. This interpretation constrained the results for the other geophysical results obtained in this study.
Regarding the passive seismic methods, the H/V technique can provide a detailed sediment thickness map also required for resolving near-surface issues. When borehole information is sparse, the H/V technique can be a good alternative for 5 obtaining a shallow subsurface model for statics calculation or field survey planning.
The array technique results helped to identify a significant velocity inversion in the first 100 m depth. The bottom of the high-velocity layer has been identified between 60 and 80 m approximately depending on the array location and method of surface waves analysis (FK or SPAC). If we compare these results with the shear-wave velocity log obtained at GW1 and 10 GW3 (Fig. 12), we observe that the bottom of the high-velocity layer would correspond to the change from a sector characterized by velocity decrease with depth to another one made by interbedded layers of high and low velocity. The influence of this type of layering in a seismic signal with a wavelength longer than the fine-scale details of velocity variations has been studied by numerous authors (e.g . Stovas and Ursin, 2007). According to Hovem (1995), for a large ratio between wavelength (λ) and thickness (d) of one cycle in the layering, the layered sector behaves as a homogeneous, 15 effective medium (Backus 1962). For low ratio, the alternate layers may be replaced by thicker single layer using a timeaverage or ray velocity for the total depth range.
In our case, seismic wavelength of the surface waves analyzed with the array technique varies from 125 m at 50 m to 500 m at 200 m according to the dispersion curve (Fig. 4). On the other hand, spatial Fourier analysis of the GW1 sonic log from 92 20 to 200 m gives two maximum at 10 and 18 m. In our case, the λ/d ratio would range between 5 at 92 m to 25 at 200 m. Since these values would correspond to different models (Backus average or time-average velocity), we have considered both models.
Firstly, a P-wave time-average velocity has been obtained from 92 to 200 m giving a value of 2348 m/s. This is in agreement 25 with a velocity inversion between the near-surface limestone and the layers below (marly limestones and marls). Since Swave velocity is not available for this sector we cannot obtain an S-wave time average velocity.
Secondly, from a Backus average point of view, we need first to define the number of layers and assign a P-and S-wave velocity characteristic of the high and low velocity layers. For a characteristic thickness of 10 m, we can use a number of 10 30 layers with alternate high and low velocity as equivalent to the interbedded layer sector (92 to 200 m). P-wave for the high velocity layers has been fixed to 4000 m/s and for the low-velocity layers to 1500 m/s. This leads to a Backus P-wave velocity of 1986 m/s. Regarding S-wave velocity, we have chosen a high-velocity of 2000 m/s confirmed by Vs calculation from sonic log where some hints of shear-wave velocity were obtained at the high-velocity layers (e.g. at around 150 m When we compare the S-wave velocity values obtained from well logging and the array technique, we observe that array velocity is always lower than S-wave logging velocity. This result is consistent with the fact that sonic logging uses highfrequency signals which travels at different velocity than low-frequency seismic signals (Box and Lowrey, 2003). Due to the 15 dispersion effects, the pulses generated with a sonic probe travel a few percent faster than the array surface waves.

Conclusions
Slim well logging completes the depth record of old well logging data with sonic data from the surface down to 400 m depth. 20 This information combined with the gamma ray log helps to refine geological interpretation of the first 400 m. P-wave velocity data delineates two significant changes with depth in the Upper Cretaceous materials: a high velocity layer overlying a zone with velocity decrease with depth, and an interbedded sector with high and low velocities. For the propagation of a seismic wave, this sector is equivalent to a low velocity layer. Below this layer, we find low-velocity sediments corresponding to Utrillas (Low Cretaceous) formation. 25 H/V technique enabled us to obtain a sediment thickness map of a complex near-surface area in a fast and affordable way.
This methodology can produce good results in areas with scarcity of boreholes.
Shear-wave velocity profile up to 800 m is obtained with the array technique. Regarding the near-surface, this profile maps a 30 velocity inversion within the Upper Cretaceous that has been correlated to the presence of the high-velocity layer over the interbedded sector. The change from Upper to Lower Cretaceous sediments has not been delineated with this technique. Solid Earth Discuss., doi:10.5194/se-2016-19, 2016 Manuscript under review for journal Solid Earth Published: 4 February 2016 c Author(s) 2016. CC-BY 3.0 License.
As a first step, passive seismic methods combined with near-surface geophysical logging constitute a fast way of characterizing near-surface complexity. Identification of velocity inversions and sediment thickness map are valuable inputs for geophysics field planning, seismic statics correction or further modelling. 5 Solid Earth Discuss., doi:10.5194/se-2016-19, 2016 Manuscript under review for journal Solid Earth Published: 4 February 2016 c Author(s) 2016. CC-BY 3.0 License.

Author contribution
S. Figueras and A. Gabàs planned the field surveys for well logging and acquired the data from GW3 borehole. They have also participated with comments and ideas for processing and interpretation. A. Macau planned the field surveys, acquired well logging data from GW1 borehole, the passive seismic datasets and processed the H/V and array data. He carried out the 5 inversion of the dispersion and autocorrelation curves and interpretation of the shear-wave velocity models. B. Benjumea planned the field surveys, acquired well logging data from GW1 borehole and passive seismic datasets. She analyzed the well logging data and processed the sonic datasets. She has interpreted the datasets and made the comparison between welllogging and passive seismic results. She prepared the manuscript with contributions from all co-authors.

Acknowledgments 10
We would like to thank Andrés Pérez-Estaún in memoriam for being the person who invited us to participate in this project.
He contributed with enthusiasm and wisdom to the interpretation of well logging data (Appendix A).   wavenumber value of kmin/2 related to the resolution limit of this array geometry (Wathelet et al., 2008). The dashed line is a constant wavenumber value of kmax that defines the aliasing limit constrained by the array geometry. The dots represent kmin (left) and kmax/2 (right) respectively. Solid Earth Discuss., doi:10.5194/se-2016-19, 2016 Manuscript under review for journal Solid Earth   Figure 6. a) P-wave velocity (blue line) and S-wave velocity (red line) for GW-1, b) Recorded sonic waveform for the farreceiver from GW-1 borehole c) natural gamma log from GW-1 and lithological interpretation resulting from zonation process and clustering as shown in Figure 7 for GW-1. e), f), g) and h) same as a), b), c) and d) respectively but for GW-3.