Using horizontal-to-vertical spectral ratios to construct shear-wave velocity profiles

Abstract. For seismic hazard assessment and earthquake hypocentre localization, detailed shear-wave velocity profiles are an important input parameter. Here, we present a method to construct a shear-wave velocity profiles for a deep unconsolidated sedimentary layer by using strong teleseismic phases and the ambient noise field. Gas extraction in the Groningen field, in the northern part of the Netherlands, is causing low-magnitude, induced seismic events. This region forms an excellent case study due to the presence of a permanent borehole network and detailed subsurface knowledge. Instead of conventional horizontal-to-vertical spectral ratios (H∕V ratios) from amplitude spectra, we calculate power spectral densities and use those as input for H∕V calculations. The strong teleseisms provide resonance recordings at low frequencies, where the seismic noise field is too weak to be recorded well with the employed geophones and accelerometers. The H∕V ratios of the ambient noise field are compared with several forward modelling approaches to quality check the teleseism-based shear-wave velocity profiles. Using the well-constrained depth of the sedimentary basin, we invert the H∕V ratios for velocity profiles. A close relationship is observed between the H∕V spectral ratios from the ambient noise field, shear-wave resonance frequencies and Rayleigh-wave ellipticity. By processing only five teleseismic events, we are able to derive shear-wave velocities for the deeper sedimentary sequence with a 7 % bias in comparison with the existing detailed velocity model for the Cenozoic sediments overlying the Groningen gas field. Furthermore, a relation between resonance frequency and unconsolidated sediment thickness is derived, to be used in other areas in the Netherlands, where detailed depth maps are not available.

localisation of the seismic events at the 200 m deep seismometer in each of the 70 stations of the Groningen borehole network . The depth of the base North Sea Group is based on lithostratigraphy of hundreds of production wells in the Groningen gas field (Van Dalfsen et al., 2006). Therefore this depth is a very reliable parameter to base the calculations of the shear wave velocities on. The shear-wave velocity model from Kruiver et al. (2017) acts as reference model to compare our model with.

2 Geological Setting
The area of interest is the northeast of the province of Groningen, with a flat topography, close to mean sea-level and a water table almost up to the surface. The sedimentary cover of interest is formed by the Cenozoic soft sediments, named the North Sea Group (NSG), including the Upper, Middle and Lower North Sea Groups ( Fig. 1) (Mulders, 2003;Vos, 2015). The top of the underlying Upper Cretaceous limestones (Chalk Group) is defined as reference horizon in previous site response studies 75  and recent hazard models Kruiver et al., 2017).
The North Sea Group has an average thicknesses of about 800m, with local variations in depth due to underlying salt diapirs and basin deepening towards the north. Paleogene and Neogene sediments mostly consist of transgressive marine clays alternated with continental siliciclastic sediments. At around 500-600m depth, a high velocity layer (Brussel Sand Member) of compacted sands is deposited during a phase of regression and subsequently partly eroded. Marine delta slope deposits 80 consisting of clays and sands are dominating the upper part of the North Sea Group, unconformibly overlain by Quaternary glacial deposits (Wong et al., 2007). The Pleistocene sedimentary sequence is characterized by fluvial, eolian, glacial and marine deposits, that are crosscut by deep (down to 200m) glacial tunnel valleys which are filled with sand and compacted clays, named the Peelo formation. Very soft Holocene tidal channel sands, marine clays and peat form the upper 10-20m of the Groningen subsurface (Meijles, 2015;Wong et al., 2007). 3 Existing velocity models Kruiver et al. (2017) have built an integrated shear-wave velocity model for the Groningen area from the surface to the base of NSG. In this 3-D velocity model, large amounts of geological, geotechnical and geophysical observations are merged. Due to the rich data availability for the upper 40m, the model is very detailed in this depth range (Noorlandt et al., 2018). At depths of 40-120m, less detail is available as the velocity model model is constrained through inversion of surface waves (recorded as 90 noise) from reflection seismic surveys. shear-wave velocities in the deepest part are derived from a pre-existing pressure-wave velocity model from reflection seismics and the conversion ratio is based on two borehole sonic logs (Zeerijp and Borgsweer).
We refer to this velocity model as the Deltares-NAM model. Mainly for the deeper section of the NSG, the uncertainties are high for the shear-wave velocities since they are derived from a single V p /V s ratio. depth constraints to overcome weak sensitivities of surface waves at these depths.
Although several high quality shear-wave velocity profiles across different scales became available over the last years, in situ velocity measurements for each seismic station in the Groningen network are missing and form the objective of our work.

