# Site characterization at Groningen gas field area through joint surface-borehole H/V analysis

Site characterization at Groningen gas field area through joint surface-borehole H/V analysis SUMMARY A new interpretation of the horizontal to vertical (H/V) spectral ratio in terms of the Diffuse Field Assumption (DFA) has fuelled a resurgence of interest in that approach. The DFA links H/V measurements to Green’s function retrieval through autocorrelation of the ambient seismic field. This naturally allows for estimation of layered velocity structure. In this contribution, we further explore the potential of H/V analysis. Our study is facilitated by a distributed array of surface and co-located borehole stations deployed at multiple depths, and by detailed prior information on velocity structure that is available due to development of the Groningen gas field. We use the vertical distribution of H/V spectra recorded at discrete depths inside boreholes to obtain shear wave velocity models of the shallow subsurface. We combine both joint H/V inversion and borehole interferometry to reduce the non-uniqueness of the problem and to allow faster convergence towards a reliable velocity model. The good agreement between our results and velocity models from an independent study validates the methodology, demonstrates the power of the method, but more importantly provides further constraints on the shallow velocity structure, which is an essential component of integrated hazard assessment in the area. Downhole methods, Interferometry, joint inversion, seismic noise, site effects 1 INTRODUCTION Due to long-term exploitation of a large onshore gas field and subsequent compaction of the reservoir at depth, the Groningen area in the northern Netherlands is subject to ground subsidence and induced earthquakes. To date, the largest induced earthquake occurred near Huizinge in August 2012 and was recorded with a local magnitude ML of 3.6 (moment magnitude M = 3.4). This event revived concerns about the hazard from induced seismicity and prompted the government and the Nederlandse Aardolie Maatschappij (NAM) to take action. As part of the response, NAM built a dense permanent borehole microseismic network covering an area of about 35 × 45 km2. The network (Fig. 1) is composed of ∼70 accelerographs at the surface that are co-located with a 200 m depth borehole in which geophones are installed at depth intervals of 50 m. The network is fully operational and has recorded continuous waveforms since late 2015. Figure 1. View largeDownload slide (A) KNMI-NAM permanent borehole sites (the G-array) in the Groningen gas field area. The inverted coloured triangles refer to the H/Vs shown in Fig 2. The blue contour depicts the outline of the gas field. White areas depict urban centres. (B) The location of the Groningen gas field (red box) in the northern part of the Netherlands. Figure 1. View largeDownload slide (A) KNMI-NAM permanent borehole sites (the G-array) in the Groningen gas field area. The inverted coloured triangles refer to the H/Vs shown in Fig 2. The blue contour depicts the outline of the gas field. White areas depict urban centres. (B) The location of the Groningen gas field (red box) in the northern part of the Netherlands. It is well known that near-surface lithology can strongly influence damage from earthquake shaking by increasing the amplitude and duration of shaking, and by responding nonlinearly to incident seismic waves (e.g. Olsen 2000). For this reason, site characterization is of great importance to seismic hazard analysis (Bommer et al. 2017). Such characterization requires a knowledge of geomechanical properties of the stratigraphy at a site, which can be often challenging to obtain. Over the last few decades, ambient seismic field seismology has emerged as a valuable tool to characterize shear wave (VS) velocity models over all distance scales. For site characterization, these techniques include the inversion of dispersion curves obtained by ambient seismic field correlation using small-aperture arrays such as the spatial autocorrelation (SPAC; Aki 1957) or frequency–wavenumber methods (f-k; Capon 1969; Lacoss et al. 1969). Depending on the aperture of the array, these techniques can recover the local 1-D velocity structure from surface to a depth of up to several hundred metres. Ambient seismic field cross-correlation tomography with high-density arrays is an emerging method to obtain high-resolution 3-D models of the shallow surface (e.g. Nakata et al. 2015; Roux et al. 2016); however, these kinds of experiments are very expensive, and require significant investment for design, deployment, and operation. The ratio of the horizontal to vertical (H/V) components of the ambient seismic field (sometimes referred to as microtremor) is a widely used method to determine a simple velocity model of the subsurface (i.e. one layer over a half-space) and from that the frequency-dependent site response. The reason this method gives reliable results, however, has long been controversial. This was due to the absence of a clear theoretical basis for the measurements, which leads, inevitably, to lack of clarity in its interpretation. It is not even clear which waves comprise the noise field that generates the H/V peak frequencies (Nakamura 2000; Malischewsky & Scherbaum 2004; Bonnefoy-Claudet et al. 2008). The possibility to retrieve the elastodynamic Green’s function between two stations embedded in an elastic medium from the average time-domain cross-correlation of ambient field records (Weaver 2005) led Sánchez-Sesma et al. (2011a) to advance a new theory for H/V. Based on the diffuse field assumption (DFA), Sánchez-Sesma et al. (2011a) linked the H/V of the autocorrelated signal to the ratio of the imaginary parts of the Green’s functions. The observed H/V can be compared to its theoretical counterpart, and provides a basis for estimating the subsurface structure. It represents an opportunity to obtain the local subsurface velocity model with only one short three-component measurement of the ambient field (e.g. only 20 min in Spica et al. (2015)). This capability is of particular interest in geotechnical engineering, seismic exploration, and engineering seismology. To date it has only been applied to limited data sets and geological settings (Salinas et al. 2014; Spica et al. 2015; Lontsi et al. 2015; Piña-Flores 2015; Rivet et al. 2015; Lontsi et al. 2016; García-Jerez et al. 2016; Piña-Flores et al. 2016; Perton et al. 2017). While these initial results are promising, there is still a need to explore the potential of the method and to validate the results against independent information. In this contribution, we pursue the idea of Lontsi et al. (2015) to use multiple H/V measurements obtained from receivers at depth in a borehole to estimate complex geological structure at a site. The Groningen area, is a well-studied area and it has ∼70 borehole sites, such that it represents an ideal natural laboratory to test the method. We start with an overview of the geological horizons of interest that are sampled by the measurements and some background on constraints on structure from previous studies (Section 2). We then outline the theoretical background and seismic processing necessary to compute H/V (Section 4), and the discrete wave number (DWN; Bouchon 1981) method used to compute the H/V ratios for the forward problem in Section 4.2. In Section 5, we describe the estimation procedure for the velocity models at several borehole sites. Finally, we discuss results and compare them with velocity models obtained independently (Section 6). In this paper we only present a few inversions because our focus is primarily on implementation of the method and the validation of the results. 2 GEOLOGICAL SETTING AND EXISTING VELOCITY MODEL The Groningen area is characterized by flat, low-lying topography with altitude close to mean sea level. The horizons of interest (i.e. those sampled by our measurements) are the Paleogene, Neogene and younger deposits overlying the North Sea Supergroup (De Mulder et al. 2003; Vos 2015). This ∼800 m thick layer of unconsolidated sediments contains a large degree of vertical and lateral heterogeneity due to the influence of the last three ice ages and associated sea level fluctuations. The deposits range from fluvial braid plain sands to shallow marine (intertidal) and terrestrial deposits of soft clays to distinct organic-rich peat formations. The uppermost sedimentary sequence is characterized by a succession of fluvial, glacial, and marine deposits that are crosscut by deep subglacial features (‘tunnel valleys’), which are filled with sands and clays and buried under younger sediments. The North Sea Supergroup consists primarily of alternating marine grey sands, sandstones, and clays. Our study benefits from important, independent information in the form of an integrated VS model from the surface to the base of the North Sea Supergroup (Kruiver et al. 2017). This 3-D VS model is a synthesis of three different VS models obtained through a comprehensive set of geological, geotechnical, and geophysical observations; each of them having a different depth range and spatial resolution. The shallowest part of the model extends from the surface to ∼50 m below sea level, and is based on a high-resolution 3-D geological model GeoTOP (Stafleu et al. 2011; Stafleu & Dubelaar 2016), combined with VS distributions with depth for the sediments (Kruiver et al. 2017). At depths of ∼40 m to ∼120 m below sea level, the model is constrained through inversion of Rayleigh wave observations from spatially extensive reflection seismic surveys. The deepest part of the model is based on the Pre-Stack Depth Migration velocity model derived from sonic logs. It ranges from ∼70 m below sea level to the base of the North Sea Supergroup at (∼800 m). These three distinct models were spliced into a single model that covers the full depth range for site response analyses. In what follows, we refer to this velocity model as the Deltares-NAM model. Other shallow velocity models were provided by Noorlandt et al. (In review). We use 1-D VS profiles of the Deltares-NAM model at the borehole sites to validate our analysis. It is important to keep in mind, however, that the resolution and sensitivity of the methods used to construct the Deltares-NAM model and the H/V velocity profiles differ, as discussed in Section 4.3. 3 DATA We used continuous data from the shallow borehole network installed in the Groningen area (Fig. 1; the G-array). Each borehole site includes one accelerograph at the surface that is collocated with four geophones at 50 m depth intervals (−50, −100, −150, −200 m). The continuous data are freely available on the Royal Dutch Meteorological Institute (KNMI) website. We downsample the continuous data to 50 sps, remove instrument response and convert the time-series to velocity prior to processing. 4 H/V FOR A DIFFUSE WAVEFIELD By definition, the H/V spectral ratio corresponds to the square root of the spectral energy ratio of the horizontal amplitudes (with indices 1, and 2) over the vertical direction (index 3; Arai & Tokimatsu 2004) :   $$\frac{\text{H}}{\text{V}}(\mathbf {x},\omega )=\sqrt{\frac{E_1(\mathbf {x},\omega )+E_2(\mathbf {x},\omega )}{E_3(\mathbf {x},\omega )}}.$$ (1)Perton et al. (2009) showed that at an observer location (surface or depth), these spectral energies (i.e. the directional energy densities, as defined in Perton et al. (2009)) are proportional to the average autocorrelations of the diffuse wavefield components, which in turn are proportional to the imaginary parts of the Green’s function ($$Im[\mathcal {G}]$$) (Sánchez-Sesma et al. 2011b):   $$E_i(\mathbf {x},\omega ) = \rho \omega ^2\langle u_i(\mathbf {x},\omega )u_i^*(\mathbf {x},\omega )\rangle \propto -\omega \mbox{Im} [ \mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega )].$$ (2)where ω is the angular frequency, ui(x, ω) is the displacement field in the i direction at a point x, Im[ ] stands for the imaginary part, $$\mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega )$$ is the displacement Green’s function in the direction i at a point x due to the application of a unit point force in the same direction applied at the same point. The symbol * stands for the complex conjugate operator and the product $$u_i(\mathbf {x},\omega )u_i^*(\mathbf {x},\omega )$$ corresponds in the frequency domain to the autocorrelation of the displacement field in the time domain. The brackets 〈 〉 represent averaging over time. In what follows, the dependence on ω and x is implicit. The ∝ means that the expressions are proportional by a factor that is independent of ω and x. Eq. (2) represents the same case as classic ambient seismic field cross-correlations (e.g. Weaver & Lobkis 2004) but for the special case when the source and receiver are the same. Within this theoretical framework, Sánchez-Sesma et al. (2011a) proposed a theoretical description of H/V ratios and suggested that the H/V spectral ratio recorded at a receiver could also be computed in terms of the imaginary part of the GF:   \begin{eqnarray} \displaystyle\frac{\text{H}}{\text{V}}(\mathbf {x},\omega )&=\sqrt{\displaystyle\frac{\left\langle | u_1(\mathbf {x},\omega )|^2 \right\rangle +\left\langle | u_2(\mathbf {x},\omega )|^2 \right\rangle }{\left\langle | u_3(\mathbf {x}, \omega )|^2 \right\rangle }} = \sqrt{\displaystyle\frac{ \mbox{Im}( \mathcal {G}_{11}+\mathcal {G}_{22}) }{ \mbox{Im}(\mathcal {G}_{33})}}. \end{eqnarray} (3)Eq. (3) links the average energy densities (i.e. the ambient field measurements; see Section 4.1) with the Green’s function (i.e. the theoretical counterpart; see Section 4.2) and treats the H/V spectral ratio as an intrinsic property of the medium. It naturally allows for the inversion (see Section 5) of H/V that includes contributions of the full wavefield—that is, including both surface and body waves. 4.1 Observed H/V In the context of the DFA, eq. (2) is only valid when the seismic wave field is equipartitioned, that is, all the incident waves are (e.g. P, S, or Rayleigh waves) have the same energies (Perton et al. 2016). This assumption is difficult to verify and unlikely to be true in the majority of ambient field data. Therefore, some signal processing must be applied to enhance the equipartitioning of the seismic wavefield, just as for traditional ambient seismic field cross-correlation (e.g. Bensen et al. 2007). As in Spica et al. (2015) and Perton et al. (2017), we apply spectral whitening, which corresponds to source deconvolution. Because several sources can act in different frequency bands, the operation consists of normalizing the signals by the source energies computed in each time window and across several frequency bands as $$\tilde{v}_i(\mathbf {x},\omega )=v_i(\mathbf {x},\omega )/\sqrt{\sum _{i=1}^3|v_i(\mathbf {x},\Delta \omega )|^2}$$; where Δω is a frequency band of 2.5 Hz width centred on ω. We work with the particle velocity vj(x, ω) = iωuj(x, ω) to preserve the link to energy. To remove only the spectral envelope, the bandwidth has to be much larger than the oscillations in the spectra. The H/V spectral ratio is therefore computed in terms of the wavefield autocorrelations as:   \begin{eqnarray} \frac{\text{H}}{\text{V}}(\mathbf {x},\omega )&=\sqrt{\displaystyle\frac{\left\langle | v_1(\mathbf {x},\omega )|^2 \right\rangle +\left\langle | v_2(\mathbf {x},\omega )|^2 \right\rangle }{\left\langle | v_3(\mathbf {x},\omega )|^2 \right\rangle }}. \end{eqnarray} (4)Eq. (4) requires that the averaging is performed separately for each component. In that sense, it is different from the calculation of the usual H/V spectral ratio (Nakamura 1989), which corresponds to the average of the ratios: H/V = 〈Hw/Vw〉. The average is computed on the spectra obtained over one day of continuous data that is windowed into sections of 100 s duration with an overlap of 20 per cent. Each time window is demeaned, detrended and bandpass filtered from 0.1 to 10 Hz. The attenuation is generally a severe limitation for a field to become diffuse and then equipartitioned. Several studies have shown this limitation for cross-correlation with large interstation distance (e.g. Lawrence & Prieto 2011). However, here, the interstation distance is null and the attenuation has a similar effect on waves in all directions. Provided there is sufficient nearby sources of noise, attenuation should not be a limitation. We have also shown theoretically that in presence of attenuation the autocorrelation is proportional to the $$Im[\mathcal {G}]$$ for the same media but without attenuation (Perton & Sánchez-Sesma 2014). A possible limitation of the method is, however, the presence of strong anisotropy in the shallow sediment. This is not taken into account in this contribution. Both the evaluation of anisotropy in the area and its assessment through H/V spectral ratio are possible future research directions. Fig. 2 shows different H/V curves computed for the NW–SE line of boreholes highlighted in Fig. 1. By far, the strongest variations in the H/V curves are observed for the surface measurements (black lines in Fig. 2). We attribute the higher amplitude of the surface H/V ratios to a stronger impedance contrast the air–solid interface and a stronger contribution from surface waves. The H/V curves at depth have much lower amplitude variations, but slight changes in their shape reflect structural changes at depth. These are expressions of shallow structural variability from northwest to the southeast in the study area. Figure 2. View largeDownload slide H/V at different depth levels computed along the NW–SE line of stations shown as inverted coloured triangles in Fig. 1. The dashed grey lines in each panel is the $$\sqrt{2}$$ level (see Section 4.3). The H/V in the green frames (G03, G23, G29 and G56) are shown in Fig. 5 along with their inversion results. Figure 2. View largeDownload slide H/V at different depth levels computed along the NW–SE line of stations shown as inverted coloured triangles in Fig. 1. The dashed grey lines in each panel is the $$\sqrt{2}$$ level (see Section 4.3). The H/V in the green frames (G03, G23, G29 and G56) are shown in Fig. 5 along with their inversion results. In order to test that the observed H/V is not biased due to lack of diffusivity of the ambient field, we analyse the variability between $$\frac{\text{H}_{\text{E-W}}}{\text{V}}$$ and $$\frac{\text{H}_{\text{N-S}}}{\text{V}}$$ for stations at the surface (Supporting Information Fig. S1). The fact that these measurements are similar at almost all the sites, suggests a near isotropic illumination by the noise sources (Perton et al. 2017). Equipartition of the wavefield is also discussed in Section 4.3. 4.2 Theoretical H/V In the second term of eq. (3), the $$\text{Im}[\mathcal {G}]$$ components must be related to a geometry and material properties that explain the data. Perton et al. (2017) explored the case of strong lateral heterogeneities along a crater cliff; however, the case of Groningen does not have strong topographic variation of the surface, which allows us to use a 3-D unbounded layered elastic media and smooth variation of the subsoil in the vicinity of the stations as a reasonable starting assumption. Several methods have been proposed to model the $$\mathcal {G}_{ii}$$ components of H/V in a layered medium under the DFA (Sánchez-Sesma et al. 2011b; García-Jerez et al. 2016; Perton & Sánchez-Sesma 2016; Lontsi et al. 2015). Here, we use the Discrete Wave Number method (DWN; Bouchon 2003, 1981). The DWN method is efficient and suitable for solving the forward problem for H/V for stations located at the surface (Sánchez-Sesma et al. 2011a; Perton et al. 2017). In what follows, we discuss its efficiency when several stations are located in a single geometrical and material configuration at different vertical and/or horizontal positions. As in Sánchez-Sesma et al. (2011a), the wavefield is decomposed as a sum of plane waves according to the horizontal component of the wavenumber vector in the radial direction kr. The Green’s function is the sum of two contributions, one due to the P and SV waves and another due to SH waves. A compact form of the displacement P–SV Green’s function component $$\mathcal {G}^n_{ii}(\mathbf {x})$$ inside the layer n, where receiver and point source are superimposed at x = {x, z} and where the vector source is oriented in the same direction i is:   \begin{eqnarray} \mathcal {G}^n_{ii}(\mathbf {x}) &=& \mathcal {G}^{\text{FS}n}_{ii}(\mathbf {x}) \nonumber\\ && +\, \int _{-\infty }^{+\infty } \big( S^{\nearrow Pn} u_i^{\nearrow Pn} + S^{\nearrow SVn} u_i^{\nearrow SVn} +S^{\searrow Pn} u_i^{\searrow Pn} \nonumber\\ && +\, S^{\searrow SVn} u_i^{\searrow SVn} \big) \text{d}k_r \end{eqnarray} (5)For an extended description of the integrand, the interested reader can consult the textbook of Aki & Richards (2002). $$\mathcal {G}^{\text{FS}n}$$ is the incident field as in a full space (FS) having the same properties as the layer n. The wave amplitudes S are unknowns that are solved for at discrete values of kr by enforcing the continuity of displacement and traction at each interface. Superscripts P and SV are for P and SV waves. We write the system of equations in matrix form [A][S] = [B]; where [S] = {..., S↗Pn, S↘Pn, S↗SVn, S↘SVn, ...} is the unknown vector and [B] is the vector associated with the stress and displacement components of the source (i.e. to $$\mathcal {G}^{\text{FS}n}_{ii}$$). [A] is comprised of the amplitude coefficients of the boundary conditions and has dimensions equal to [4(N − 1))*(4(N − 1)] with N being the number of layers. As several sources are considered at different depths and for different directions, the system of equations should be solved several times (i.e. for each source); however, since we consider a single configuration (i.e. the matrix [A] is the same), it is more efficient to calculate [A]−1 just once. Therefore, the wave amplitudes related to each source are calculated as the result of the multiplication [S] = [A]−1[B], rather through inversion by Gaussian elimination for each vector [B]. In this way all the Green’s function’ components $$\mathcal {G}^n_{ii}$$ can be efficiently evaluated at multiple receiver positions and for the several directions. The same procedure applies for the Green’s function associated with the SH contribution. This allows consideration of all H/V for the same borehole site simultaneously during the joint inversion. 4.3 Forward modelling using a complex velocity model We show in Fig. 3 the Deltares-NAM VS model at borehole site G03. In the first 50 m of the model, there are two strong velocity contrasts (C1 and C2 in Fig. 3). The deeper velocity dependence is mainly expressed by a velocity gradient (dashed black line in Fig. 3). Observed (black lines) and computed H/V (blue lines) at each depth level are shown on the right part of the figure. Figure 3. View largeDownload slide Theoretical H/V computed from the Deltares-NAM model (blue lines) and observed H/V (black lines) at borehole G03. The letters C1 and C2 refer to possible strong velocity contrasts in the upper 50 m depth range, and the dashed line refers to the gradient-like part of the velocity model. Figure 3. View largeDownload slide Theoretical H/V computed from the Deltares-NAM model (blue lines) and observed H/V (black lines) at borehole G03. The letters C1 and C2 refer to possible strong velocity contrasts in the upper 50 m depth range, and the dashed line refers to the gradient-like part of the velocity model. The discrepancies between the observed and synthetic H/V suggest that the inversion of the observed H/V will result in a different velocity model, but with certain common characteristics. While the model of Kruiver et al. (2017) is site-specific and has a resolution that changes with depth (see Section 2), the H/V are sensitive to elastic properties over a larger horizontal area that depends on the wavelength of the waves that contribute to the observations (Perton et al. 2009; Piña-Flores et al. 2016). The joint inversion result will therefore give an averaged (vertical and horizontal) velocity model at a borehole site. The H/V at the surface show one clear frequency peak, while H/V at depth shows small oscillations superimposed on a flat trend. As recently detailed by Piña-Flores et al. (2016), H/V is sensitive to both surface and body waves. The H/V spectra from receivers at the free surface mainly originate from the strong contribution of surface waves and show important localized frequency peaks of high amplitude. These frequency peaks are dependent on the existence of strong impedance contrasts at depth. In the case of a simple velocity model of one layer over a half-space (with shear velocity β1 and thickness h1), the main peak frequency $$f_{\beta _1}$$ is well approximated by (Yamanaka et al. 1994; Piña-Flores et al. 2016):   $$f_{\beta _1} = \frac{\beta _1}{4h_1}.$$ (6)When two layers are considered on top of a half space, as suggested by the Deltares-NAM model (i.e. C1 and C2 in Fig. 3), H/V can present two peaks depending on the thickness of each layer and the impedance contrasts at their boundaries. If the thickness of the deepest layer is greater than the shallower layer and a sufficient impedance contrast between these three materials exists, the two main peaks are well separated (Bonnefoy-Claudet et al. 2006; Field & Jacob 1995; Bard & SESAME-team 2004; Piña-Flores et al. 2016). Therefore, if β2 ≫ β1 and h1 ≫ h2, the shear velocity of the superficial layer can also be evaluated using eq. (6) (here $$f_{\beta _1}$$ is the high frequency peak) and the approximate frequency associated with the peak at low frequencies can be evaluated from   $$f_{\beta _{12}}=\displaystyle\frac{1}{4\left(\displaystyle\frac{h_1}{\beta _1}+\frac{h_2}{\beta _2}\right)},$$ (7)where β2 and h2 are the thickness and shear velocity of the second layer. While more complex formulae that reflect several major effects of the model on the resonance frequency exist (Tuan et al. 2016), modelling of eq. (7) appears to be a good approximation for layered structure (e.g. Piña-Flores et al. 2016); however, since H/V at the free surface in Fig. 3 does not show such character (i.e. a double frequency peak), we would not expect the inverted velocity model to develop two strong velocity contrasts in the first 50 m at this site. Because surface wave energy quickly decreases with depth, H/V in boreholes becomes more sensitive to body waves. While surface waves propagate in 2-D space and are generally not strongly reflected by lateral heterogeneity, body waves propagate in 3-D space and are reflected by the free surface and also by possible strong impedance contrasts at depth. As shown in Perton et al. (2009) for a half space, the waves that travel vertically up and down interfere and result in spectral oscillation periods in the energy density components (E1, E2, E3). The amplitude of these oscillations decays with frequency and depth. In a half-space, the dependence on frequency and depth are similar. The $$H=\sqrt{E_1+E_2}$$ is mainly sensitive to the shear wave velocity while the $$V=\sqrt{E_3}$$ is mainly sensitive to the compressional wave velocity. Also, as a requirement of equipartition, the three energy densities tend to the same value with depth. Therefore, from eq. (3), H/V in a half-space should present these oscillations and tend to $$\sqrt{2}$$ at high frequency and at depth (as observed in Fig. 2). The gradient-like velocity structure suggested by the Deltares-NAM model does not modify these conclusions since the observed H/V at depth expresses these features. The latter suggests that we might be able to evaluate the velocity at each sensor depth; however, the absence of strong impedance contrast at depth reduces our ability to recover a wave reflected away from the buried station and H/V is only weakly sensitive to structure deeper than the deepest sensor (i.e. ≥200 m). Other examples of observed H/V based on the complex velocity model provided by Deltares-NAM are shown with the inversion results below. This allows us to assess how close the site-specific model is from the H/V measurements and to assess the prospects for improvement. 5 INVERSION As discussed in Piña-Flores et al. (2016), consideration of H/V at surface alone is insufficient to characterize shallow properties uniquely since velocities and thicknesses trade-off and lead to a similar H/V. Additionally, the forward problem is highly non-linear and depends on several uncorrelated parameters (Piña-Flores et al. 2016; García-Jerez et al. 2016). Constraining the inversion by adding observations from sensors at depth significantly reduces the possible range of parameters (Lontsi et al. 2015). Starting with an accurate first guess—as, for example, a guess obtained from an independent measurement such as a dispersion curve (Scherbaum et al. 2003; Piña-Flores 2015; Lontsi et al. 2016)—allows more rapid convergence. Although the Deltares-NAM model exists, our goal is to provide independent measurements in our analysis to test the validity of the method. 5.1 Starting models from borehole interferometry We use ambient field borehole interferometry between pairs of adjacent overlying sensors (e.g. Miyazawa et al. 2008) to obtain the mean shear wave traveltime. In Fig. 4(A), we show the geometry of the borehole stations and the sensor pairs used for interferometry. The cross-correlation of the ambient seismic field v(zS) and v(zr) in the frequency domain is given by:   $$C_{\text{AB}}(\omega )= \langle v(z_S)^{\ast } v(z_r) \rangle ,$$ (8)where zS and zr are the depths of the virtual source and of the receiver, respectively, which are any of the consecutive sensors. Then, the imaginary part of the Green’s function between pairs of sensors ($$\text{Im}[\mathcal {G}]_{(z_S,z_r)}$$) is obtained from the average correlation using   $${}\text{Im}[\mathcal {G}]_{(z_S,z_r)} \propto \frac{\langle v(z_S)^{\ast } v(z_r) \rangle }{|S(\omega )|^2},$$ (9)where |S(ω)|2 is the power spectrum of the noise (e.g. Wapenaar & Fokkema 2006). Figure 4. View largeDownload slide (A) Geometry of the boreholes where the grey stars represent the virtual sources (S) and triangles represent boreholes (rb) or surface receivers (rs). Ambient field cross-correlations are performed for every depth level independently. (B) Obtained correlation functions for the north–south and east–west components. The picked traveltime velocities (orange dots) are marked for every depth level and the mean velocity for the entire well is marked in red. Note that similar traveltimes are picked in both horizontal components, suggesting weak horizontal anisotropy at this site. (C) Resulting average velocity model (init) used as starting model for inversion, compared to the Deltares-NAM velocity model harmonically averaged over 50 m intervals. Figure 4. View largeDownload slide (A) Geometry of the boreholes where the grey stars represent the virtual sources (S) and triangles represent boreholes (rb) or surface receivers (rs). Ambient field cross-correlations are performed for every depth level independently. (B) Obtained correlation functions for the north–south and east–west components. The picked traveltime velocities (orange dots) are marked for every depth level and the mean velocity for the entire well is marked in red. Note that similar traveltimes are picked in both horizontal components, suggesting weak horizontal anisotropy at this site. (C) Resulting average velocity model (init) used as starting model for inversion, compared to the Deltares-NAM velocity model harmonically averaged over 50 m intervals. The cross-correlations are computed for both east-west and north-south component pairs between the wave motion observed at a depth zS and another directly overlying sensor zr. Only ambient noise is used and small earthquakes are removed (e.g. Miyazawa et al. 2008) based on the KNMI earthquake catalogue. Correlations are computed over one month of continuous data for 25 s windows. Fig. 4(b) shows the resulting bandpass filtered (1.5–8 Hz) stacked correlations in the time domain. An impulsive arrival is observed on both causal and anti-causal parts of the Green’s function, corresponding to upward and downward propagating S-waves. We estimate the arrival time as in Nakata & Snieder (2012); that is, by seeking the three adjacent samples with the largest amplitude values and applying a quadratic interpolation to find the time at which the wave has maximum amplitude. This time is taken as the traveltime for a shear wave that propagates between the borehole and overlying sensor. We use the average of estimated traveltime on the east-west and north-south components to estimate an average velocity layer. The obtained velocity model at site G43 is shown in Fig. 4(C) and is compared to the Deltares-NAM velocity model harmonically averaged over 50 m intervals. These results might be biased by multiple arrivals caused by shear wave reflection inside the 50 m section due to strong geological contrasts. Besides, if the noise is not fully diffuse, Green’s functions are not well retrieved and some small bias can be observed in the traveltimes (Tsai 2009). In such case, the interferometry picked traveltime would be smaller than the expected one, which results in higher velocities. Because the average velocity models tend to underestimate the Deltares-NAM velocity model at most of the sites (i.e. Fig. 4C and Supporting Information Fig. S2), and because the analysis performed on the H/Vs in Supporting Information Fig. S1, this is not likely to be a strong effect. The discrepancy observed in Supporting Information Fig. S2 is also partly explained by the fact that the interferometry picked traveltimes are directly translated to average velocity. Contrary to the traveltime, the velocity estimated by the interferometry should not be exactly equal to the average of the velocities in the 50 m section since it corresponds to a harmonic average. 5.2 Parametrization, misfit function and inversion The soil properties and the layer thicknesses inside a borehole site are assessed by inverting jointly H/V at several depths, z. The objective function is defined as the root mean square difference between observed and predicted H/V computed with the DWN method:   $$\epsilon = \sqrt{\frac{1}{N_z N_{\omega }}\sum _{z_i}\sum _{\omega } \Bigg ( \frac{H}{V}^{\text{exp}} - \frac{H}{V}^{\text{DWN}}\Bigg )^2; }$$ (10)where Nz = 5 is the number of depths at which zi considered in the joint inversion, and Nω = 14 is the number of points taken in the spectra from 1. to 8. Hz. Fluctuations in the H/V spectra at higher frequency are associated with small, local heterogeneities, which are of minor importance here (Piña-Flores et al. 2016). The only free parameters considered during the inversion are the shear wave velocities (β). To reduce the number of parameters, the compressional velocity (α) and the mass density are assumed to be related to the shear velocity through polynomial relationships (Brocher 2005; Berteussen 1977). We impose that one velocity layer in each section of 50 m (i.e. between each sensor pair) is constrained by the interferometry picked traveltime and by the other velocities of the section. In other words, if i is the layer index of height hi and shear velocity βi, and j0 ≤ j ≤ j1 are the indices of all the layers belonging to the 50 m section number k then βi should be comprised between [0.85 − 1.15]hi/ti with $$t_i=t^{corr}_k-\sum _{j=\lbrace j_0:j_1\rbrace ,j\ne i} h^*_j/\beta _j$$. For j0 and j1, only part of the thicknesses that belongs to the section is considered. Because this constraint strongly relies on assignment of the average velocity from interferometry, we allow the constrained velocity to vary by 15 per cent during the inversion. Also, because of this constraint, the inverted velocity model will always be closer to the starting velocity model than to the Deltares-NAM model. The thicknesses are assumed constant during the inversion, and the model is refined iteratively (e.g. Spica et al. 2016, 2017). The layers that are refined between two iterations are chosen based on the sensitivity of the misfit to a small velocity change. We use a Pattern-Search method for the iterating inversion because of the strong nonlinearity of the problem and because of its efficiency (e.g. Audet & Dennis 2002). All the inversions presented in this paper have a misfit value (ε) lower than 0.01. 6 RESULTS AND DISCUSSION We show examples of joint inversion results at borehole sites G03, G18, G34 and G56 in Fig. 5. The left panels depict observed H/V (black lines) along with the best fits after inversion (green lines) and H/V computed from the complex Deltares-NAM velocity model (blue lines). The right panel depicts the best velocity model (green line) along with initial model obtained from borehole interferometry (dashed black line) and the Deltares-NAM model (blue lines). Results at the other borehole sites are shown in the Supporting Information (Fig. S3). Figure 5. View largeDownload slide Inversion results at sites G03, G18, G34 and G56. The left panels of each subplot depict the observed H/Vs (black lines) along with the inverted H/Vs (green lines), the theoretical H/Vs corresponding to the complex Deltares-NAM velocity model (blue lines) and the H/V spectral ratio uncertainty range (pink dashed lines). The upper and lower bounds of the autocorrelations are represented by black dashed lines. They are obtained as the maximum and minimum at each frequency from all the autocorrelations computed with 50 windows numbers, which allows convergence. The related velocity models are shown in the right panels of each subplot. Figure 5. View largeDownload slide Inversion results at sites G03, G18, G34 and G56. The left panels of each subplot depict the observed H/Vs (black lines) along with the inverted H/Vs (green lines), the theoretical H/Vs corresponding to the complex Deltares-NAM velocity model (blue lines) and the H/V spectral ratio uncertainty range (pink dashed lines). The upper and lower bounds of the autocorrelations are represented by black dashed lines. They are obtained as the maximum and minimum at each frequency from all the autocorrelations computed with 50 windows numbers, which allows convergence. The related velocity models are shown in the right panels of each subplot. The inverted velocity models are in good agreement with the results of Kruiver et al. (2017). For the first 50 m section, the inverted velocity models at G03, G18 and G34 show the same trends as the Deltares-NAM velocity model, that is, some high contrast with a higher local velocity at approximately the same depth; however, the velocity ratio between these two models is sometimes close to 2. This difference could be explained by the use of the relationship between β, α and ρ, which could be incorrect at shallow depth where the soil is likely fluid-saturated and where shear velocity might be expected to drop dramatically. Allowing the compressional velocity to become a free parameter, however, would result in an unmanageable increase in the number of unknown parameters for the inversion. Below 100 m depth, the Deltares-NAM and inverted H/V models agree well, and they are concordant with the gradient-like trend. This good agreement suggests that the different volumes sampled by the two studies have similar properties. The use of H/V at depth in the inversion problem appears to reduce the non-uniqueness significantly and, potentially, to successfully reconstruct an accurate velocity model. As mentioned in Section 4.3, some of the discrepancies observed between the different sets of velocity models may be attributable to the different sensitivities of the two approaches. An important difference is the greater number, and strength of velocity reductions with increasing depth. This could be due to the strong constraints imposed by the interferometry picked traveltime during the inversion. On the other hand, we note that the theoretical H/V computed from the Deltares-NAM model fits the observed H/V fairly well, especially at shallow depths. This may suggest that we are converging to an accurate velocity model and hence an accurate site characterization. 7 CONCLUSIONS We present new theoretical and empirical results on the computation and the inversion of H/V spectral ratios based on the diffuse field assumption for combinations of receivers at surface and depth. First, we obtained the mean shear traveltimes by applying ambient field interferometry between adjacent sensors inside the boreholes. By cross-correlating the ambient seismic field we were able to retrieve upward and downward propagating S waves between adjacent borehole seismic stations. We used this mean velocity model as a starting velocity model for the joint inversion of the depth-dependent H/V. We found it especially useful to constrain the inversion for velocities deeper than 100 m. The use of five independent H/V measurements at different depths, together with the average shear wave traveltimes, helps reduce the range of acceptable parameters and helps speed convergence of the solution. In this inversion, the theoretical H/V counterparts (i.e. the forward problem) were computed using the DWN method, which allowed efficient modelling of the wave field for different sources at the surface and at depth. We successfully obtained complex VS velocity profiles of the shallow sub-surface at different borehole sites in the Groningen area. Velocity models are globally in good agreement with previous site characterization for the region (Kruiver et al. 2017). Our approach has the potential to reduce uncertainty in modelling the response of the shallow crust, which is an important component of probabilistic seismic hazard analysis at Groningen. Our results also validate the 3D character of the ambient field sources and motivate measurements of H/V at depth. They also demonstrate the power and reliability of the method, which could be applied elsewhere, where other constraints on shallow structure are lacking. The methodology presented here allows reliable recovery of layered velocity structure at boreholes sites. ACKNOWLEDGEMENTS We would like to thank Pauline Kruiver from Deltares for providing us with their velocity models at borehole sites (Kruiver et al. 2017). We thank the KNMI to make data accessible. We appreciate the contribution of two anonymous reviewers. All the figures have been plotted with matplotlib (Hunter 2007) and some of the data processing steps have been performed using obspy (Beyreuther et al. 2010) and pyrocko (available at: http://pyrocko.org/). Research was funded by Shell Global Solutions International B.V. REFERENCES Aki K., 1957. Space and time spectra of stationary stochastic waves, with spectral reference to microtremors, Bull. Earthq. Res. Inst. , 35, 415– 456. Aki K., Richards P., 2002. Quantitative Seismology , University Science Book. Arai H., Tokimatsu K., 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum, Bull. seism. Soc. Am. , 94( 1), 53– 63. Google Scholar CrossRef Search ADS   Audet C., Dennis J.E. Jr, 2002. Analysis of generalized pattern searches, SIAM J. Opti. , 13( 3), 889– 903. Google Scholar CrossRef Search ADS   Bard P.-Y. & SESAME-team 2004. ‘Guidelines for the implementation of the H/V spectral ratio technique on ambient vibrations measurements, processing and interpretations, SESAME European research project EVG1-CT-2000–00026, deliverable d23.12’. Available at: ftp://ftp.geo.uib.no/pub/seismo/SOFTWARE/SESAME/USER-GUIDELINES/SESAME-HV-User-Guidelines.pdf, last accessed 23 October 2017. Bensen G.D., Ritzwoller M.H., Barmin M.P., Levshin A.L., Lin F., Moschetti M.P., Shapiro N.M., Yang Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169( 3), 1239– 1260. Google Scholar CrossRef Search ADS   Berteussen K.A., 1977. Moho depth determinations based on spectral-ratio analysis of NORSAR long-period P waves, Phys. Earth planet. Inter. , 15( 1), 13– 27. Google Scholar CrossRef Search ADS   Beyreuther M., Barsch R., Krischer L., Megies T., Behr Y., Wassermann J., 2010. ObsPy: A Python toolbox for seismology, Seismol. Res. Lett. , 81( 3), 530– 533. Google Scholar CrossRef Search ADS   Bommer J.J. et al.  , 2017. Framework for a ground-motion model for induced seismic hazard and risk analysis in the Groningen gas field, the Netherlands, Earthq. Spectra , 33, 481– 498. Google Scholar CrossRef Search ADS   Bonnefoy-Claudet S., Cornou C., Bard P.-Y., Cotton F., Moczo P., Kristek J., Fah D., 2006. H/V ratio: a tool for site effects evaluation. results from 1-D noise simulations, Geophys. J. Int. , 167( 2), 827– 837. Google Scholar CrossRef Search ADS   Bonnefoy-Claudet S., Köhler A., Cornou C., Wathelet M., Bard P.-Y., 2008. Effects of Love waves on microtremor H/V ratio, Bull. seism. Soc. Am. , 98( 1), 288– 300. Google Scholar CrossRef Search ADS   Bouchon M., 1981. A simple method to calculate Green’s functions for elastic layered media, Bull. seism. Soc. Am. , 71( 4), 959– 971. Bouchon M., 2003. A review of the discrete wavenumber method, Pure appl. Geophys. , 160( 3-4), 445– 465. Google Scholar CrossRef Search ADS   Brocher T.M., 2005. Empirical relations between elastic wavespeeds and density in the Earth’s crust, Bull. seism. Soc. Am. , 95( 6), 2081– 2092. Google Scholar CrossRef Search ADS   Capon J., 1969. High-resolution frequency-wavenumber spectrum analysis, Proc. IEEE , 57( 8), 1408– 1418. Google Scholar CrossRef Search ADS   De Mulder E.F., Geluk M.C., Ritsema I., Westerhoff W.E., Wong T.E., 2003. De ondergrond van Nederland, Geologie van Nederland, deel . Field E.H., Jacob K.H., 1995. A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent, Bull. seism. Soc. Am. , 85( 4), 1127– 1143. García-Jerez A., Piña-Flores J., Sánchez-Sesma F.J., Luzón F., Perton M., 2016. A computer code for forward calculation and inversion of the H/V spectral ratio under the diffuse field assumption, Comput. Geosci. , 97, 67– 78. Google Scholar CrossRef Search ADS   Hunter J.D., 2007. Matplotlib: A 2D graphics environment, Comput. Sci. Eng. , 9( 3), 90– 95. Google Scholar CrossRef Search ADS   Kruiver P.P. et al.  , 2017. An integrated shear-wave velocity model for the Groningen gas field, the Netherlands, Bull. Earthq. Eng. , 15 3555– 3580. Google Scholar CrossRef Search ADS   Lacoss R., Kelly E., Toksöz M., 1969. Estimation of seismic noise structure using arrays, Geophysics , 34( 1), 21– 38. Google Scholar CrossRef Search ADS   Lawrence J.F., Prieto G.A., 2011. Attenuation tomography of the western United States from ambient seismic noise, J. geophys. Res , 116, B06. Lontsi A.M., Sánchez-Sesma F.J., Molina-Villegas J.C., Ohrnberger M., Krüger F., 2015. Full microtremor H/V (z, f) inversion for shallow subsurface characterization, Geophys. J. Int. , 202( 1), 298– 312. Google Scholar CrossRef Search ADS   Lontsi A.M., Ohrnberger M., Krüger F., Sánchez-Sesma F.J., 2016. Combining surface-wave phase-velocity dispersion curves and full microtremor horizontal-to-vertical spectral ratio for subsurface sedimentary site characterization, Interpretation , 4( 4), SQ41– SQ49. Google Scholar CrossRef Search ADS   Malischewsky P.G., Scherbaum F., 2004. Love’s formula and H/V-ratio (ellipticity) of Rayleigh waves, Wave Motion , 40( 1), 57– 67. Google Scholar CrossRef Search ADS   Miyazawa M., Snieder R., Venkataraman A., 2008. Application of seismic interferometry to extract P-and S-wave propagation and observation of shear-wave splitting from noise data at Cold lake, Alberta, Canada, Geophysics , 73( 4), D35– D40. Google Scholar CrossRef Search ADS   Nakamura Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, Railway Technical Research Institute, Quarterly Reports , Vol. 30, No. 1. Nakamura Y., 2000. Clear identification of fundamental idea of Nakamura’s technique and its applications., in Proceedings of the 12th World Conference on Earthquake Engineering , Auckland, New Zealand. Nakata N., Snieder R., 2012. Estimating near-surface shear wave velocities in Japan by applying seismic interferometry to kik-net data, J. geophys. Res.: Solid Earth , 117( B1). Nakata N., Chang J.P., Lawrence J.F., Boué P., 2015. Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry, J. geophys. Res. , 120( 2), 1159– 1173. Google Scholar CrossRef Search ADS   Noorlandt R.P. et al.  , (In review), Characterisation of ground-motion recording stations in the Groningen gas field, J. Seismol.  Olsen K., 2000. Site amplification in the Los Angeles basin from three-dimensional modeling of ground motion, Bull. seism. Soc. Am. , 90( 6B), S77– S94. Google Scholar CrossRef Search ADS   Perton M., Sánchez-Sesma F.J., 2014. Normalization during the processing of noise correlations, in IUGG , Merida. Perton M., Sánchez-Sesma F.J., 2016. Green’s function calculation from equipartition theorem, J. acoust. Soc. Am. , 140( 2), 1309– 1318. Google Scholar CrossRef Search ADS PubMed  Perton M., Sánchez-Sesma F.J., Rodríguez-Castellanos A., Campillo M., Weaver R.L., 2009. Two perspectives on equipartition in diffuse elastic fields in three dimensions, J. acoust. Soc. Am. , 126( 3), 1125– 1130. Google Scholar CrossRef Search ADS PubMed  Perton M., Contreras-Zazueta M.A., Sánchez-Sesma F.J., 2016. Indirect boundary element method to simulate elastic wave propagation in piecewise irregular and flat regions, Geophys. J. Int. , 205( 3), 1832– 1842. Google Scholar CrossRef Search ADS   Perton M., Spica Z., Caudron C., 2017. Inversion of the horizontal to vertical spectral ratio in presence of strong lateral heterogeneity, Geophys. J. Int. , in press, doi:10.1093/gji/ggx458. Piña-Flores J., 2015. Cálculo e inversión del cociente H/V a partir de ruido ambiental, MSc thesis , Universidad Nacional Autónoma de México, Mexico city. Piña-Flores J., Perton M., García-Jerez A., Carmona E., Luzón F., Molina-Villegas J.C., Sánchez-Sesma F.J., 2016. The inversion of spectral ratio H/V in a layered system using the diffuse field assumption (dfa), Geophys. J. Int. , doi:10.1093/gji/ggw416. Rivet D., Campillo M., Sanchez-Sesma F., Shapiro N.M., Singh S.K., 2015. Identification of surface wave higher modes using a methodology based on seismic noise and coda waves, Geophys. J. Int. , 203( 2), 856– 868. Google Scholar CrossRef Search ADS   Roux P., Moreau L., Lecointre A., Hillers G., Campillo M., Ben-Zion Y., Zigone D., Vernon F., 2016. A methodological approach towards high-resolution surface wave imaging of the san Jacinto fault zone using ambient-noise recordings at a spatially dense array, Geophys. J. Int. , 206( 2), 980– 992. Google Scholar CrossRef Search ADS   Salinas V. et al.  , 2014. Using diffuse field theory to interpret the H/V spectral ratio from earthquake records in Cibeles seismic station, Mexico city, Bull. seism. Soc. Am. , 104( 2), 995– 1001. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J. et al.  , 2011a. A theory for microtremor H/V spectral ratio: application for a layered medium, Geophys. J. Int. , 186( 1), 221– 225. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J., Weaver R.L., Kawase H., Matsushima S., Luzón F., Campillo M., 2011b. Energy partitions among elastic waves for dynamic surface loads in a semi-infinite solid, Bull. seism. Soc. Am. , 101( 4), 1704– 1709. Google Scholar CrossRef Search ADS   Scherbaum F., Hinzen K.-G., Ohrnberger M., 2003. Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations, Geophys. J. Int. , 152( 3), 597– 612. Google Scholar CrossRef Search ADS   Spica Z. et al.  , 2015. Velocity models and site effects at Kawah Ijen volcano and Ijen caldera (Indonesia) determined from ambient noise cross-correlations and directional energy density spectral ratios, J. Volcanol. Geotherm. Res. , 302, 173– 189. Google Scholar CrossRef Search ADS   Spica Z., Perton M., Calò M., Legrand D., Córdoba-Montiel F., Iglesias A., 2016. 3-D shear wave velocity model of Mexico and south US: bridging seismic networks with ambient noise cross-correlations (C1) and correlation of coda of correlations (C3), Geophys. J. Int. , 206( 3), 1795– 1813. Google Scholar CrossRef Search ADS   Spica Z., Perton M., Legrand D., 2017. Anatomy of the Colima volcano magmatic system, Mexico, Earth planet. Sci. Lett. , 459, 1– 13. Google Scholar CrossRef Search ADS   Stafleu J., Dubelaar C.W., 2016. Product specification Subsurface model GeoTOP. Vol. 10133. version 1.3, Tech. Rep, 2016 . Available at: https://www.dinoloket.nl/en/subsurface-models, last accessed 23 October 2017. Stafleu J., Maljers D., Gunnink J., Menkovic A., Busschers F., 2011. 3D modelling of the shallow subsurface of Zeeland, the Netherlands, Neth. J. Geosci. , 90( 04), 293– 310. Tsai V.C., 2009. On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, Geophys. J. Int. , 178 ( 3), 1555– 1564. Tuan T.T., Vinh P.C., Ohrnberger M., Malischewsky P., Aoudia A., 2016. An improved formula of fundamental resonance frequency of a layered half-space model used in H/V ratio technique, Pure appl. Geophys. , 173( 8), 2803– 2812. Google Scholar CrossRef Search ADS   Vos P.C., 2015. Origin of the Dutch Coastal Landscape. Long-Term Landscape Evolution of the Netherlands during the Holocene, Described and Visualized in National, Regional and Local Palaeogeographical Map Series , Utrecht University. Wapenaar K., Fokkema J., 2006. Green’s function representations for seismic interferometry, Geophysics , 71( 4), SI33– SI46. Google Scholar CrossRef Search ADS   Weaver R.L., 2005. Information from seismic noise, Science , 307( 5715), 1568– 1569. Google Scholar CrossRef Search ADS PubMed  Weaver R.L., Lobkis O.I., 2004. Diffuse fields in open systems and the emergence of the Green’s function (l), J. acoust. Soc. Am. , 116( 5), 2731– 2734. Google Scholar CrossRef Search ADS   Yamanaka H., Takemura M., Ishida H., Niwa M., 1994. Characteristics of long-period microtremors and their applicability in exploration of deep sedimentary layers, Bull. seism. Soc. Am. , 84( 6), 1831– 1841. SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1.$$\frac{H_{E-W}}{V}$$ and $$\frac{H_{N-S}}{V}$$ for stations at the surface. Same stations as in Figs 1 and 2. Figure S2. Comparison between the Deltares-NAM velocity model and the initial model form interferometry. Most of the initial models have lower velocities in the first 100 m. Figure S3. Inversion results at all borehole sites highlighted in Fig. 1. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Site characterization at Groningen gas field area through joint surface-borehole H/V analysis

, Volume 212 (1) – Jan 1, 2018
10 pages

/lp/ou_press/site-characterization-at-groningen-gas-field-area-through-joint-NlvahuK0IT
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx426
Publisher site
See Article on Publisher Site

### Abstract

SUMMARY A new interpretation of the horizontal to vertical (H/V) spectral ratio in terms of the Diffuse Field Assumption (DFA) has fuelled a resurgence of interest in that approach. The DFA links H/V measurements to Green’s function retrieval through autocorrelation of the ambient seismic field. This naturally allows for estimation of layered velocity structure. In this contribution, we further explore the potential of H/V analysis. Our study is facilitated by a distributed array of surface and co-located borehole stations deployed at multiple depths, and by detailed prior information on velocity structure that is available due to development of the Groningen gas field. We use the vertical distribution of H/V spectra recorded at discrete depths inside boreholes to obtain shear wave velocity models of the shallow subsurface. We combine both joint H/V inversion and borehole interferometry to reduce the non-uniqueness of the problem and to allow faster convergence towards a reliable velocity model. The good agreement between our results and velocity models from an independent study validates the methodology, demonstrates the power of the method, but more importantly provides further constraints on the shallow velocity structure, which is an essential component of integrated hazard assessment in the area. Downhole methods, Interferometry, joint inversion, seismic noise, site effects 1 INTRODUCTION Due to long-term exploitation of a large onshore gas field and subsequent compaction of the reservoir at depth, the Groningen area in the northern Netherlands is subject to ground subsidence and induced earthquakes. To date, the largest induced earthquake occurred near Huizinge in August 2012 and was recorded with a local magnitude ML of 3.6 (moment magnitude M = 3.4). This event revived concerns about the hazard from induced seismicity and prompted the government and the Nederlandse Aardolie Maatschappij (NAM) to take action. As part of the response, NAM built a dense permanent borehole microseismic network covering an area of about 35 × 45 km2. The network (Fig. 1) is composed of ∼70 accelerographs at the surface that are co-located with a 200 m depth borehole in which geophones are installed at depth intervals of 50 m. The network is fully operational and has recorded continuous waveforms since late 2015. Figure 1. View largeDownload slide (A) KNMI-NAM permanent borehole sites (the G-array) in the Groningen gas field area. The inverted coloured triangles refer to the H/Vs shown in Fig 2. The blue contour depicts the outline of the gas field. White areas depict urban centres. (B) The location of the Groningen gas field (red box) in the northern part of the Netherlands. Figure 1. View largeDownload slide (A) KNMI-NAM permanent borehole sites (the G-array) in the Groningen gas field area. The inverted coloured triangles refer to the H/Vs shown in Fig 2. The blue contour depicts the outline of the gas field. White areas depict urban centres. (B) The location of the Groningen gas field (red box) in the northern part of the Netherlands. It is well known that near-surface lithology can strongly influence damage from earthquake shaking by increasing the amplitude and duration of shaking, and by responding nonlinearly to incident seismic waves (e.g. Olsen 2000). For this reason, site characterization is of great importance to seismic hazard analysis (Bommer et al. 2017). Such characterization requires a knowledge of geomechanical properties of the stratigraphy at a site, which can be often challenging to obtain. Over the last few decades, ambient seismic field seismology has emerged as a valuable tool to characterize shear wave (VS) velocity models over all distance scales. For site characterization, these techniques include the inversion of dispersion curves obtained by ambient seismic field correlation using small-aperture arrays such as the spatial autocorrelation (SPAC; Aki 1957) or frequency–wavenumber methods (f-k; Capon 1969; Lacoss et al. 1969). Depending on the aperture of the array, these techniques can recover the local 1-D velocity structure from surface to a depth of up to several hundred metres. Ambient seismic field cross-correlation tomography with high-density arrays is an emerging method to obtain high-resolution 3-D models of the shallow surface (e.g. Nakata et al. 2015; Roux et al. 2016); however, these kinds of experiments are very expensive, and require significant investment for design, deployment, and operation. The ratio of the horizontal to vertical (H/V) components of the ambient seismic field (sometimes referred to as microtremor) is a widely used method to determine a simple velocity model of the subsurface (i.e. one layer over a half-space) and from that the frequency-dependent site response. The reason this method gives reliable results, however, has long been controversial. This was due to the absence of a clear theoretical basis for the measurements, which leads, inevitably, to lack of clarity in its interpretation. It is not even clear which waves comprise the noise field that generates the H/V peak frequencies (Nakamura 2000; Malischewsky & Scherbaum 2004; Bonnefoy-Claudet et al. 2008). The possibility to retrieve the elastodynamic Green’s function between two stations embedded in an elastic medium from the average time-domain cross-correlation of ambient field records (Weaver 2005) led Sánchez-Sesma et al. (2011a) to advance a new theory for H/V. Based on the diffuse field assumption (DFA), Sánchez-Sesma et al. (2011a) linked the H/V of the autocorrelated signal to the ratio of the imaginary parts of the Green’s functions. The observed H/V can be compared to its theoretical counterpart, and provides a basis for estimating the subsurface structure. It represents an opportunity to obtain the local subsurface velocity model with only one short three-component measurement of the ambient field (e.g. only 20 min in Spica et al. (2015)). This capability is of particular interest in geotechnical engineering, seismic exploration, and engineering seismology. To date it has only been applied to limited data sets and geological settings (Salinas et al. 2014; Spica et al. 2015; Lontsi et al. 2015; Piña-Flores 2015; Rivet et al. 2015; Lontsi et al. 2016; García-Jerez et al. 2016; Piña-Flores et al. 2016; Perton et al. 2017). While these initial results are promising, there is still a need to explore the potential of the method and to validate the results against independent information. In this contribution, we pursue the idea of Lontsi et al. (2015) to use multiple H/V measurements obtained from receivers at depth in a borehole to estimate complex geological structure at a site. The Groningen area, is a well-studied area and it has ∼70 borehole sites, such that it represents an ideal natural laboratory to test the method. We start with an overview of the geological horizons of interest that are sampled by the measurements and some background on constraints on structure from previous studies (Section 2). We then outline the theoretical background and seismic processing necessary to compute H/V (Section 4), and the discrete wave number (DWN; Bouchon 1981) method used to compute the H/V ratios for the forward problem in Section 4.2. In Section 5, we describe the estimation procedure for the velocity models at several borehole sites. Finally, we discuss results and compare them with velocity models obtained independently (Section 6). In this paper we only present a few inversions because our focus is primarily on implementation of the method and the validation of the results. 2 GEOLOGICAL SETTING AND EXISTING VELOCITY MODEL The Groningen area is characterized by flat, low-lying topography with altitude close to mean sea level. The horizons of interest (i.e. those sampled by our measurements) are the Paleogene, Neogene and younger deposits overlying the North Sea Supergroup (De Mulder et al. 2003; Vos 2015). This ∼800 m thick layer of unconsolidated sediments contains a large degree of vertical and lateral heterogeneity due to the influence of the last three ice ages and associated sea level fluctuations. The deposits range from fluvial braid plain sands to shallow marine (intertidal) and terrestrial deposits of soft clays to distinct organic-rich peat formations. The uppermost sedimentary sequence is characterized by a succession of fluvial, glacial, and marine deposits that are crosscut by deep subglacial features (‘tunnel valleys’), which are filled with sands and clays and buried under younger sediments. The North Sea Supergroup consists primarily of alternating marine grey sands, sandstones, and clays. Our study benefits from important, independent information in the form of an integrated VS model from the surface to the base of the North Sea Supergroup (Kruiver et al. 2017). This 3-D VS model is a synthesis of three different VS models obtained through a comprehensive set of geological, geotechnical, and geophysical observations; each of them having a different depth range and spatial resolution. The shallowest part of the model extends from the surface to ∼50 m below sea level, and is based on a high-resolution 3-D geological model GeoTOP (Stafleu et al. 2011; Stafleu & Dubelaar 2016), combined with VS distributions with depth for the sediments (Kruiver et al. 2017). At depths of ∼40 m to ∼120 m below sea level, the model is constrained through inversion of Rayleigh wave observations from spatially extensive reflection seismic surveys. The deepest part of the model is based on the Pre-Stack Depth Migration velocity model derived from sonic logs. It ranges from ∼70 m below sea level to the base of the North Sea Supergroup at (∼800 m). These three distinct models were spliced into a single model that covers the full depth range for site response analyses. In what follows, we refer to this velocity model as the Deltares-NAM model. Other shallow velocity models were provided by Noorlandt et al. (In review). We use 1-D VS profiles of the Deltares-NAM model at the borehole sites to validate our analysis. It is important to keep in mind, however, that the resolution and sensitivity of the methods used to construct the Deltares-NAM model and the H/V velocity profiles differ, as discussed in Section 4.3. 3 DATA We used continuous data from the shallow borehole network installed in the Groningen area (Fig. 1; the G-array). Each borehole site includes one accelerograph at the surface that is collocated with four geophones at 50 m depth intervals (−50, −100, −150, −200 m). The continuous data are freely available on the Royal Dutch Meteorological Institute (KNMI) website. We downsample the continuous data to 50 sps, remove instrument response and convert the time-series to velocity prior to processing. 4 H/V FOR A DIFFUSE WAVEFIELD By definition, the H/V spectral ratio corresponds to the square root of the spectral energy ratio of the horizontal amplitudes (with indices 1, and 2) over the vertical direction (index 3; Arai & Tokimatsu 2004) :   $$\frac{\text{H}}{\text{V}}(\mathbf {x},\omega )=\sqrt{\frac{E_1(\mathbf {x},\omega )+E_2(\mathbf {x},\omega )}{E_3(\mathbf {x},\omega )}}.$$ (1)Perton et al. (2009) showed that at an observer location (surface or depth), these spectral energies (i.e. the directional energy densities, as defined in Perton et al. (2009)) are proportional to the average autocorrelations of the diffuse wavefield components, which in turn are proportional to the imaginary parts of the Green’s function ($$Im[\mathcal {G}]$$) (Sánchez-Sesma et al. 2011b):   $$E_i(\mathbf {x},\omega ) = \rho \omega ^2\langle u_i(\mathbf {x},\omega )u_i^*(\mathbf {x},\omega )\rangle \propto -\omega \mbox{Im} [ \mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega )].$$ (2)where ω is the angular frequency, ui(x, ω) is the displacement field in the i direction at a point x, Im[ ] stands for the imaginary part, $$\mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega )$$ is the displacement Green’s function in the direction i at a point x due to the application of a unit point force in the same direction applied at the same point. The symbol * stands for the complex conjugate operator and the product $$u_i(\mathbf {x},\omega )u_i^*(\mathbf {x},\omega )$$ corresponds in the frequency domain to the autocorrelation of the displacement field in the time domain. The brackets 〈 〉 represent averaging over time. In what follows, the dependence on ω and x is implicit. The ∝ means that the expressions are proportional by a factor that is independent of ω and x. Eq. (2) represents the same case as classic ambient seismic field cross-correlations (e.g. Weaver & Lobkis 2004) but for the special case when the source and receiver are the same. Within this theoretical framework, Sánchez-Sesma et al. (2011a) proposed a theoretical description of H/V ratios and suggested that the H/V spectral ratio recorded at a receiver could also be computed in terms of the imaginary part of the GF:   \begin{eqnarray} \displaystyle\frac{\text{H}}{\text{V}}(\mathbf {x},\omega )&=\sqrt{\displaystyle\frac{\left\langle | u_1(\mathbf {x},\omega )|^2 \right\rangle +\left\langle | u_2(\mathbf {x},\omega )|^2 \right\rangle }{\left\langle | u_3(\mathbf {x}, \omega )|^2 \right\rangle }} = \sqrt{\displaystyle\frac{ \mbox{Im}( \mathcal {G}_{11}+\mathcal {G}_{22}) }{ \mbox{Im}(\mathcal {G}_{33})}}. \end{eqnarray} (3)Eq. (3) links the average energy densities (i.e. the ambient field measurements; see Section 4.1) with the Green’s function (i.e. the theoretical counterpart; see Section 4.2) and treats the H/V spectral ratio as an intrinsic property of the medium. It naturally allows for the inversion (see Section 5) of H/V that includes contributions of the full wavefield—that is, including both surface and body waves. 4.1 Observed H/V In the context of the DFA, eq. (2) is only valid when the seismic wave field is equipartitioned, that is, all the incident waves are (e.g. P, S, or Rayleigh waves) have the same energies (Perton et al. 2016). This assumption is difficult to verify and unlikely to be true in the majority of ambient field data. Therefore, some signal processing must be applied to enhance the equipartitioning of the seismic wavefield, just as for traditional ambient seismic field cross-correlation (e.g. Bensen et al. 2007). As in Spica et al. (2015) and Perton et al. (2017), we apply spectral whitening, which corresponds to source deconvolution. Because several sources can act in different frequency bands, the operation consists of normalizing the signals by the source energies computed in each time window and across several frequency bands as $$\tilde{v}_i(\mathbf {x},\omega )=v_i(\mathbf {x},\omega )/\sqrt{\sum _{i=1}^3|v_i(\mathbf {x},\Delta \omega )|^2}$$; where Δω is a frequency band of 2.5 Hz width centred on ω. We work with the particle velocity vj(x, ω) = iωuj(x, ω) to preserve the link to energy. To remove only the spectral envelope, the bandwidth has to be much larger than the oscillations in the spectra. The H/V spectral ratio is therefore computed in terms of the wavefield autocorrelations as:   \begin{eqnarray} \frac{\text{H}}{\text{V}}(\mathbf {x},\omega )&=\sqrt{\displaystyle\frac{\left\langle | v_1(\mathbf {x},\omega )|^2 \right\rangle +\left\langle | v_2(\mathbf {x},\omega )|^2 \right\rangle }{\left\langle | v_3(\mathbf {x},\omega )|^2 \right\rangle }}. \end{eqnarray} (4)Eq. (4) requires that the averaging is performed separately for each component. In that sense, it is different from the calculation of the usual H/V spectral ratio (Nakamura 1989), which corresponds to the average of the ratios: H/V = 〈Hw/Vw〉. The average is computed on the spectra obtained over one day of continuous data that is windowed into sections of 100 s duration with an overlap of 20 per cent. Each time window is demeaned, detrended and bandpass filtered from 0.1 to 10 Hz. The attenuation is generally a severe limitation for a field to become diffuse and then equipartitioned. Several studies have shown this limitation for cross-correlation with large interstation distance (e.g. Lawrence & Prieto 2011). However, here, the interstation distance is null and the attenuation has a similar effect on waves in all directions. Provided there is sufficient nearby sources of noise, attenuation should not be a limitation. We have also shown theoretically that in presence of attenuation the autocorrelation is proportional to the $$Im[\mathcal {G}]$$ for the same media but without attenuation (Perton & Sánchez-Sesma 2014). A possible limitation of the method is, however, the presence of strong anisotropy in the shallow sediment. This is not taken into account in this contribution. Both the evaluation of anisotropy in the area and its assessment through H/V spectral ratio are possible future research directions. Fig. 2 shows different H/V curves computed for the NW–SE line of boreholes highlighted in Fig. 1. By far, the strongest variations in the H/V curves are observed for the surface measurements (black lines in Fig. 2). We attribute the higher amplitude of the surface H/V ratios to a stronger impedance contrast the air–solid interface and a stronger contribution from surface waves. The H/V curves at depth have much lower amplitude variations, but slight changes in their shape reflect structural changes at depth. These are expressions of shallow structural variability from northwest to the southeast in the study area. Figure 2. View largeDownload slide H/V at different depth levels computed along the NW–SE line of stations shown as inverted coloured triangles in Fig. 1. The dashed grey lines in each panel is the $$\sqrt{2}$$ level (see Section 4.3). The H/V in the green frames (G03, G23, G29 and G56) are shown in Fig. 5 along with their inversion results. Figure 2. View largeDownload slide H/V at different depth levels computed along the NW–SE line of stations shown as inverted coloured triangles in Fig. 1. The dashed grey lines in each panel is the $$\sqrt{2}$$ level (see Section 4.3). The H/V in the green frames (G03, G23, G29 and G56) are shown in Fig. 5 along with their inversion results. In order to test that the observed H/V is not biased due to lack of diffusivity of the ambient field, we analyse the variability between $$\frac{\text{H}_{\text{E-W}}}{\text{V}}$$ and $$\frac{\text{H}_{\text{N-S}}}{\text{V}}$$ for stations at the surface (Supporting Information Fig. S1). The fact that these measurements are similar at almost all the sites, suggests a near isotropic illumination by the noise sources (Perton et al. 2017). Equipartition of the wavefield is also discussed in Section 4.3. 4.2 Theoretical H/V In the second term of eq. (3), the $$\text{Im}[\mathcal {G}]$$ components must be related to a geometry and material properties that explain the data. Perton et al. (2017) explored the case of strong lateral heterogeneities along a crater cliff; however, the case of Groningen does not have strong topographic variation of the surface, which allows us to use a 3-D unbounded layered elastic media and smooth variation of the subsoil in the vicinity of the stations as a reasonable starting assumption. Several methods have been proposed to model the $$\mathcal {G}_{ii}$$ components of H/V in a layered medium under the DFA (Sánchez-Sesma et al. 2011b; García-Jerez et al. 2016; Perton & Sánchez-Sesma 2016; Lontsi et al. 2015). Here, we use the Discrete Wave Number method (DWN; Bouchon 2003, 1981). The DWN method is efficient and suitable for solving the forward problem for H/V for stations located at the surface (Sánchez-Sesma et al. 2011a; Perton et al. 2017). In what follows, we discuss its efficiency when several stations are located in a single geometrical and material configuration at different vertical and/or horizontal positions. As in Sánchez-Sesma et al. (2011a), the wavefield is decomposed as a sum of plane waves according to the horizontal component of the wavenumber vector in the radial direction kr. The Green’s function is the sum of two contributions, one due to the P and SV waves and another due to SH waves. A compact form of the displacement P–SV Green’s function component $$\mathcal {G}^n_{ii}(\mathbf {x})$$ inside the layer n, where receiver and point source are superimposed at x = {x, z} and where the vector source is oriented in the same direction i is:   \begin{eqnarray} \mathcal {G}^n_{ii}(\mathbf {x}) &=& \mathcal {G}^{\text{FS}n}_{ii}(\mathbf {x}) \nonumber\\ && +\, \int _{-\infty }^{+\infty } \big( S^{\nearrow Pn} u_i^{\nearrow Pn} + S^{\nearrow SVn} u_i^{\nearrow SVn} +S^{\searrow Pn} u_i^{\searrow Pn} \nonumber\\ && +\, S^{\searrow SVn} u_i^{\searrow SVn} \big) \text{d}k_r \end{eqnarray} (5)For an extended description of the integrand, the interested reader can consult the textbook of Aki & Richards (2002). $$\mathcal {G}^{\text{FS}n}$$ is the incident field as in a full space (FS) having the same properties as the layer n. The wave amplitudes S are unknowns that are solved for at discrete values of kr by enforcing the continuity of displacement and traction at each interface. Superscripts P and SV are for P and SV waves. We write the system of equations in matrix form [A][S] = [B]; where [S] = {..., S↗Pn, S↘Pn, S↗SVn, S↘SVn, ...} is the unknown vector and [B] is the vector associated with the stress and displacement components of the source (i.e. to $$\mathcal {G}^{\text{FS}n}_{ii}$$). [A] is comprised of the amplitude coefficients of the boundary conditions and has dimensions equal to [4(N − 1))*(4(N − 1)] with N being the number of layers. As several sources are considered at different depths and for different directions, the system of equations should be solved several times (i.e. for each source); however, since we consider a single configuration (i.e. the matrix [A] is the same), it is more efficient to calculate [A]−1 just once. Therefore, the wave amplitudes related to each source are calculated as the result of the multiplication [S] = [A]−1[B], rather through inversion by Gaussian elimination for each vector [B]. In this way all the Green’s function’ components $$\mathcal {G}^n_{ii}$$ can be efficiently evaluated at multiple receiver positions and for the several directions. The same procedure applies for the Green’s function associated with the SH contribution. This allows consideration of all H/V for the same borehole site simultaneously during the joint inversion. 4.3 Forward modelling using a complex velocity model We show in Fig. 3 the Deltares-NAM VS model at borehole site G03. In the first 50 m of the model, there are two strong velocity contrasts (C1 and C2 in Fig. 3). The deeper velocity dependence is mainly expressed by a velocity gradient (dashed black line in Fig. 3). Observed (black lines) and computed H/V (blue lines) at each depth level are shown on the right part of the figure. Figure 3. View largeDownload slide Theoretical H/V computed from the Deltares-NAM model (blue lines) and observed H/V (black lines) at borehole G03. The letters C1 and C2 refer to possible strong velocity contrasts in the upper 50 m depth range, and the dashed line refers to the gradient-like part of the velocity model. Figure 3. View largeDownload slide Theoretical H/V computed from the Deltares-NAM model (blue lines) and observed H/V (black lines) at borehole G03. The letters C1 and C2 refer to possible strong velocity contrasts in the upper 50 m depth range, and the dashed line refers to the gradient-like part of the velocity model. The discrepancies between the observed and synthetic H/V suggest that the inversion of the observed H/V will result in a different velocity model, but with certain common characteristics. While the model of Kruiver et al. (2017) is site-specific and has a resolution that changes with depth (see Section 2), the H/V are sensitive to elastic properties over a larger horizontal area that depends on the wavelength of the waves that contribute to the observations (Perton et al. 2009; Piña-Flores et al. 2016). The joint inversion result will therefore give an averaged (vertical and horizontal) velocity model at a borehole site. The H/V at the surface show one clear frequency peak, while H/V at depth shows small oscillations superimposed on a flat trend. As recently detailed by Piña-Flores et al. (2016), H/V is sensitive to both surface and body waves. The H/V spectra from receivers at the free surface mainly originate from the strong contribution of surface waves and show important localized frequency peaks of high amplitude. These frequency peaks are dependent on the existence of strong impedance contrasts at depth. In the case of a simple velocity model of one layer over a half-space (with shear velocity β1 and thickness h1), the main peak frequency $$f_{\beta _1}$$ is well approximated by (Yamanaka et al. 1994; Piña-Flores et al. 2016):   $$f_{\beta _1} = \frac{\beta _1}{4h_1}.$$ (6)When two layers are considered on top of a half space, as suggested by the Deltares-NAM model (i.e. C1 and C2 in Fig. 3), H/V can present two peaks depending on the thickness of each layer and the impedance contrasts at their boundaries. If the thickness of the deepest layer is greater than the shallower layer and a sufficient impedance contrast between these three materials exists, the two main peaks are well separated (Bonnefoy-Claudet et al. 2006; Field & Jacob 1995; Bard & SESAME-team 2004; Piña-Flores et al. 2016). Therefore, if β2 ≫ β1 and h1 ≫ h2, the shear velocity of the superficial layer can also be evaluated using eq. (6) (here $$f_{\beta _1}$$ is the high frequency peak) and the approximate frequency associated with the peak at low frequencies can be evaluated from   $$f_{\beta _{12}}=\displaystyle\frac{1}{4\left(\displaystyle\frac{h_1}{\beta _1}+\frac{h_2}{\beta _2}\right)},$$ (7)where β2 and h2 are the thickness and shear velocity of the second layer. While more complex formulae that reflect several major effects of the model on the resonance frequency exist (Tuan et al. 2016), modelling of eq. (7) appears to be a good approximation for layered structure (e.g. Piña-Flores et al. 2016); however, since H/V at the free surface in Fig. 3 does not show such character (i.e. a double frequency peak), we would not expect the inverted velocity model to develop two strong velocity contrasts in the first 50 m at this site. Because surface wave energy quickly decreases with depth, H/V in boreholes becomes more sensitive to body waves. While surface waves propagate in 2-D space and are generally not strongly reflected by lateral heterogeneity, body waves propagate in 3-D space and are reflected by the free surface and also by possible strong impedance contrasts at depth. As shown in Perton et al. (2009) for a half space, the waves that travel vertically up and down interfere and result in spectral oscillation periods in the energy density components (E1, E2, E3). The amplitude of these oscillations decays with frequency and depth. In a half-space, the dependence on frequency and depth are similar. The $$H=\sqrt{E_1+E_2}$$ is mainly sensitive to the shear wave velocity while the $$V=\sqrt{E_3}$$ is mainly sensitive to the compressional wave velocity. Also, as a requirement of equipartition, the three energy densities tend to the same value with depth. Therefore, from eq. (3), H/V in a half-space should present these oscillations and tend to $$\sqrt{2}$$ at high frequency and at depth (as observed in Fig. 2). The gradient-like velocity structure suggested by the Deltares-NAM model does not modify these conclusions since the observed H/V at depth expresses these features. The latter suggests that we might be able to evaluate the velocity at each sensor depth; however, the absence of strong impedance contrast at depth reduces our ability to recover a wave reflected away from the buried station and H/V is only weakly sensitive to structure deeper than the deepest sensor (i.e. ≥200 m). Other examples of observed H/V based on the complex velocity model provided by Deltares-NAM are shown with the inversion results below. This allows us to assess how close the site-specific model is from the H/V measurements and to assess the prospects for improvement. 5 INVERSION As discussed in Piña-Flores et al. (2016), consideration of H/V at surface alone is insufficient to characterize shallow properties uniquely since velocities and thicknesses trade-off and lead to a similar H/V. Additionally, the forward problem is highly non-linear and depends on several uncorrelated parameters (Piña-Flores et al. 2016; García-Jerez et al. 2016). Constraining the inversion by adding observations from sensors at depth significantly reduces the possible range of parameters (Lontsi et al. 2015). Starting with an accurate first guess—as, for example, a guess obtained from an independent measurement such as a dispersion curve (Scherbaum et al. 2003; Piña-Flores 2015; Lontsi et al. 2016)—allows more rapid convergence. Although the Deltares-NAM model exists, our goal is to provide independent measurements in our analysis to test the validity of the method. 5.1 Starting models from borehole interferometry We use ambient field borehole interferometry between pairs of adjacent overlying sensors (e.g. Miyazawa et al. 2008) to obtain the mean shear wave traveltime. In Fig. 4(A), we show the geometry of the borehole stations and the sensor pairs used for interferometry. The cross-correlation of the ambient seismic field v(zS) and v(zr) in the frequency domain is given by:   $$C_{\text{AB}}(\omega )= \langle v(z_S)^{\ast } v(z_r) \rangle ,$$ (8)where zS and zr are the depths of the virtual source and of the receiver, respectively, which are any of the consecutive sensors. Then, the imaginary part of the Green’s function between pairs of sensors ($$\text{Im}[\mathcal {G}]_{(z_S,z_r)}$$) is obtained from the average correlation using   $${}\text{Im}[\mathcal {G}]_{(z_S,z_r)} \propto \frac{\langle v(z_S)^{\ast } v(z_r) \rangle }{|S(\omega )|^2},$$ (9)where |S(ω)|2 is the power spectrum of the noise (e.g. Wapenaar & Fokkema 2006). Figure 4. View largeDownload slide (A) Geometry of the boreholes where the grey stars represent the virtual sources (S) and triangles represent boreholes (rb) or surface receivers (rs). Ambient field cross-correlations are performed for every depth level independently. (B) Obtained correlation functions for the north–south and east–west components. The picked traveltime velocities (orange dots) are marked for every depth level and the mean velocity for the entire well is marked in red. Note that similar traveltimes are picked in both horizontal components, suggesting weak horizontal anisotropy at this site. (C) Resulting average velocity model (init) used as starting model for inversion, compared to the Deltares-NAM velocity model harmonically averaged over 50 m intervals. Figure 4. View largeDownload slide (A) Geometry of the boreholes where the grey stars represent the virtual sources (S) and triangles represent boreholes (rb) or surface receivers (rs). Ambient field cross-correlations are performed for every depth level independently. (B) Obtained correlation functions for the north–south and east–west components. The picked traveltime velocities (orange dots) are marked for every depth level and the mean velocity for the entire well is marked in red. Note that similar traveltimes are picked in both horizontal components, suggesting weak horizontal anisotropy at this site. (C) Resulting average velocity model (init) used as starting model for inversion, compared to the Deltares-NAM velocity model harmonically averaged over 50 m intervals. The cross-correlations are computed for both east-west and north-south component pairs between the wave motion observed at a depth zS and another directly overlying sensor zr. Only ambient noise is used and small earthquakes are removed (e.g. Miyazawa et al. 2008) based on the KNMI earthquake catalogue. Correlations are computed over one month of continuous data for 25 s windows. Fig. 4(b) shows the resulting bandpass filtered (1.5–8 Hz) stacked correlations in the time domain. An impulsive arrival is observed on both causal and anti-causal parts of the Green’s function, corresponding to upward and downward propagating S-waves. We estimate the arrival time as in Nakata & Snieder (2012); that is, by seeking the three adjacent samples with the largest amplitude values and applying a quadratic interpolation to find the time at which the wave has maximum amplitude. This time is taken as the traveltime for a shear wave that propagates between the borehole and overlying sensor. We use the average of estimated traveltime on the east-west and north-south components to estimate an average velocity layer. The obtained velocity model at site G43 is shown in Fig. 4(C) and is compared to the Deltares-NAM velocity model harmonically averaged over 50 m intervals. These results might be biased by multiple arrivals caused by shear wave reflection inside the 50 m section due to strong geological contrasts. Besides, if the noise is not fully diffuse, Green’s functions are not well retrieved and some small bias can be observed in the traveltimes (Tsai 2009). In such case, the interferometry picked traveltime would be smaller than the expected one, which results in higher velocities. Because the average velocity models tend to underestimate the Deltares-NAM velocity model at most of the sites (i.e. Fig. 4C and Supporting Information Fig. S2), and because the analysis performed on the H/Vs in Supporting Information Fig. S1, this is not likely to be a strong effect. The discrepancy observed in Supporting Information Fig. S2 is also partly explained by the fact that the interferometry picked traveltimes are directly translated to average velocity. Contrary to the traveltime, the velocity estimated by the interferometry should not be exactly equal to the average of the velocities in the 50 m section since it corresponds to a harmonic average. 5.2 Parametrization, misfit function and inversion The soil properties and the layer thicknesses inside a borehole site are assessed by inverting jointly H/V at several depths, z. The objective function is defined as the root mean square difference between observed and predicted H/V computed with the DWN method:   $$\epsilon = \sqrt{\frac{1}{N_z N_{\omega }}\sum _{z_i}\sum _{\omega } \Bigg ( \frac{H}{V}^{\text{exp}} - \frac{H}{V}^{\text{DWN}}\Bigg )^2; }$$ (10)where Nz = 5 is the number of depths at which zi considered in the joint inversion, and Nω = 14 is the number of points taken in the spectra from 1. to 8. Hz. Fluctuations in the H/V spectra at higher frequency are associated with small, local heterogeneities, which are of minor importance here (Piña-Flores et al. 2016). The only free parameters considered during the inversion are the shear wave velocities (β). To reduce the number of parameters, the compressional velocity (α) and the mass density are assumed to be related to the shear velocity through polynomial relationships (Brocher 2005; Berteussen 1977). We impose that one velocity layer in each section of 50 m (i.e. between each sensor pair) is constrained by the interferometry picked traveltime and by the other velocities of the section. In other words, if i is the layer index of height hi and shear velocity βi, and j0 ≤ j ≤ j1 are the indices of all the layers belonging to the 50 m section number k then βi should be comprised between [0.85 − 1.15]hi/ti with $$t_i=t^{corr}_k-\sum _{j=\lbrace j_0:j_1\rbrace ,j\ne i} h^*_j/\beta _j$$. For j0 and j1, only part of the thicknesses that belongs to the section is considered. Because this constraint strongly relies on assignment of the average velocity from interferometry, we allow the constrained velocity to vary by 15 per cent during the inversion. Also, because of this constraint, the inverted velocity model will always be closer to the starting velocity model than to the Deltares-NAM model. The thicknesses are assumed constant during the inversion, and the model is refined iteratively (e.g. Spica et al. 2016, 2017). The layers that are refined between two iterations are chosen based on the sensitivity of the misfit to a small velocity change. We use a Pattern-Search method for the iterating inversion because of the strong nonlinearity of the problem and because of its efficiency (e.g. Audet & Dennis 2002). All the inversions presented in this paper have a misfit value (ε) lower than 0.01. 6 RESULTS AND DISCUSSION We show examples of joint inversion results at borehole sites G03, G18, G34 and G56 in Fig. 5. The left panels depict observed H/V (black lines) along with the best fits after inversion (green lines) and H/V computed from the complex Deltares-NAM velocity model (blue lines). The right panel depicts the best velocity model (green line) along with initial model obtained from borehole interferometry (dashed black line) and the Deltares-NAM model (blue lines). Results at the other borehole sites are shown in the Supporting Information (Fig. S3). Figure 5. View largeDownload slide Inversion results at sites G03, G18, G34 and G56. The left panels of each subplot depict the observed H/Vs (black lines) along with the inverted H/Vs (green lines), the theoretical H/Vs corresponding to the complex Deltares-NAM velocity model (blue lines) and the H/V spectral ratio uncertainty range (pink dashed lines). The upper and lower bounds of the autocorrelations are represented by black dashed lines. They are obtained as the maximum and minimum at each frequency from all the autocorrelations computed with 50 windows numbers, which allows convergence. The related velocity models are shown in the right panels of each subplot. Figure 5. View largeDownload slide Inversion results at sites G03, G18, G34 and G56. The left panels of each subplot depict the observed H/Vs (black lines) along with the inverted H/Vs (green lines), the theoretical H/Vs corresponding to the complex Deltares-NAM velocity model (blue lines) and the H/V spectral ratio uncertainty range (pink dashed lines). The upper and lower bounds of the autocorrelations are represented by black dashed lines. They are obtained as the maximum and minimum at each frequency from all the autocorrelations computed with 50 windows numbers, which allows convergence. The related velocity models are shown in the right panels of each subplot. The inverted velocity models are in good agreement with the results of Kruiver et al. (2017). For the first 50 m section, the inverted velocity models at G03, G18 and G34 show the same trends as the Deltares-NAM velocity model, that is, some high contrast with a higher local velocity at approximately the same depth; however, the velocity ratio between these two models is sometimes close to 2. This difference could be explained by the use of the relationship between β, α and ρ, which could be incorrect at shallow depth where the soil is likely fluid-saturated and where shear velocity might be expected to drop dramatically. Allowing the compressional velocity to become a free parameter, however, would result in an unmanageable increase in the number of unknown parameters for the inversion. Below 100 m depth, the Deltares-NAM and inverted H/V models agree well, and they are concordant with the gradient-like trend. This good agreement suggests that the different volumes sampled by the two studies have similar properties. The use of H/V at depth in the inversion problem appears to reduce the non-uniqueness significantly and, potentially, to successfully reconstruct an accurate velocity model. As mentioned in Section 4.3, some of the discrepancies observed between the different sets of velocity models may be attributable to the different sensitivities of the two approaches. An important difference is the greater number, and strength of velocity reductions with increasing depth. This could be due to the strong constraints imposed by the interferometry picked traveltime during the inversion. On the other hand, we note that the theoretical H/V computed from the Deltares-NAM model fits the observed H/V fairly well, especially at shallow depths. This may suggest that we are converging to an accurate velocity model and hence an accurate site characterization. 7 CONCLUSIONS We present new theoretical and empirical results on the computation and the inversion of H/V spectral ratios based on the diffuse field assumption for combinations of receivers at surface and depth. First, we obtained the mean shear traveltimes by applying ambient field interferometry between adjacent sensors inside the boreholes. By cross-correlating the ambient seismic field we were able to retrieve upward and downward propagating S waves between adjacent borehole seismic stations. We used this mean velocity model as a starting velocity model for the joint inversion of the depth-dependent H/V. We found it especially useful to constrain the inversion for velocities deeper than 100 m. The use of five independent H/V measurements at different depths, together with the average shear wave traveltimes, helps reduce the range of acceptable parameters and helps speed convergence of the solution. In this inversion, the theoretical H/V counterparts (i.e. the forward problem) were computed using the DWN method, which allowed efficient modelling of the wave field for different sources at the surface and at depth. We successfully obtained complex VS velocity profiles of the shallow sub-surface at different borehole sites in the Groningen area. Velocity models are globally in good agreement with previous site characterization for the region (Kruiver et al. 2017). Our approach has the potential to reduce uncertainty in modelling the response of the shallow crust, which is an important component of probabilistic seismic hazard analysis at Groningen. Our results also validate the 3D character of the ambient field sources and motivate measurements of H/V at depth. They also demonstrate the power and reliability of the method, which could be applied elsewhere, where other constraints on shallow structure are lacking. The methodology presented here allows reliable recovery of layered velocity structure at boreholes sites. ACKNOWLEDGEMENTS We would like to thank Pauline Kruiver from Deltares for providing us with their velocity models at borehole sites (Kruiver et al. 2017). We thank the KNMI to make data accessible. We appreciate the contribution of two anonymous reviewers. All the figures have been plotted with matplotlib (Hunter 2007) and some of the data processing steps have been performed using obspy (Beyreuther et al. 2010) and pyrocko (available at: http://pyrocko.org/). Research was funded by Shell Global Solutions International B.V. REFERENCES Aki K., 1957. Space and time spectra of stationary stochastic waves, with spectral reference to microtremors, Bull. Earthq. Res. Inst. , 35, 415– 456. Aki K., Richards P., 2002. Quantitative Seismology , University Science Book. Arai H., Tokimatsu K., 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum, Bull. seism. Soc. Am. , 94( 1), 53– 63. Google Scholar CrossRef Search ADS   Audet C., Dennis J.E. Jr, 2002. Analysis of generalized pattern searches, SIAM J. Opti. , 13( 3), 889– 903. Google Scholar CrossRef Search ADS   Bard P.-Y. & SESAME-team 2004. ‘Guidelines for the implementation of the H/V spectral ratio technique on ambient vibrations measurements, processing and interpretations, SESAME European research project EVG1-CT-2000–00026, deliverable d23.12’. Available at: ftp://ftp.geo.uib.no/pub/seismo/SOFTWARE/SESAME/USER-GUIDELINES/SESAME-HV-User-Guidelines.pdf, last accessed 23 October 2017. Bensen G.D., Ritzwoller M.H., Barmin M.P., Levshin A.L., Lin F., Moschetti M.P., Shapiro N.M., Yang Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169( 3), 1239– 1260. Google Scholar CrossRef Search ADS   Berteussen K.A., 1977. Moho depth determinations based on spectral-ratio analysis of NORSAR long-period P waves, Phys. Earth planet. Inter. , 15( 1), 13– 27. Google Scholar CrossRef Search ADS   Beyreuther M., Barsch R., Krischer L., Megies T., Behr Y., Wassermann J., 2010. ObsPy: A Python toolbox for seismology, Seismol. Res. Lett. , 81( 3), 530– 533. Google Scholar CrossRef Search ADS   Bommer J.J. et al.  , 2017. Framework for a ground-motion model for induced seismic hazard and risk analysis in the Groningen gas field, the Netherlands, Earthq. Spectra , 33, 481– 498. Google Scholar CrossRef Search ADS   Bonnefoy-Claudet S., Cornou C., Bard P.-Y., Cotton F., Moczo P., Kristek J., Fah D., 2006. H/V ratio: a tool for site effects evaluation. results from 1-D noise simulations, Geophys. J. Int. , 167( 2), 827– 837. Google Scholar CrossRef Search ADS   Bonnefoy-Claudet S., Köhler A., Cornou C., Wathelet M., Bard P.-Y., 2008. Effects of Love waves on microtremor H/V ratio, Bull. seism. Soc. Am. , 98( 1), 288– 300. Google Scholar CrossRef Search ADS   Bouchon M., 1981. A simple method to calculate Green’s functions for elastic layered media, Bull. seism. Soc. Am. , 71( 4), 959– 971. Bouchon M., 2003. A review of the discrete wavenumber method, Pure appl. Geophys. , 160( 3-4), 445– 465. Google Scholar CrossRef Search ADS   Brocher T.M., 2005. Empirical relations between elastic wavespeeds and density in the Earth’s crust, Bull. seism. Soc. Am. , 95( 6), 2081– 2092. Google Scholar CrossRef Search ADS   Capon J., 1969. High-resolution frequency-wavenumber spectrum analysis, Proc. IEEE , 57( 8), 1408– 1418. Google Scholar CrossRef Search ADS   De Mulder E.F., Geluk M.C., Ritsema I., Westerhoff W.E., Wong T.E., 2003. De ondergrond van Nederland, Geologie van Nederland, deel . Field E.H., Jacob K.H., 1995. A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent, Bull. seism. Soc. Am. , 85( 4), 1127– 1143. García-Jerez A., Piña-Flores J., Sánchez-Sesma F.J., Luzón F., Perton M., 2016. A computer code for forward calculation and inversion of the H/V spectral ratio under the diffuse field assumption, Comput. Geosci. , 97, 67– 78. Google Scholar CrossRef Search ADS   Hunter J.D., 2007. Matplotlib: A 2D graphics environment, Comput. Sci. Eng. , 9( 3), 90– 95. Google Scholar CrossRef Search ADS   Kruiver P.P. et al.  , 2017. An integrated shear-wave velocity model for the Groningen gas field, the Netherlands, Bull. Earthq. Eng. , 15 3555– 3580. Google Scholar CrossRef Search ADS   Lacoss R., Kelly E., Toksöz M., 1969. Estimation of seismic noise structure using arrays, Geophysics , 34( 1), 21– 38. Google Scholar CrossRef Search ADS   Lawrence J.F., Prieto G.A., 2011. Attenuation tomography of the western United States from ambient seismic noise, J. geophys. Res , 116, B06. Lontsi A.M., Sánchez-Sesma F.J., Molina-Villegas J.C., Ohrnberger M., Krüger F., 2015. Full microtremor H/V (z, f) inversion for shallow subsurface characterization, Geophys. J. Int. , 202( 1), 298– 312. Google Scholar CrossRef Search ADS   Lontsi A.M., Ohrnberger M., Krüger F., Sánchez-Sesma F.J., 2016. Combining surface-wave phase-velocity dispersion curves and full microtremor horizontal-to-vertical spectral ratio for subsurface sedimentary site characterization, Interpretation , 4( 4), SQ41– SQ49. Google Scholar CrossRef Search ADS   Malischewsky P.G., Scherbaum F., 2004. Love’s formula and H/V-ratio (ellipticity) of Rayleigh waves, Wave Motion , 40( 1), 57– 67. Google Scholar CrossRef Search ADS   Miyazawa M., Snieder R., Venkataraman A., 2008. Application of seismic interferometry to extract P-and S-wave propagation and observation of shear-wave splitting from noise data at Cold lake, Alberta, Canada, Geophysics , 73( 4), D35– D40. Google Scholar CrossRef Search ADS   Nakamura Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, Railway Technical Research Institute, Quarterly Reports , Vol. 30, No. 1. Nakamura Y., 2000. Clear identification of fundamental idea of Nakamura’s technique and its applications., in Proceedings of the 12th World Conference on Earthquake Engineering , Auckland, New Zealand. Nakata N., Snieder R., 2012. Estimating near-surface shear wave velocities in Japan by applying seismic interferometry to kik-net data, J. geophys. Res.: Solid Earth , 117( B1). Nakata N., Chang J.P., Lawrence J.F., Boué P., 2015. Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry, J. geophys. Res. , 120( 2), 1159– 1173. Google Scholar CrossRef Search ADS   Noorlandt R.P. et al.  , (In review), Characterisation of ground-motion recording stations in the Groningen gas field, J. Seismol.  Olsen K., 2000. Site amplification in the Los Angeles basin from three-dimensional modeling of ground motion, Bull. seism. Soc. Am. , 90( 6B), S77– S94. Google Scholar CrossRef Search ADS   Perton M., Sánchez-Sesma F.J., 2014. Normalization during the processing of noise correlations, in IUGG , Merida. Perton M., Sánchez-Sesma F.J., 2016. Green’s function calculation from equipartition theorem, J. acoust. Soc. Am. , 140( 2), 1309– 1318. Google Scholar CrossRef Search ADS PubMed  Perton M., Sánchez-Sesma F.J., Rodríguez-Castellanos A., Campillo M., Weaver R.L., 2009. Two perspectives on equipartition in diffuse elastic fields in three dimensions, J. acoust. Soc. Am. , 126( 3), 1125– 1130. Google Scholar CrossRef Search ADS PubMed  Perton M., Contreras-Zazueta M.A., Sánchez-Sesma F.J., 2016. Indirect boundary element method to simulate elastic wave propagation in piecewise irregular and flat regions, Geophys. J. Int. , 205( 3), 1832– 1842. Google Scholar CrossRef Search ADS   Perton M., Spica Z., Caudron C., 2017. Inversion of the horizontal to vertical spectral ratio in presence of strong lateral heterogeneity, Geophys. J. Int. , in press, doi:10.1093/gji/ggx458. Piña-Flores J., 2015. Cálculo e inversión del cociente H/V a partir de ruido ambiental, MSc thesis , Universidad Nacional Autónoma de México, Mexico city. Piña-Flores J., Perton M., García-Jerez A., Carmona E., Luzón F., Molina-Villegas J.C., Sánchez-Sesma F.J., 2016. The inversion of spectral ratio H/V in a layered system using the diffuse field assumption (dfa), Geophys. J. Int. , doi:10.1093/gji/ggw416. Rivet D., Campillo M., Sanchez-Sesma F., Shapiro N.M., Singh S.K., 2015. Identification of surface wave higher modes using a methodology based on seismic noise and coda waves, Geophys. J. Int. , 203( 2), 856– 868. Google Scholar CrossRef Search ADS   Roux P., Moreau L., Lecointre A., Hillers G., Campillo M., Ben-Zion Y., Zigone D., Vernon F., 2016. A methodological approach towards high-resolution surface wave imaging of the san Jacinto fault zone using ambient-noise recordings at a spatially dense array, Geophys. J. Int. , 206( 2), 980– 992. Google Scholar CrossRef Search ADS   Salinas V. et al.  , 2014. Using diffuse field theory to interpret the H/V spectral ratio from earthquake records in Cibeles seismic station, Mexico city, Bull. seism. Soc. Am. , 104( 2), 995– 1001. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J. et al.  , 2011a. A theory for microtremor H/V spectral ratio: application for a layered medium, Geophys. J. Int. , 186( 1), 221– 225. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J., Weaver R.L., Kawase H., Matsushima S., Luzón F., Campillo M., 2011b. Energy partitions among elastic waves for dynamic surface loads in a semi-infinite solid, Bull. seism. Soc. Am. , 101( 4), 1704– 1709. Google Scholar CrossRef Search ADS   Scherbaum F., Hinzen K.-G., Ohrnberger M., 2003. Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations, Geophys. J. Int. , 152( 3), 597– 612. Google Scholar CrossRef Search ADS   Spica Z. et al.  , 2015. Velocity models and site effects at Kawah Ijen volcano and Ijen caldera (Indonesia) determined from ambient noise cross-correlations and directional energy density spectral ratios, J. Volcanol. Geotherm. Res. , 302, 173– 189. Google Scholar CrossRef Search ADS   Spica Z., Perton M., Calò M., Legrand D., Córdoba-Montiel F., Iglesias A., 2016. 3-D shear wave velocity model of Mexico and south US: bridging seismic networks with ambient noise cross-correlations (C1) and correlation of coda of correlations (C3), Geophys. J. Int. , 206( 3), 1795– 1813. Google Scholar CrossRef Search ADS   Spica Z., Perton M., Legrand D., 2017. Anatomy of the Colima volcano magmatic system, Mexico, Earth planet. Sci. Lett. , 459, 1– 13. Google Scholar CrossRef Search ADS   Stafleu J., Dubelaar C.W., 2016. Product specification Subsurface model GeoTOP. Vol. 10133. version 1.3, Tech. Rep, 2016 . Available at: https://www.dinoloket.nl/en/subsurface-models, last accessed 23 October 2017. Stafleu J., Maljers D., Gunnink J., Menkovic A., Busschers F., 2011. 3D modelling of the shallow subsurface of Zeeland, the Netherlands, Neth. J. Geosci. , 90( 04), 293– 310. Tsai V.C., 2009. On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, Geophys. J. Int. , 178 ( 3), 1555– 1564. Tuan T.T., Vinh P.C., Ohrnberger M., Malischewsky P., Aoudia A., 2016. An improved formula of fundamental resonance frequency of a layered half-space model used in H/V ratio technique, Pure appl. Geophys. , 173( 8), 2803– 2812. Google Scholar CrossRef Search ADS   Vos P.C., 2015. Origin of the Dutch Coastal Landscape. Long-Term Landscape Evolution of the Netherlands during the Holocene, Described and Visualized in National, Regional and Local Palaeogeographical Map Series , Utrecht University. Wapenaar K., Fokkema J., 2006. Green’s function representations for seismic interferometry, Geophysics , 71( 4), SI33– SI46. Google Scholar CrossRef Search ADS   Weaver R.L., 2005. Information from seismic noise, Science , 307( 5715), 1568– 1569. Google Scholar CrossRef Search ADS PubMed  Weaver R.L., Lobkis O.I., 2004. Diffuse fields in open systems and the emergence of the Green’s function (l), J. acoust. Soc. Am. , 116( 5), 2731– 2734. Google Scholar CrossRef Search ADS   Yamanaka H., Takemura M., Ishida H., Niwa M., 1994. Characteristics of long-period microtremors and their applicability in exploration of deep sedimentary layers, Bull. seism. Soc. Am. , 84( 6), 1831– 1841. SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1.$$\frac{H_{E-W}}{V}$$ and $$\frac{H_{N-S}}{V}$$ for stations at the surface. Same stations as in Figs 1 and 2. Figure S2. Comparison between the Deltares-NAM velocity model and the initial model form interferometry. Most of the initial models have lower velocities in the first 100 m. Figure S3. Inversion results at all borehole sites highlighted in Fig. 1. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial