Ambient vibration analysis on seismic arrays to investigate the properties of the upper crust: an example from Herdern in Switzerland

SUMMARY The difﬁculty and the high cost to assess the subsurface properties led to the development of several geophysical techniques. Generally, the focus of a site study is the reconstruction of the S -wave velocity proﬁle down to few tens to hundreds of metres (e.g. 30–300 m), but not the investigation of deeper structures, such as the transition to the crystalline basement. However, such deeper structures are of interest when seismic hazard products have to relate to a reference rock-velocity proﬁle, for example in regional seismic hazard assessment and microzonation studies. To estimate the S -wave velocity proﬁles down to several kilometres, we study the potential of Rayleigh and Love waves at low (down to 0.1 Hz) and high (up to 20 Hz) frequencies using two seismic arrays of increasing size. The small array, with a maximum inter-station distance of 900 m and a recording time of 3 hr, was aimed at constraining the shallow subsurface down to about 350–400 m, while the big one, with a maximum interstation distance of more than 29 km and 23 hr of recording had the goal to constrain the deeper structure. The arrays were deployed in northern Switzerland (east of the village of Herdern) within the Swiss Molasse basin, a sedimentary basin north of the Alps stretching from Lake Constance to Lake Geneva; its thickness increases from 800 to 900 m in the northeast to more than 5 km in the southwest. The seismic data recorded by the two arrays were analysed using the techniques developed for the analysis of small-aperture arrays. The results were inverted for the S -wave velocity proﬁle in two steps: ﬁrst, the Rayleigh and Love wave phase dispersion curves were inverted together. Secondly, the previous dispersion curves were jointly inverted with the measured Rayleigh wave ellipticity angle. The resulting S -wave velocity proﬁles are similar and show agreement with the available geological and geophysical data, conﬁrming the potential of surface waves to investigate deep structures. Moreover, our analysis proves the feasibility of site characterization techniques to large arrays and the possibility to estimate the P - and S -wave velocity proﬁles down to 5 km, deeper than the contrast between Molasse basin and crystalline rock at around 2.1 km.


I N T RO D U C T I O N
The Swiss Molasse basin is the northern foreland of the Alps located between the Alpine chain to the south and the Jura Mountains to the northwest (Fig. 1a). The basin is northeast/southwest oriented and its sedimentary cover decreases towards the northeast (Sommaruga et al. 2012). The Swiss Molasse basin was largely investigated in the past for oil and gas exploration and radioactive waste disposal purposes (Nagra, http://www.nagra.ch) using active (seismic lines) and invasive (borehole drilling) techniques. One of the drilling sites, located in Herdern, was selected as a study site to test the performance of passive seismic techniques. These analyse the ambient seismic vibrations (natural and anthropogenic) propagating through the subsurface. In comparison with active techniques and tomographic studies, the main advantages of passive seismic techniques include the low installation costs and the short time needed for the field experiments. The main disadvantage is the impossibility to identify lateral variations occurring in the subsurface.
The collected seismic data are generally analysed using both single-station and array processing techniques. Single-station methods identify properties of the subsurface underneath the sensor, such as the site fundamental frequency (H/V - Nogoshi & Igarashi 1971;Nakamura 1989) or the Rayleigh wave ellipticity (RayDec -Hobiger et al. 2009), while the array-processing techniques analyse the seismic sensors of an array together, making it possible to determine the velocity and azimuth of seismic waves crossing the array and to identify Rayleigh and Love waves. Different array techniques were developed in the last decades: Aki (1957) introduced SPAC (SPatial AutoCorrelation) calculating the phase velocities in sedimentary layers using circular shape arrays for vertical-component sensors. Lacoss et al. (1969) and Capon (1969) developed the classical and the high-resolution beam-forming techniques, respectively, where the signals from different sensors are delayed and summed in order to identify the velocity and azimuth of the waves. These techniques are in general used for single components and are mostly used on the vertical component. The horizontal components are more Downloaded from https://academic.oup.com/gji/article-abstract/222/1/526/5823134 by guest on 28 May 2020 difficult to analyse. Fäh et al. (2008) extended the high-resolution beam-forming to the horizontal components allowing the identification of Love wave phase velocities on the transverse component and the Rayleigh waves on the vertical and radial components. Poggi & Fäh (2010) proposed the estimation of Rayleigh wave ellipticity together with the dispersion curve analysis. A different approach was developed by Maranò et al. (2012); the authors estimated the wave parameters and the particle motion of Love and Rayleigh waves using a maximum likelihood approach. The final goal of a site characterization analysis is the estimation of the S-wave velocity profile. To achieve this result, the phase velocity dispersion curves, H/V curves and/or the ellipticity information are inverted. Some authors inverted the Rayleigh wave dispersion curves jointly with the H/V information to better identify the bedrock depth and the layers thickness (Arai & Tokimatsu 2005;Picozzi et al. 2005;Parolai et al. 2005);Boxberger et al. (2011) inverted the Rayleigh and Love wave dispersion curves and the H/V information, improving the shallow layer reconstruction. Hobiger et al. (2013) reconstructed the S-wave velocity profiles using the Rayleigh wave dispersion and ellipticity curves.
The above-mentioned methods are generally used to estimate the S-wave velocities from small and medium-size arrays with array aperture of less than a kilometre Michel et al. 2014). In this work, we extend the application of such techniques to large scale arrays to investigate the deeper structures focusing on the transition between the sedimentary rock cover and the crystalline basement at several thousand of metres of depth. To achieve this goal, two arrays of three-component seismic sensors were deployed in northern Switzerland, east of the village of Herdern (Fig. 1a). The area was chosen for three reasons: (1) geological and geophysical information at depth were collected by SEAG in the 1980s with the Herdern 1 borehole (Roth et al. 2010), (2)  The GeoMol model presents the geology of the Swiss Molasse Basin extending from Lake Constance to Lake Geneva in 3-D (Allenbach et al. 2017). The model, realized in cooperation with five partners, was obtained by the interpretation of topographic and geological maps, geological profiles, deep boreholes (more than 60) and seismic lines (more than 5000 km in total). Most of these data were collected by oil and gas companies for subsoil exploration (Figs b and c).
In the following, we present the deployed seismic arrays, the analysis performed using single-station and array processing techniques, an overview of the picked Rayleigh and Love wave dispersion curves and the ellipticities. To reconstruct the velocity profiles, we interpreted the computed curves and inverted them using the dinver toolbox of the Geopsy package (Wathelet 2008). Finally, the inverted P-and S-wave velocity profiles are compared with the available geological and geophysical information.

M E A S U R E M E N T S
Two arrays of seismic sensors were deployed around the village of Herdern, with array centres close to the location of the borehole Herdern 1 (Fig. 2). The small array was deployed in October 2018 and consisted of 16 short-period sensors (Lennartz 5 s) connected to 12 Centaur digitizers. The sensors were installed on metallic tripods which were either placed in holes dug into the ground for a better coupling at soft locations or placed directly on the ground. The deployed configuration consisted of a central station and five rings of increasing radii (planned to be 12.8, 32, 80, 200 and 500 m) with three sensors each. Each ring was rotated clockwise compared to the inner one, giving the final array configuration a spiral shape. The ring radius increase and the rotation angle give a wide range of analysed interstation distances and yield a good array response (Fig. 2 -bottom left-hand side). The minimum and maximum interstation distances were 12.7 and 900 m, respectively, and the recording time was 3 hr. One seismic sensor in the inner deployed ring showed unusual peaks on the vertical component and was not used for the array analysis.
The large array measurement was performed in May 2018, 13 broad-band sensors (Trillium Compact 120 s) were deployed and connected to the same number of Centaur digitizers. Due to the broadband stabilization time, the first 2 hr of signal after the installation were discarded and a total of 23 hr used for the array analysis. To protect the sensors from external influences, they were covered during the entire recording time using protective bells. The deployed configuration consisted of a central sensor, a few metres apart from the centre of the small array and 12 stations arranged in four rings of increasing sizes; this time, mostly for logistical reasons, the rings were not rotated, resulting in an array of three lines of sensors. This array had maximum and minimum interstation distances of 793 m and 20.85 km, respectively. One of the 13 installed sensors, located in the inner ring, showed connection problems between the sensor and the digitizer, probably due to water infiltration during the night and was not used.
Additional seismic noise recordings were provided by three permanently installed short-period stations: TRULL, STEIN and WEIN. These, consisting of three Lennartz 5 s sensors, are part of the permanent Swiss Seismic Network. The inclusion of these permanent stations improved the array response of the large array.
The deployed array configurations were chosen to cover the widest possible frequency ranges and to achieve a good overlap between the array limits of the small and the big arrays. Due to technical field problems (two sensors not fully working during the two measurements), a perfect overlap was not obtained.
The array responses of the two recording arrays are shown in Fig. 2 (bottom row). In Table 1, the information related to the deployed arrays are reported with the relative k min and k max values. These were calculated using the stations recording over the entire recording time and used for the dispersion curve calculation. For high-resolution array methods, we can expect to resolve seismic waves with wavenumbers between k min 2 and k max . k min and k max were obtained by analysing the array responses of the seismic arrays defined by where N is the number of sensors and r is the sensor positions . k min is identified as the diameter of the central lobe and k max as the largest radius for which the theoretical array response (A) is lower than 0.5. The two arrays have wide array responses (k max /k min ): 27.48 for the small and 21.08 for the big. The recorded seismic data were analysed using two single-station  The array analysis was performed using 15 temporary stations for the small array and 12 temporary stations plus 3 permanent stations for the big array. Two deployed sensors (one for the small and one for the big array), not continuously recording over the entire measurement time (3 and 23 hr, respectively), were only used for the single-station analysis. In the following sections, we present all obtained results.

Horizontal-to-vertical spectral ratio (H/V)
The H/V curves were computed by subdividing the entire recording in time windows of 500 s each, filtering them using a high-pass Butterworth filter, tapering with a 5 per cent Tukey window, computing the absolute spectra of each component and smoothing them using a Konno-Ohmachi window (b = 40 -Konno & Ohmachi 1998). The H/V curve was computed using the formula where N(f), E(f) and Z(f) are the amplitude spectra of the north, east and vertical components, respectively. Using the assumption that Rayleigh and Love waves contribute equally to the horizontal components, the curves were divided by √ 2 (Fäh et al. 2001). The averaged H/V curve was finally computed using the geometric mean method, averaging the curves computed for each window.
The H/V curves computed for the small array are shown in Fig. 3(a). The blue and red stars in both plots represent the peaks selected analysing the H/V curves as described in Bard (2004) and the shape of the amplitude spectra (Castellaro 2016), respectively. Using both methods, two peaks were identified in most of the H/V curves: the first one at around 0.9 Hz and the second one between 6 and 7 Hz. The first one has amplitudes between 3 and 4; it was clearly identified only by the stations located in the first three inner rings and not by the outer two. The second peak, instead, has higher amplitude values of up to 5 and was recognized in all curves.
Using the same window length and a wider frequency range (0.01-20 Hz), the H/V curves were computed for the big array ( Fig. 3b). The shape of the H/V curves computed for the big array shows higher variability than the curves calculated for the small one in particular at low and high frequencies. The same picking methods as the small array were applied, and three peaks were identified: one below 0.11 Hz, one at 0.7-1.15 Hz and one between 1.5 and 15 Hz. The biggest difference between the two methods is the first peak at low frequency; it was identified by the first method on the H/V curve plateau (blue stars) in Fig. 3 and at the intersection between the flat area and the beginning of the steep flank by the second technique (red stars). For both methods, the peak amplitude values vary from 10 to 50 for the blue stars and from 6 to 40 for the red ones. The shapes and the high amplitude values are related to the increase at low frequency of the instrumental noise for the three components. The second peak (0.7-1.15 Hz) was identified by both methods; the first method clearly recognized this peak at five sites with amplitude values ranging from 2 to 11, while the second allowed the peak identification at almost all sites. Two amplitude clusters were identified for this peak: one at amplitudes of around 2 and one between 3 and 4.5. The third peak is caused by the shallow subsurface. It varies over wide frequency and amplitude ranges (2-15 Hz and between 2 and 5).

Rayleigh wave ellipticity (RayDec)
The RayDec technique (Hobiger et al. 2009) is a single-station method to determine the Rayleigh wave ellipticity using the random decrement technique (Asmussen 1997). The method emphasizes the Rayleigh waves with respect to other waves. Starting at each zero-crossing of the vertical signal, time windows of a given length are extracted. The two horizontal components are then combined to a single horizontal component in a way to maximize the coherence with the vertical component, taking the quarter-period time shift between the vertical and horizontal components typical for Rayleigh waves into account. The resulting horizontal and vertical component signals for all time windows are summed and the energy of both components is calculated. The Rayleigh wave ellipticity is then obtained as the square root of the ratio of the horizontal and vertical energies. As a single-station technique, RayDec does not determine the wave propagation direction and can thus not distinguish between retrograde and prograde particle motion. Furthermore, different Rayleigh wave modes cannot be separated.
The ellipticity curves were calculated using a 500 s window between 0.2 and 20 Hz for the small array ( Fig. 3c) and 0.01 and 20 Hz for the big array ( Fig. 3d). Most of the curves calculated for the small array have similar shape and show two peaks: a first one at 0.5-1.5 Hz with amplitude values between 1 and 2 and a second one, wider, at 3.2-12.1 Hz with values in the range 1.5-3. The ellipticity curves computed for the big array show, moving from low to high frequency, an unclear trend with different ellipticity values (0.01-0.07 Hz), a steep flank (0.07-0.17 Hz) and a wide peak (0.7-0.8 Hz), similar to the H/V curves (Fig. 3b). In the frequency range 1.2-4.35 Hz, six stations located in the centre and in the eastern parts of the big array show a linear trend and the ellipticity values increase with frequency similar to the curves of the small array (dark green curves). The other seven curves do not show any recognizable trend confirming the variability of the shallow layers at big scale (black curves).

Array analysis
For the array analysis, three methods were used: (1) threecomponent high-resolution fk, (2) WaveDec and (3) MSPAC. In the following section, the dispersion curves, manually picked, are shown for each method using the small and the big array, respectively.

Three-component high-resolution fk (3CFK)
The three-component high-resolution fk analysis (Fäh et al. 2008) extends the original high-resolution algorithm (Capon 1969) to the seismic recordings of a three-component sensor. For the vertical signal, which is related with the Rayleigh wave motion, the processing is the same as for Capon's method. For the two horizontal components, the analysis is angle-dependent. The two horizontal signals of each time window are projected in different angles with a 5 • step, generating radial and transverse signals. These signals are then analysed using Capon's method and identifying the power output with the signal energy in the respective direction. The azimuth resulting in the highest signal energy is identified for both the radial and transverse components independently. Using the energy information of the radial and vertical signals at the same time, the ellipticity of the identified wave can also be estimated (Poggi & Figure 3. H/V curves computed for the small (a) and the big (b) arrays and ellipticity curves computed for the small (c) and big (d) arrays using RayDec. The red stars in (a) and (b) represent the peaks identified using the shape of the amplitude spectra, while the blue stars show the location of the peak directly picked in the H/V curve. The six dark green curves in plot (d) are located in the centre and in the eastern sector of the big array as shown in Fig. 2 Fäh 2010). The information concerning the Rayleigh waves is contained in the vertical and radial components, while the Love waves are identified on the transverse component.
The computed and picked phase velocity dispersion curves are shown for the small array on the left side of Fig. 4. On the transverse component, which is associated with Love waves (a), a unique dispersion curve was picked in the frequency range 1.3-10.4 Hz. The vertical (b) and the radial (c) components correspond to Rayleigh waves. For the vertical component, two dispersion curves are picked; the one at lower frequencies was picked between 1.37 and 4.2 Hz, while the other one starts at 3.2 Hz and ends close to the upper array limit at 18.5 Hz. For the radial component, a single dispersion curve is identified within the array limits (3.4-12.5 Hz). It shows a similar trend as the right curve of the vertical component.
For the analysis of the big array, the picked phase velocity dispersion curves are reported in the right column of Fig. 4. For Love waves, two branches of dispersion curves were identified on the transverse component (d); the mode at lower velocity covers almost the entire frequency range from 0.14 to 1.69 Hz, while the other one is visible in the frequency range 0.31-0.72 Hz with velocities between 3500 and 4000 m s -1 . For the vertical component, two modes were picked (Fig. 4e). The first one is visible in the entire frequency range defined by the array limits, while the second mode is identified in the frequency range 0.24-0.86 Hz. For the radial component, two modes were identified (Fig. 4f). The first one (0.12-0.17 Hz) is similar to the observation on the vertical component. The second curve (0.29-1.12 Hz) has higher velocities than the observed curve on the vertical component.
For the three phase velocity dispersion curves identified using the small array (2 modes on the vertical and 1 on the radial components) and the four curves picked for the big array (2 modes on the vertical and 2 on the radial components), the ellipticity curves were derived following Poggi & Fäh (2010), as the spectral ratio between the same signal on the radial and vertical components. All ellipticity curves are rather flat without strong peaks or troughs (Fig. 5).

Wavefield decomposition (WaveDec)
The Wavefield decomposition method (Maranò et al. 2012) estimates, for a given time window, the particle motion of the predominant wave, either a Rayleigh or a Love wave, using a maximum likelihood (ML) approach. The contribution of this wave to the wave field is then subtracted from the signals and the parameters of the next wave are estimated, adjusting the previously identified waves subsequently. In this way, the parameters of several waves can be estimated at the same time. For Rayleigh waves, the ellipticity angle is also determined and retrograde or prograde particle motions are distinguished. For retrograde particle motion, the ellipticity angle has negative values, and positive ones for prograde. The ellipticity angle is linked to the classical ellipticity by where ξ is the ellipticity angle and the ellipticity. The ellipticity, the ratio between the horizontal and vertical amplitude, ranges between 0 and ∞, and does not contain information on particle motion. Compared to other array techniques, WaveDec improves the estimation of Love and Rayleigh wave contributions of weaker signals. However, the signals on the different sensors have to be more similar and WaveDec works best on homogeneous underground structures. Different parameters were tested for the window length, the maximum number of waves modelled simultaneously, and the gamma value, which controls the model complexity and tunes the algorithm. For γ = 0, the estimation follows a pure Maximum Likelihood approach, for γ = 1, a pure Bayesian Information Criterion approach. Values between 0 and 1 correspond to mixtures of both approaches.
For the small array, the best result was achieved using γ = 0, a maximum number of modelled waves equal to 5 and a window length of 20 seconds. One mode was picked for the Love waves between 3.5 and 9.1 Hz (Fig. 6a). For the Rayleigh waves (Fig. 6b), two segments were picked: a first one between 2.95 and 4.8 Hz, and a second one between 6.15 and 13.9 Hz. The ellipticity angle of the first segment is negative (Fig. 6c), corresponding to a retrograde Rayleigh wave particle motion, while for the second segment, the ellipticity angle curve is positive and mostly flat between 6.0 and 11.5 Hz. Above 11.5 Hz, a strong decrease of the ellipticity angle is observed (Fig. 6d). The particle motion of the second segment changes from prograde to retrograde. The orange box in Fig. 6(d) shows a possible ellipticity angle curve not picked and with negative values. The possible co-existence of two ellipticity angle segments at the same frequency range will be discussed later.
The same analysis was performed for the big array using γ = 0.1 (mainly maximum-likelihood algorithm) for time windows of 200 seconds each, and a number of modelled waves equal to 3. The picked phase velocity dispersion curves are shown in Figs 6(e) and (f); one mode was picked for the Love waves between 0.15 and 0.66 Hz and two for the Rayleigh waves between 0.12 and 0.89 Hz and between 3.45 and 7 Hz. The ellipticity angle of the first mode is negative over the whole picked frequency range, corresponding to retrograde particle motion (Fig. 6g); for the second mode, the ellipticity angle is always positive and the particle motion is prograde (Fig. 6h).

Modified SPatial AutoCorrelation (MSPAC)
The Modified SPatial AutoCorrelation (Bettig et al. 2001) is an implementation of the SPatial AutoCorrelation (SPAC) method of Aki (1957). The initial formulation expected seismic stations located on a circle around a central station. The stations, at equal distance from the centre, are spaced with regular azimuth. The SPatial AutoCorrelation method calculates the cross-correlation function of each station pair and groups the results in rings of increasing size based on the inter-station distance. For each ring the averaged cross-correlation is computed, using at least three station pairs. The averaged cross-correlations correspond to Bessel functions, which are subsequently inverted for the Rayleigh wave dispersion curves. The MSPAC technique uses the same theory as SPAC but gives more flexibility to the station placement, grouping the station pairs in rings defined by minimum and maximum radii. This technique takes the advantage of the SPAC method in investigating lower frequency ranges than the wavenumber-frequency methods and has the possibility to be more easily applied in urbanized areas. In many cases, the computed dispersion curves do not correspond to a unique mode but to an apparent curve. The seismic data recorded using the two arrays were analysed using the MSPAC technique implemented in Geopsy (Köhler et al. 2007) for the vertical component.
The MSPAC technique was applied for both arrays. Fig. 7 shows the Rayleigh wave dispersion curve obtained for the small array; it was computed using 30 rings with sizes from 12.0 to around 872 m. The picked curve extends from 1.45 to 6.84 Hz.
For the big array, due to the bad shape of the computed Bessel functions, we decided not to pick any phase velocity dispersion curve. The bad quality of the cross-correlation functions, in this case, might be related to the size of the array and the investigation of too heterogeneous portions of the subsurface.