Data set
To conduct the present study, the KNMI shallow borehole network (further referred to as G-network) was chosen, which has 115 both borehole and surface seismographs. The G-network installed on top of the Groningen gasfield consists of 69 stations  and Fig. 2 shows the distribution of the seismic stations in Groningen. Each station is equipped with threecomponent, 4.5 Hz seismometers at 50m depth intervals (50, 100, 150, 200 m) and an accelerometer at the surface. The stations are continuously recording since 2015 and the data is available via the data portal of Royal Netherlands Meteorological Institute (KNMI, 1993). We use a large three-component data set of ambient noise and five strong teleseimic events (Table 1) measured 120 by the G-network. Several strong teleseismic events are required to obtain an estimate of the H/V ratios and the standard deviation thereof. We refer to 'station' for the entire string with 1 accelerometer and 4 geophones, and refer to 'seismometer' for a single sensor measurement at a certain depth. For calibration purposes we use one broadband seismometer which is nearly co-located with one of the geophones.

125
Multiple strong (M>7.0) global teleseismic phases are recorded within the G-network. With their arrivals, body waves resonate between the free surface and the large impedance contrast between the soft sediment layer and the seismic bedrock, resulting in both vertical and horizontal motions. Teleseismic S-phases are transmitted into the NSG with angles smaller than 5 degrees with vertical. The size of the underlying impedance contrast determines how strongly the waves are reflected and thus how much of the wave energy is trapped. The interference of multiple reverberations within the soft layer leads to a resonance 130 pattern in which certain frequencies are amplified and others interfere destructively.
From five teleseismic arrivals, 1000 s time windows are taken, starting just prior the S-phase onset. The direct S phase, the later coda and the subsequent weaker shear-wave phases are included in this window. Figure 3 shows a seismogram for the Mexico M8.1 event for a soft-sedimentary site in Groningen and a hard-rock site in the south of the Netherlands. Both stations have a comparable distance to the event (9359 vs 9378 km). A time window around the P-phase and S-phase is shown. Their amplitude is larger at (a) than at (b) mainly due to amplification over the soft sediments. Also, the basin setting for (a) leads to stronger coda. The response is band-pass filtered between 0.03 and 2 Hz. This is a low frequency band for the equipped geophones and accelerometers. The lower part of this band, from 0.03 to about 0.2 Hz, is indeed most of the time dominated by instrumental noise. Strong teleseismic arrivals, however, are still well recorded down to 0.03 Hz. Their amplitudes are up 140 till a factor 1000 higher than the microseism. Events with magnitudes above 7.0 and with distances up till 90 degrees, produce S-phases with sufficient energy to compute H/V curves that contain resonances of the entire sedimentary sequence.