Overview of the results
The results of the picked Rayleigh and Love wave dispersion curves are shown in Fig. 8; no information regarding the uncertainties is shown. Using three different array processing techniques, the Rayleigh and Love waves were computed both for the small and the big arrays. The Love waves were picked using the WaveDec and 3CFK methods. The results show a continuous curve in the frequency range 0.14-10.5 Hz; a second curve was picked using 3CFK. This curve, stretching between 0.3 and 0.7 Hz, is at higher velocities. The Rayleigh waves were computed using WaveDec and 3CFK (both arrays) and MSPAC (small array only). The phase velocity dispersion curves of the different methods are in good agreement. Due to the numerous gaps, it is not possible to pick a continuous curve over the entire analysed frequency range for the Rayleigh waves. The dispersion curves computed using WaveDec show higher scatter than the other results.
A summary reporting the ellipticity curves and the ellipticity angles is shown in Fig. 9; it contains the ellipticity curves computed using RayDec for the central station of each array, the curves picked on the vertical and radial components of both arrays and the ellipticity curves derived from the WaveDec method. The Rayleigh wave ellipticity curves computed using the RayDec method were calculated for the central stations of the small (dashed black line) and the big (solid black line) arrays (Fig. 9a). Both show three peaks in the frequency range 0.1-20 Hz; for the big array, these peaks are wide (0.45-1 Hz, 1.7-8 Hz and 9-20 Hz), while for the small array, with the exception of the peak at low frequency, they are narrower and shifted towards higher frequencies (2.5-10.4 Hz and 12.9-20 Hz). The ellipticity curves computed using the 3CFK for the vertical (red) and the radial (green) components are in agreement with the small-array RayDec curve. The WaveDec curves were derived from the ellipticity angles picked for the big and the small arrays. The curves for the big array show a minimum at 0.22 Hz, which is similar to the ellipticity minimum of the RayDec curve computed for the big array but shifted towards higher frequencies and lower ellipticity values. As the particle motion of this ellipticity curve is retrograde over the whole frequency range, the trough does not correspond to a singularity. Two curves were picked for the small array. The first one with frequencies between 2.15 and 3.35 Hz showed retrograde particle motion and is in good agreement with the Ray-Dec curve and with the curve picked for the vertical component using 3CFK for the small array. The second part, between 6.0 and 13.6 Hz, corresponds to prograde particle motion in agreement with the RayDec curve of the small array and with the two curves for the vertical and radial components of 3CFK, picked using the small array. In the frequency range 11.3-13.35 Hz, an ellipticity angle curve with retrograde particle motion was recognized but not picked. It limits, where the low-frequency limit corresponds to k max of the big array and the high-frequency limit to k min 2 of the small array. Figure 9. Summary of the ellipticity curves (left) and the corresponding ellipticity angles (right) using single-station and array techniques for the small (dashed lines) and the big (solid lines) arrays. As RayDec and 3CFK cannot distinguish between retrograde and prograde particle motion, the ellipticity angles computed from these are shown twice, once for each sense of rotation. The RayDec ellipticities are in black, the picked 3CFK curves for the vertical and radial components are in red and green, respectively, and the picked WaveDec curves are shown in blue. extends over the same frequency range as the right ellipticity angle curve picked for the small array and is shown with an orange box in Fig. 6d. The superposition of two curves within the same frequency range will be discussed in the next section.

I N V E R S I O N
To estimate the P-and S-wave velocity profiles, we inverted the dispersion curves and the ellipticity angle curves using the Conditional Neighbourhood Algorithm implemented in the dinver toolbox of the Geopsy package (Wathelet et al. 2004). To investigate how deep we can actually resolve the subsurface, several tests were performed changing the input dispersion and ellipticity angle curves, increasing the number of layers from 10 to 14, and the maximum investigated depth from 500 to 10 000 m. For each selected depth value and using an equal number of layers, several inversion runs were performed to check the stability of our results. P and S waves, bottom depth and density values were attributed to each layer of the velocity profile. The density values were fixed and increase from 2000 to 3400 g cm -3 . Ranges were defined for the other three parameters. As we expect the velocity profile to have increasing velocity with depth, we forced the velocity values to increase with depth. The Conditional Neighbourhood Algorithm generates a random initial set of models; for each model, the misfit between the input data and the generated curves is calculated using where N is the total number of data points, D and M are measured and modelled data points, respectively, and σ i is the uncertainty of the measurement. i indicates the data points used for the inversion. A new generation of models is generated in the neighbourhood of the best models (Wathelet et al. 2004). The selection of the best final model takes into account the computed misfit value.