Identification of resonance frequencies using teleseismic phases
Conventional H/V curves are computed based on Fourier amplitude spectra (FAS) (e.g. in Bard (2002), or with response spectra (Zhu et al., 2020). Here we present how we calculate H/V curves from power spectrum densities (PSDs). Computation of a 145 PSD starts with computing a FAS, which is the absolute value of the discrete Fourier transform: |F |. Computing the PSD amounts to where ∆t is the sample duration and T the duration of the time window that was used to compute the FAS. The factor of two only needs to be applied when merely positive frequencies are used. The PSD is a frequency-normalized version of the power 150 spectrum, which normalization makes the spectrum less sensitive to input data duration T (Havskov and Alguacil, 2004). Our From a 1000s window, containing the P-and S-phase of the teleseimic arrival, PSDs are computed using time windows of 102.4s with a 75% overlap (Peterson, 1993). The two horizontal components are averaged by vector summation and subsequently the horizontal-to-vertical ratio is taken from the PSDs, resulting in the body wave H/V ratio (HVBW): where the subscripts E , N and Z denote the east, north and vertical component, respectively. When substituting Eq. 1 into the 160 above equation, one finds that the PSD-based HVWB is consistent with a FAS-based H/V computation: in which the FAS of both horizontal components is combined through vector summation. No smoothing is applied since smoothing will affect the picking of the exact fundamental resonance peak. For sites with multiple significant peaks, smoothing could result in a shift to higher frequencies and a bias in identification of the fundamental resonance peak (Zhu et al.,165 2020). Albarello and Lunedei (2013) present several alternative procedures to average the horizontal components for computing ambient noise H/V ratios, and conclude that vector summation is the most physically logic averaging method, therefore this averaging procedure has been applied on the two horizontal-component PSDs. Nevertheless we tested the different methods and found no difference in resulting peak frequencies for our dataset, but biases may arise between the different averaging methods when evaluating the amplitude of the H/V peak (Albarello and Lunedei, 2013). interference of the notch with f 0 leads also to a shift in frequency. Since the frequency perturbation is less than 0.01 Hz, it is decided to calculate the mean HVBW over all depth levels to better cancel out noise and to obtain a more stable result for for a vertically incident shear wave, recorded at the 5 different depth levels. Both the fundamental mode and first higher mode of the entire sedimentary column can be seen. Also, waves resonate between the free surface and the high velocity layer at 450m depth (Brussels Sands).
The resulting fundamental mode interferes with the first higher mode of the entire NSG and results in the elevated peak just below 0.5 Hz.
We compute the HVBWs (Eq. 2) for five teleseismic arrivals (Table 1) and Fig. 5 shows example HVBW curves for station G24. Note that, though the amplitudes are slightly different, the frequency of the main peak (f 0 ) is consistent from event to 180 event. Frequencies above 0.7 Hz are omitted since the S-phase provides insufficient energy above this frequency threshold.
By applying the H/V method on teleseismic phases, the source and path effect are largely suppressed to focus on the site effect. However, this suppression is always incomplete. Furthermore, we assume vertical incidence but in reality this is not obtained. These two facts explain the difference in amplitudes of the resonance peak for each event, despite the uniform subsurface. Some HVBW curves show multiple peaks near the f 0 (Fig. 5a,d), causing a challenge to distinguish the peak 185 corresponding to the resonance of the entire North Sea Group. Comparing HVBW curves from different events helps to pick the peak belonging to the fundamental resonance frequency and by averaging over multiple events, the error in the mean f 0 is minimized.

Calculating shear-wave velocities
From the HVBW we interpret the prominent low-frequency peak as the fundamental shear-wave resonance frequency (f 0 ) of 190 the complete NSG package. From each HVBW peak (f 0 ) and the well established depth of the base of the unconsolidated sedimentary layer (NSG) (Van Dalfsen et al., 2006), shear-wave velocities can be calculated based on (Tsai, 1970), assuming vertical incident body-wave contribution: where f 0 represents the fundamental resonance frequency, V s is the (harmonically averaged) shear-wave velocity from sur-

Calculating an empirical bedrock-depth function
In the previous section we took advantage of the well constrained depth of the unconsolidated sediments to invert for S-wave velocity. In other areas in the Netherlands, similar sedimentary fills exist, but the depth of the NSG is less well constrained.
Some areas might do have constraints from a few deep wells, but many les than in Groningen. Other areas rely on reflection seismic data and again other areas have neither of the two. For these regions outside Groningen, one could either use equation 4 220 together with an average NSG shear-wave velocity as found in the previous section. Alternatively, one could find an empirical relation connecting depth d with resonance frequency f 0 in the form as introduced by Ibs-von Seht and Wohlenberg (1999). Comparing Eqs. 4 and 6 one can see that if the lossless single-resonating wave assumption is taken (Eq. 4) the factor a and b in Eq. 6 correspond to V s /4 and -1, respectively. 225 Fig. 9 shows the f 0 -depth distribution and least-squares fit of Eq. 6, yielding as a mean model for the NSG sediments: and a depth standard deviation of 73 m. Using the mean depth of the dataset (814 m), this corresponds to a standard deviation of 9 %. In comparison, the model that Ibs-von Seht and Wohlenberg (1999) found in the Aachen area (Germany) is shown.
Clearly, using the latter model would yield a too deep sedimentary fill for the Netherlands. Here we show that the sedimentary    (2004), but instead of a PDF of the PSD, we compute a PDF of the ambient-noise H/V ratios. This distribution we further refer to as PDF HVAN. Figure 10 shows examples. In a conventional H/V plot, there would be a mean H/V value for each frequency, or a mean value and a confidence zone, based on a normally distributed error model. Instead, 250 the PDF shows the complete digitized H/V probability distribution for each frequency. From this distribution, the mean, mode, median and different percentiles, can be taken. Generally, H/V calculations are based on noise measurements of largely varying duration: for less than 1 hour e.g. (Fäh et al., 2001;Scherbaum et al., 2003), a few days (Ruigrok et al., 2012) or up to a month (Parolai et al., 2002). Due to the permanent deployment of the network in Groningen, noise has been recorded for almost three years. This creates an excellent opportunity 255 to test what minimum window of noise data yields a stable H/V curve. Furthermore, by measuring a long period of noise, the HVAN can be expressed as a PDF for any desired period for each seismometer in the network. Figure 10 shows distribution 15 https://doi.org/10.5194/se-2020-86 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License.
plots of respectively one week, one month and one year of the calculated PDF HVAN. A stable peak around 0.4 Hz is created by at least one month of recording of the ambient noise field. This one month is the preferred duration because of manageable data handling instead of a year of data. Therefore, for all stations at 50 m depth, one month of ambient noise data is used to 260 compute the PDF HVAN. From this distribution, the mean H/V ratio (which we further refer to as HVAN) is extracted. This HVAN is likely descriptive of the site, rather than a remnant of source or path effects.
For the 4.5 Hz seismometers in the G-network, instrument noise is dominating frequencies below 0.2 Hz, where also the fundamental mode of the NSG resonance is expected. Due to large storms across the year, in which the microseism noise 265 surpasses the instrumental noise, the computed PDF HVAN for one year (Fig. 10c) shows a peak around the fundamental resonance frequency (0.15 Hz) of the NSG. In case of an energetic ocean state, the amplitudes in the fundamental mode are similar to the higher mode peak. However, over the year the probability is low and picking the corresponding peak frequency is unreliable. In August 2018, 4 broadband seismometers were installed in Groningen at 100m depth. With this broadband data we are able to confirm the recordings of the fundamental mode of the NSG. Figure 11a) shows the PDF HVAN for 270 broadband station G82B which does have a clearly resolved peak for the fundamental frequency. Due to the unreliability of the fundamental mode peak in the G-network geophones and accelerometers, in the following we focus on the subsequent notch, between 0.3 and 0.4 Hz (Fig. 10). Triggered by a study by Lontsi et al. (2015), the depth dependency of the H/V curve across a downhole array is evaluated in accelerometer, nearly all ground motion falls outside its dynamic range, below 0.3 Hz. The instrumental noise has the same level on all 3 components. As a consequence, the H/V ratio is close to 1 below 0.3 Hz, whereas the geophones still have a small sensitivity down to 0.2 Hz. The remaining pattern of the H/V curves in Fig. 11 are consistent with Fig. 3b, were the body-wave H/V amplitude is decreasing as function of depth. Here, we use the seismometer at 50 meter depth, where anthropogenic noise is considerably less than at the Earth's surface and hence a more stable HVAN can be computed. The software (Geopsy) used 280 for forward modeling assumes measurements at the surface, therefore the 50 m depth station is preferred over the deeper levels.
7 Surface-and body wave forward modeling Commonly, the ambient seismic field is dominated by surface waves (Aki and Richards, 2002;Bonnefoy-Claudet et al., 2006a, b). Ocean-generated seismic noise (microseism) is measured on the Groningen seismometers from 0.1-1.0 Hz. Here, the low-frequency ambient seismic field is composed of microseisms induced by Atlantic-Ocean gravity waves, while the higher perspective on the composition of the microseismic ambient noise field, we perform both body-wave forward modeling (shear wave resonance) and surface-wave forward modeling (Rayleigh-wave ellipticity).

Body-wave forward modeling
Body-wave forward modeling is performed using OpenHVSR, developed by Bignardi et al. (2016). This software computes theoretical transfer functions of layered soil models based on the fast recursive algorithm from Tsai (1970). For the forward 295 model input files we use the subsurface shear-wave velocity profiles as described in Section 5. To the 5-layer over a halfspace velocity models, we add also a 5-layer seismic quality factor. Q-values of resp. 10, 50, 50, 50 and 100 are determined from the borehole stations. Pressure-wave velocities are derived from existing velocity models (Van Dalfsen et al., 2006). The bodywave forward model (Fig. 12) shows modeled shear-and pressure-wave resonance spectra and is combined with the HVAN curve and the HVBW from station G39. The fundamental frequency for the HVBW has quite a good fit with the modeled 300 body-wave resonance.  (Table 1) .

Ellipticity Forward Modeling
Next to body-wave forward modeling, we apply a forward modeling tool for surface waves (Rayleigh waves) to further investigate the wave field composition. Rayleigh waves exhibit an elliptical particle motion confined to the radial-vertical plane (Aki and Richards, 2002). Many studies have shown that the H/V ratio of ambient noise is related to the amplitude ratio of 305 the radial component of the Rayleigh wave over its vertical component (Arai and Tokimatsu, 2004;Bard, 1999;Bonnefoy-Claudet et al., 2006a;Konno and Ohmachi, 1998;Lachetl and Bard, 1994;Maranò et al., 2012Maranò et al., , 2017 and that the H/V ratio does not depend on the source or the path of the Rayleigh wave, but it depends only on the structure beneath the receiving station (Ferreira and Woodhouse, 2007a, b). Rayleigh-wave ellipticity is frequency dependent in case of a layered and more complex subsurface. In case of a strong impedance contrast of factor 2.5-3.0, the ellipticity curve shows singularities where 310 either the vertical or horizontal component vanishes (Wathelet, 2005). In theory, this results in a singular peak at the fundamental resonance frequency and a singular trough at higher frequency (Bard, 2002;Malischewsky and Scherbaum, 2004;Malischewsky et al., 2006). In reality, no full singularities are obtained since the recordings rarely contain pure Rayleigh waves.
The acquired subsurface shear-wave velocity profiles (Section 5) based on teleseismic arrivals, are used as input for Rayleigh-315 wave forward modelling, using the module Gpell in Geopsy (Wathelet et al., 2008). Velocities from the half-space (Chalk Group) are derived by averaging sonic logs from the wells Zeerijp and Borgsweer , which are assumed to be constant over the area. We perform several tests to investigate the sensitivity for properties of the reference horizon (Malischewsky and Scherbaum (2004)), shear-and pressure wave velocity variations, single -and multiple velocity layers over a half space, and the presence of a high velocity layer representing the Brussel Sand Member. These variations are of minor impact 320 on the resulting ellipticity curves.
Since instrument noise is dominating in the frequencies below 0.2 Hz, the fundamental peak of the ellipticity curves cannot be compared to the HVAN of the seismometers in the G-network. According to Maranò et al. (2012), the ellipticity ratio H/V of the fundamental mode does not only exhibit a peak at the fundamental resonance frequency (f 0 ), but also a trough at higher 325 frequencies corresponding to the vanishing of the horizontal component. In the HVAN curve, a trough is well expressed and can be linked to the trough of the ellipticity fundamental mode while the first and second higher modes are hidden in the wide peak around 0.4 Hz and cannot be distinguished as single peaks in the HVAN (Fig.13a).
The misfit between the troughs of the fundamental-mode HVAN and modeled ellipticity is presented in Fig. 13(b). The ma-330 jority of the stations are within the ≤ 0.03 Hz error range. The mean misfit is about zero, which indicates that the computation of HVBW and the further extraction of velocity profiles (Section 5) has not introduced a bias. From the overall good fit between the HVAN curves and ellipticity troughs of the fundamental mode, we can conclude that the input subsurface shear-wave velocity model is on average consistent with the measured ellipticity troughs of HVAN. Rayleigh-wave ellipticity.

Ambient noise field composition
Body-wave (Fig. 12) and surface-wave ( Fig. 13) forward modelling results are compared with HVAN to determine the ambient wave-field composition. It can be noted that there is a small offset in peak frequency between the shear-wave resonance peak for 345 body waves (close to 0.19 Hz) and the ellipticity peak for surface waves (around 0.2 Hz). With an impedance contrast around factor 2.5-3.0 between the NSG and Chalk Group, the offset between the fundamental resonance frequency and ellipticity notch is consistent with the observations from Wathelet (2005). The HVAN (black line in Fig. 12) shows large similarity with the HVBW observed at G39 (blue line) between 0.4 and 1.0 Hz. This indicates that there is a prominent body-wave presence in the ambient noise, for the investigated frequency range, which corresponds with findings of Spica et al. (2018a). At the same 350 time, there is a surface-wave presence, because Fokker and Ruigrok (2019) retrieved Rayleigh waves in this band, using noise recordings from the G-network.
The trough between 0.2 and 0.4 Hz can be seen both on the HVAN. This trough is not well developed in the modeled shear-wave resonance (gray dashed line in Fig. 12). It shows more similarities with the notch in the ellipticity (blue line in Fig.   13). With RayDec it was confirmed that there is indeed a Rayleigh-wave presence in this frequency band. This observation 355 is supported by Kimman et al. (2012), who used a broad-band array in the Groningen area and found well-pickable phase velocities for Rayleigh waves between 0.2 and 0.4 Hz. Interestingly, a similar ellipticity trough can be seen in the HVAN (Fig.   12). This means that there must also be a (minor) contribution of Rayleigh waves in the teleseismic phases, in order to explain a similar notch. Direct Rayleigh waves will be largely attenuated at 0.3 Hz over teleseismic distances. Thus, the Rayleigh wave presence is likely in the coda of S-phases, due to conversions in or near the sedimentary sequence. The new shear-wave velocity profile for each station location is compared with the Deltares-NAM model (Fig. 14), which is constructed using fixed V p /V s ratios over the Groningen area and thus local lithological variations that affect the V s differently than V p are not accounted for. By using direct measurements of teleseismic phases per station, local subsurface variations are directly translated into average shear-wave velocities, which can explain the differences in velocities between the Deltares-NAM model and our model. The mean shear-wave velocity difference between the two models is -48 m/s. This negative mean 375 difference means that HVBW derived velocities are in general higher than the Deltares-NAM velocities. The majority of the stations have less than 100 m/s difference in shear-wave velocity between the two models. Stations G58 and G59 represent the highest misfit in velocities because of the exceptional high V N SG s that we found. As for most unconsolidated sediments, in Groningen the V s increases with depth, which results in surface waves showing their normal dispersive behavior, meaning that lower frequencies propagate faster. When Eq. 4 is used for finding the depth 380 of these normally dispersive sediments, Bignardi (2017) showed that the depth is underestimated. Also we found, when modelling the resonance of a complex Groningen NSG model (Fig. 4), that f 0 is larger than in case of a single layer with the same harmonic average as the sequence of layers. In the normally dispersive multilayer case, the part of the wave field that already reflects before reaching the free surface, interferes with the free surface reflection, leading to an increase of the combined f 0 .
The shift of f 0 is small and highly dependent on the near-surface model that is assumed. Nevertheless, using Eq. 4, a larger f 0 385 results in underestimation of the depth. Similarly, using Eq. 4 for finding an average V s when the depth is known, a larger f 0 leads to an overestimation of the velocity. This finding is consistent with the -on average-overestimation of the velocity that we obtain with respect to the Deltares-NAM model.
Other assumptions that underlie Eq. 4 are that there are no losses and vertically incident body waves. Though f 0 goes up 390 for a more realistic layering, f 0 goes down again when realistic losses are added to the model. Moreover, f 0 would further be reduced if a correction was made for the non-vertical incidence of the teleseismic phases, which were used for computing the HVBW. Since the exact corrections are unknown, it is not straightforward to make these corrections. On average we measure higher shear-wave velocities than the Deltares-NAM model, and this inconsistency might imply the dominant effect of multilayer interference that is not incorporated in this model.

395
Finally, it has been assumed that a local 1D approximation can be made. The main impedance contrast is the base of the NSG, which is laterally smoothly varying over the region. This supports a local 1D assumption for each individual borehole station. Two stations at flanks of salt domes have been excluded from the analysis because of possible 3D effects.
In the near future, the updated V s model will be tested for hypocentre localisation with the Groningen network. Because of 400 the above considerations, it is expected that there is a bias in the absolute values, but that the lateral variations of V s are quite well resolved.

Conclusions
In this paper, we have shown how average shear-wave velocities (V N SG s ) of the soft sediments of the North Sea Group (NSG) can be retrieved from processing five teleseismic events for almost all stations of the Groningen network in the Netherlands. The shear-wave velocities derived with teleseismic body-wave horizontal-to-vertical spectral ratio (HVBW) approximate the velocities in the Deltares-NAM model and show that this method is suitable for a first estimate of a shear-wave velocity 410 model for an entire sedimentary sequence. Furthermore the new model has the advantage that it is derived from 63 direct site measurements and incorporates local subsurface variations. However, we are aware of limitations in our approach since we did not include a detailed layered subsurface model or attenuation and assumed vertical wave incidence. Therefore, the lateral variation we found is likely close to the actual variability over the boreholes. However, there could remain a bias in the absolute numbers. More generally, since many seismic properties are combined and partly interfering into a single spectrum, the inversion of H/V curves is highly non-unique, no matter which underlying assumptions one puts in. Therefore, its main use is limited to the exploration stage.
Besides, the retrieved H/V curves from a month of ambient noise recordings were combined with Rayleigh ellipticity forward models. These forward models were based on shear-wave velocity profiles constructed from teleseismic phases. We have shown that the retrieved ellipticity curves from these input velocities are in good agreement with the horizontal-to-vertical spectral 420 ratio of ambient noise (HVAN) curves.
In order to interpret stable resonance peaks, we computed the H/V ratios from the power spectrum densities, and subsequently computed probability density functions of the H/V ratios. With the permanent deployment of the Groningen network we were able to select the optimal duration of ambient noise recordings. We found that one month of data is sufficient to find a stable distribution for the resonance frequencies observed below 1.0 Hz.

425
For better understanding of the ambient wave-field composition, we performed Rayleigh-wave ellipticity forward modelling and shear-wave resonance forward modelling. By comparison of the HVAN and HVBW curves from the real data with the synthetic models, we found the presence of both body wave and surface waves in the ambient noise field in the low frequency range.
The method to derive average shear-wave velocities from teleseismic arrivals, from ambient noise, or a combination of the two 430 as demonstrated in this paper, is confirmed as a powerful geophysical tool for exploring shear-wave velocities in a sedimentary layer. We have shown that strong teleseismic arrivals are useful in extending the frequency range in which resonance spectra can be found. Furthermore, we used the accurate mapping of sediment thickness over Groningen, together with the observed resonance frequencies, to find an empirical relation between the two. This relation can be used in other areas in the Netherlands where detailed depth maps are not available.

435
Code and data availability. Supplementary information is available on shear wave velocities for each borehole site and corresponding Matlab code to plot the velocity profile for each site Author contributions. Janneke van Ginkel: constructed the manuscript with input from all co-authors, executed the H/V analysis and forward modeling. Elmer Ruigrok: daily advisor, developed PSD and PDF method, performed synthetic analysis, performed text input. Rien Herber, promotor, advisor, performed text corrections