H/V analysis
The H/V curves computed for the small array show similar shapes up to 10 Hz, indicating a rather homogeneous subsurface below the deployed array (Fig. 3a). The higher variability occurring above 10 Hz is linked to variability in the shallow subsurface (down to about 10 m). The curves computed for the big array, instead, show high variability both at low (below 0.2 Hz) and high (above 1.5 Hz) frequencies (Fig. 3b). At high frequency, the variability is explained taking into account the first hundreds of metres, while for the low frequencies several explanations are possible. We do not think the variability is related to the variation of the deep structures but rather to the short recording time and to the selection of the installation sites. For the purpose of this work and taking into account the size of the big deployed array, we do not use the frequency content below 0.1 Hz. Depending on the investigated site, the H/V curves at low frequency present either an almost flat plateau and a steep trend or a wide peak. In both cases, the high amplitude values (higher than 10) are not reasonable and related to the increase of instrumental noise at low frequency. More stable results and lower amplitude values at low frequency might be obtained by burying the broadband sensors and recording the seismic noise from few days to several weeks.

Mode attribution
In order to run the inversion and find the velocity profile for the investigated area, the mode attribution of the picked dispersion curves plays a crucial role. We start from the assumption that, in any case, both for Love and Rayleigh waves, the fundamental mode is always visible and has lower velocities than for the higher modes. According to Fig. 8(a), the definition of the fundamental mode for Love waves (L0) is straightforward; it consists of the two phase velocity dispersion curves picked using 3CFK method for the big (0.14-1.75 Hz) and small (2.03-10.12 Hz) arrays. The other two dispersion curves picked using WaveDec were not chosen because of the higher scattering.
The Rayleigh wave fundamental mode (R0) was defined taking into account the dispersion curves at lower velocities for the small and big arrays. They do not match as perfect as the Love waves, but show similar slopes. The fundamental mode is given by the left curve picked for the big array (0.1-1.3 Hz) and the two curves picked for the small array. While the left curve picked for the small array is entirely attributed to the fundamental mode (1.9-3.9 Hz), the right one was interpreted as a combination of fundamental mode (8.5-17.45 Hz) and first higher mode (3.2-8.5 Hz). According to the WaveDec results, the ellipticity angle of the left curve picked using the big array has a retrograde particle motion (Fig. 6g) in agreement with the given interpretation (Maranò et al. 2012(Maranò et al. , 2017a; the right curve, instead, has prograde particle motion and it is interpreted as the first higher mode (Fig. 6h). The curves picked using the small array show retrograde motion for the left mode (in agreement with the given interpretation) and prograde for the right one. The left mode is interpreted as fundamental mode, while the right one corresponds to the first higher mode. Over the same frequency range as the right curve, we identified a second ellipticity angle with negative values. It probably corresponds to the extension of the fundamental mode towards higher frequency. The co-existence of two curves confirms the overlap/proximity of two modes and the possibility to have a mix of two modes.
The mode interpretation was obtained using the phase velocity dispersion curves picked using several array processing techniques and the WaveDec ellipticity angle curves. These were useful to distinguish the two overlapping modes for the small array. Even if the curves picked using WaveDec provide the information about the ellipticity angle in agreement with the final mode attribution, the inverted Rayleigh and Love wave dispersion curves consist of the dispersion curves picked using 3CFK. The other curves (WaveDec and MSPAC) were not used due to their scatter and discontinuity over the analysed frequency range. However, the application of the methods allowed confirmation of our interpolation.

Maximum depth
The largest wavelength in our dispersion curves for Rayleigh waves is about 37.5 km and corresponds to a wave velocity of about 3000 m s -1 at a frequency of 0.08 Hz (Fig. 8b). We can assume to resolve a depth of about 1/3 to 1/4 of this wavelength, that is about 10 km, using the phase velocity dispersion curve information alone. Including the ellipticity and the higher mode information, deeper layers can be better characterized if they provide information at lower frequency (Arai & Tokimatsu 2005).
The best result was achieved inverting at the same time the fundamental modes of Rayleigh and Love waves and a short portion of the first higher mode for the Rayleigh waves using a maximum investigated depth of 5000 m and a parametrization consisting of 14 layers over a half-space. Fig. 10 shows the results of the inversion (top row) and the comparison of the synthetic and picked dispersion and ellipticity curves (bottom row).
The theoretical dispersion curves, in orange, were computed using the best velocity model. They show agreement with the picked dispersion curves over the entire frequency range, allow the higher modes to be identified and provide information about the investigated depth. For the Love waves ( Fig. 10 -bottom left-hand side), the mode at higher velocities might correspond to the first higher mode. For Rayleigh waves (Fig. 10 -bottom centre), the first and second higher modes of the best model of the inversion fit well with the measured dispersion curves.
The inversion results of Fig. 10 show the overlap of the picked fundamental modes of Rayleigh and Love waves (R0 and L0) with the theoretical dispersion curves over the entire frequency range. The Rayleigh wave higher mode (R1), in relation with the halfspace S-wave velocity of 3680 m s -1 , is only partially explained in terms of shape (small differences between 3 and 4 Hz) but still inside the standard deviation ranges of the picked curve.
A second inversion was performed combining the phase velocity dispersion curves used in the previous inversion with the ellipticity angle information derived for the central station of the small array (black circle in Fig. 2a). Different tests were performed inverting the ellipticity curves derived using the 3CFK, the RayDec ellipticity curve and WaveDec ellipticity angle or selected portions of them, both at low and high frequency. The best result, shown in Fig. 11, was obtained using the RayDec ellipticity angle curve between 0.22 and 15.8 Hz. Based on the retrograde particle motion Downloaded from https://academic.oup.com/gji/article-abstract/222/1/526/5823134 by guest on 28 May 2020 of the WaveDec ellipticity angle curves for the fundamental modes picked for the small and big arrays, we selected the RayDec derived ellipticity angle curve, assigning negative values and therefore retrograde particle motion. The parametrization consists of 14 layers over a half-space and a maximum investigated depth of 5000 m.
A perfect fit was simultaneously achieved for the Rayleigh wave ellipticity angle curve (R0) and the Rayleigh wave first higher mode (R1) over the entire frequency range. Some small disagreements can be seen for the fundamental mode of Rayleigh waves (R0) between 0.8 and 1.5 Hz, in the frequency range of the H/V peak and the right descending flank. This problem might be related to the fact that the vertical component of the Rayleigh wave fundamental mode is small at the H/V spectral peak and therefore the dispersion curve difficult to measure (Scherbaum et al. 2003). Compared to the results of Fig. 10, also the fit of the fundamental mode of Love waves (L0) decreases. The overlap is less perfect than in the first inversion, but the shapes and the phase velocities are still comparable. The joint inversion of the ellipticity and dispersion curves provides additional constraints to the S-wave velocity profile. The ellipticity curve puts additional constraints mainly in the thickness of the layers, while the dispersion curves constrain the shear-wave velocities (Scherbaum et al. 2003).
The best P-and S-wave velocity models for the two inversions are plotted in Fig. 12 in comparison with the available geological Figure 11. Results of the inversion using the fundamental and first higher Rayleigh wave dispersion curves, the Love wave dispersion curve and the Rayleigh wave ellipticity angle as inversion targets: P-and S-wave velocity profiles and ellipticity angle (first row), dispersion curves of the fundamental mode and first higher mode for Rayleigh waves and fundamental mode dispersion curve of Love waves (second row) and comparison of the picked dispersion and ellipticity curves with the generated curves (orange curves) obtained using the best velocity model resulting from the inversion (third and fourth rows, respectively). The grey curves in the first and second rows correspond to the model with the lowest misfit. The labels show the mode attribution for the inverted curves (first and second rows -bottom left-hand corner) and the picked ones (third row). The colours and the line styles used in the bottom row plots are the same as the ones used in Figs 8 and 9. and geophysical information. The thick dark red and orange velocity profiles represent the S-wave velocity profiles with the lowest misfit without and with the ellipticity angle used in the inversions, respectively; the black and green velocity profiles represent the two corresponding P-wave velocity profiles. The thin dark red and orange lines are the other nine best S-wave velocity profiles from the inversion, the thin black curve is the P-wave velocity profile measured in Herdern 1 (Roth et al. 2010) and the thick dashed profile is the interpretation of the velocity profile measured in Herdern 1 (Sommaruga et al. 2012). To the right, the interpreted geology of the well Herdern 1 is shown down to the top of the crystalline basement (Allenbach et al. 2017).
The borehole Herdern 1 only yielded P-wave velocity information, but in our inversion, the S-wave velocity profiles are the main target, because the P-wave velocity profiles are in general not well resolved in dispersion curve inversions (Wathelet et al. 2004). In our case, due to the availability of only a P-wave velocity profile for Herdern 1 borehole, we compared it with the P-wave velocity profile result from in the inversion.
Using the picked phase velocity dispersion curves and the computed ellipticity information, we reconstructed the P-and S-wave 1-D velocity profiles of the upper crust using a low number of sensors installed on the topographic surface from few hours to a day. From the surface down to 1350 m of depth, the inverted P-wave velocity profile shows agreement with the P-wave velocity profile measured in Herdern; at higher depths, where the Malm deposits are (1348-1600 m), the velocity profile measured in the borehole presents a strong velocity increase followed by a low-velocity layer (1600-1995 m). In the Swiss Molasse basin model (GeoMol) as in the Herdern 1 borehole report, the contrast with the crystalline basement is located at 2129 m. In all inversions, the inverted S-wave velocity profiles do only show a slight velocity increase at the top of the crystalline basement.
The inversions consisted of a blind search, where the picked dispersion curves and the ellipticity angle curves were inverted without previous depth and velocity constraints. Based on the good agreement between the best model and the upper portion of the Herdern velocity profile, a new inversion was performed constraining the thickness and the velocities of the layers below 1300 m to the available geophysical information (Fig. 13). For these layers, as reported in the measured velocity profile, low-velocity zones were allowed. The minimum misfit values obtained from this inversion (0.60) was comparable with the values computed during the previous runs (0.83) and the fit between the best model and the input curves did not improve, meaning that a velocity model with a velocity inversion at 1600 m of depth has phase velocity dispersion curves similar to a model with a velocity gradient.

Performances of the array techniques
Three different array techniques were used to calculate the Rayleigh and Love wave dispersion curves for the small and big arrays. For the analysed work, taking into account the picked dispersion curves (Fig. 8) and the results of the inversion (Fig. 11), 3CFK is the method providing the most continuous and clearest curves for both arrays. This method is less sensitive to the subsurface variation, to the source distribution and to the selected array geometry; WaveDec and MSPAC are more sensitive to the presence of lateral variations when wide areas are investigated. WaveDec allowed the retrieval of dispersion curves for both arrays; compared to 3CFK, they show a much higher scatter. MSPAC, instead, allowed the Rayleigh wave dispersion curve computation only for the small array. The impossibility to compute the dispersion curves also for the big array is probably related to the not perfect location of the sensors in rings of increasing size, to the lack of noise correlation between different pairs of installed stations or to the relatively short recording time for such a big array. Based on our results, we suggest to compute the phase velocity dispersion curves using several array techniques, check the results stability at low and high frequencies and invert Figure 13. Results of the inversion using the fundamental and first higher Rayleigh wave dispersion curves, the Love wave dispersion curve and the Rayleigh wave ellipticity angle as inversion targets including a low velocity zone at depths between 1300 and 2000 m: P-and S-wave velocity profiles and ellipticity angle (first row), phase velocity dispersion curves of the fundamental mode and first higher mode for Rayleigh waves and fundamental mode dispersion curve of Love waves (second row) and comparison of the picked dispersion and ellipticity curves with the generated curves (orange curves) obtained using the best velocity model resulting from the inversion (third and fourth rows, respectively). The grey curves in the first and second rows correspond to the model with the lowest misfit. The labels show the mode attribution for the inverted curves (first and second rows -bottom left-hand corner) and the picked ones (third row). The colours and the line styles used in the bottom row plots are the same as the ones used in Figs 8 and 9. the more continuous and smoothest dispersion curves avoiding the curves with scatter.
As shown in Fig. 9, we calculated the ellipticity curves and the ellipticity angle curves using one single-station technique (RayDec) and two array processing methods (3CFK and WaveDec). Based on the good overlap of the RayDec ellipticity curve computed for the small array with the 3CFK ellipticity curves (vertical and radial components) and its similarity with the results of WaveDec, we performed a joint inversion of this curve, as ellipticity angle, together with the Rayleigh and Love wave dispersion curves. We used eq. (3) to transform the ellipticity curve in ellipticity angle; the curve was calculated both for the positive and negative values and the sense of particle motion assumed. To define which of these two curves to use in the inversion, we compared the resulting curves with the WaveDec ellipticity angles. The curves attributed to the fundamental modes for the small and big arrays show shapes similar to the negative RayDec ellipticity angles and a retrograde particle motion. Taking into account the information provided by WaveDec regarding the particle motion and assuming to have only information concerning the Rayleigh wave ellipticity of the fundamental mode in the selected curve, we inverted the negative RayDec ellipticity angle curve computed for the small array.

C O N C L U S I O N
In the last decades, ambient seismic noise was used to investigate the subsurface at different scales. Depending on the project goal, on the processing techniques and on the available instrumentation, the maximum investigated depth ranges from few tens of metres to tens of kilometres. In this work, we aimed at reconstructing the P-and S-wave velocity profiles of the upper crust down to the contrast between the Swiss Molasse basin and the crystalline basement, using broadband and short period sensors. The thickness of the sedimentary rocks and the depth of the contrast between the Molasse basin and the crystalline basement increases from northeast to southwest (from less than 1 km to more than 5 km). In the study area, the Molasse deposits and the basement contrast was drilled at 2129 m in the Herdern 1 borehole. In order to investigate the subsurface down to more than 2 km, a site characterization study was performed around the location of the Herdern 1 borehole with the deployment of two seismic arrays. The recorded data were analysed using singlestation and array processing techniques. The dispersion curves and the ellipticity angles resulting from the analysis were inverted to reconstruct the 1-D velocity profiles. We ran the inversion assuming a homogenous subsurface below the deployed array without lateral variations (1-D assumption). The inversion was performed using the Rayleigh and Love wave dispersion curves alone and together with the ellipticity angle information.
Based on our results and on the good agreement with the available geological and geophysical data, we confirm (1) the possibility to apply array techniques developed for small arrays to large-aperture ones, (2) the potential to successfully measure Rayleigh and Love wave dispersion curves at low frequencies and (3) the possibility to estimate the velocity profiles down to several kilometres using these data. The combination of the data recorded using the small and large arrays, together with an optimized array configuration (Kennett et al. 2015) and relatively short time recordings, allowed the estimation of the P-and S-wave velocity profiles from the surface down to 5 km. The maximum investigated depth is more than double the transition from the sedimentary rocks to the crystalline basement and half of the initial estimation of 10 km. Higher investigated depths might be achieved by inverting the Rayleigh and Love waves fundamental and higher mode dispersion curves together with the ellipticity angle of the first higher mode. In the inverted S-wave velocity profiles, the transition to the crystalline basement cannot be easily identified because no strong velocity increase is found. The velocities are close to the resulting velocity profiles measured in the Herdern 1 borehole, which was integrated and interpreted in the GeoMol project. In the same way, the measured P-wave velocity profile of Fig. 12 (thin black line) does not show any strong contrast at the transition with the basement but rather a big velocity increase at the Keuper-Muschelkalk interface.
The use of a priori information at large depths (e.g. strong velocity contrasts or velocity inversions) did not improve the quality of the final results, meaning that the measured phase velocity dispersion curves are not able to identify localized velocity inversions. The inverted velocity profile with a strong velocity increase and a velocity inversion gives similar misfit results as a velocity profile with a gradient.
The characterization of the upper crustal structure is of interest when seismic hazard products have to relate to a reference rockvelocity profile, for example in regional seismic hazard assessment and microzonation studies, where the depth range related to the wavelengths of interest in seismic hazard assessment has to be covered. This can be achieved using the same single-station and array processing techniques used in site characterization measurements. Due to the fast array installation and the short recording time, the site characterization analyses have advantages compared to a tomographic approach. In cases of much deeper structure investigations (e.g. the Moho discontinuity), the site selection has to be properly evaluated, the seismic sensors have to be buried (e.g. isolated from the atmospheric variations, covered from the rain, etc.) and the total recording time should be longer.