Ambient Vibration Analysis on Large Scale Arrays When Lateral Variations Occur in the Subsurface: A Study Case in Switzerland

The ambient vibration analysis is a non-invasive and low-cost technique used in site characterization studies to reconstruct the subsurface velocity structure. Depending on the goal of the research, the investigated depth ranges from tens to hundreds of meters. In this work, we aimed at investigating the deeper contrasts within the crust and in particular down to the sedimentary-rock basement transition located at thousands of meters of depth. To achieve this goal, three seismic arrays with minimum and maximum interstation distances of 7.9 m and 26.8 km were deployed around the village of Schafisheim. Schafisheim is located in the Swiss Molasse Basin, a sedimentary basin stretching from Lake Constance to Lake Geneva with a thickness ranging from 800 to 900 m in the north to 5 km in the south. To compute the multimodal dispersion curves for Rayleigh and Love waves and the Rayleigh wave ellipticity angles, the data were processed using two single-station and three array processing techniques. A preliminary analysis of the inversion results pointed out a good agreement with the fundamental modes of Rayleigh and Love waves used in the inversion and a quite strong disagreement with the higher modes. The impossibility to explain at the same time most of the dispersion curves was interpreted as the co-existence, within the investigated area, of portions of the subsurface with different geophysical properties. The hypothesis was confirmed by the Horizontal-to-Vertical spectral analysis (H/V) which indicated the presence of two distinguished areas. The observation allowed a new interpretation and the identification of the Rayleigh and Love wave fundamental modes and of the S-wave velocity profiles to be reconstructed for each investigated zone. It results in two S-wave velocity profiles with similar velocities down to 15 km deferring only in their shallow portions due to the occurrence of a low velocity zone at a depth of 50–150 m at the centre of the investigated area.


Introduction
The ambient seismic noise consists of a mixture of body (P-and S-waves) and surface (Love and Rayleigh waves) waves propagating through the subsurface in different directions and with different phase velocities. Surface waves are widely used in site effect evaluation (e.g. Wathelet et al. 2008;Picozzi. et al. 2009;Bergamo et al. 2011;Hobiger et al. 2016) and tomographic studies (e.g. Stehly et al. 2009;Molinari et al. 2015) with the aim of investigating the subsurface and reconstructing the velocity profiles. An important feature of the surface waves is their dispersive behaviour, or in other words, the dependence of their phase velocity with frequency. This feature allows the investigation of shallow and deep layers: the first are swept by the high-frequency content of the wavefield (e.g. generated by factories and traffic), while the second are sampled only by the long wavelengths (e.g. caused by oceans and winds). The maximum investigated depth is related with the sampled wavelengths and therefore with the array size. In site characterization and engineering studies, in general, the subsurface is investigated down to several hundreds of meters (e.g. 100-400 m). Chieppa et al. (2019) demonstrated the potential of ambient vibrations to investigate deeper structures, down to several thousands of meters.
Several single-station and array processing techniques were developed for site characterization studies. In general, they assume a homogeneous layered subsurface below the study area, investigate the geophysical properties and provide useful information for the P-and S-wave velocity profiles estimation. Single-station methods investigate the subsurface using the seismic recordings of three-Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00024-020-02516-x) contains supplementary material, which is available to authorized users. component sensors at each deployed site. The Horizontal-to-Vertical spectral ratio (H/V) is one of these techniques providing the resonance of the subsoil and the information concerning the subsurface velocity contrasts (H/V- Nogoshi and Igarashi 1971;Nakamura 1989); it is calculated using the spectra of the horizontal and vertical components together. Another single-station method is RayDec where the random decrement technique is used for the Rayleigh wave ellipticity computation (Hobiger et al. 2009); this technique enhances the Rayleigh wave contributions on the vertical and radial components, suppressing at the same time the Love and body waves. However, as a single-station technique, it does not provide information about the particle motion and direction of propagation.
Array processing methods, instead, evaluate the dispersion curves of the investigated area using the dispersive property of surface waves. Several techniques were developed in the last decades. Aki (1957) extracted the Rayleigh wave dispersion curves from the ambient noise vibrations using vertical-component sensors deployed in circular-shape arrays. Ling and Okada (1993) extended Aki's technique to arrays composed of seismic stations located on rings of increasing sizes; this method is known as Extended SPatial AutoCorrelation technique (ESAC). Bettig et al. (2001) implemented the original idea of Aki (1957) to imperfect circular arrays, making the deployment of arrays in urbanized areas easier. While the previous techniques analysed the vertical components alone, a big improvement was obtained by Köhler et al. (2007); the authors analysed the threecomponents of the seismic recordings and distinguished the Rayleigh and Love wave contributions. A different group of methods was suggested by Lacoss et al. (1969-classical beam-forming) and Capon (1969-high-resolution beam-forming). These methods identify azimuth and velocity of waves crossing the array by delaying and summing the signals according to the assumed wave speed and azimuth. Fäh et al. (2008) extended the analysis to the recordings of three-component sensors: the method, known as three-component high-resolution fk (3CFK), identifies the Rayleigh wave contributions on the vertical and radial components and the Love wave contributions on the transverse one. Poggi and Fäh (2010) extracted the ellipticity curves for the same picked vertical and radial dispersion curves. A more recent approach was suggested by Maranò et al. (2012); the technique, known as Wavefield Decomposition (WaveDec), extracts the Rayleigh and Love waves using a Maximum Likelihood approach. For each picked Rayleigh wave, the ellipticity angle is also estimated providing information about the sense of rotation.
The results of the single-station and the array analysis are inverted for the velocity profiles. As mentioned before, in site characterization studies the subsurface is assumed to be homogeneous and the inverted velocity profiles are usually attributed to the centre of the deployed array. The picked Rayleigh wave dispersion curves are inverted alone (Horike 1985;Tokimatsu et al. 1992) or jointly with the H/V curves to improve the constraint on the layers thickness Tokimatsu 2004, 2005;Picozzi et al. 2005;Parolai et al. 2005). In order to better constrain the properties of the shallow layers, Boxberger et al. (2011) inverted the H/V curves together with the Rayleigh and Love wave dispersion curves. Hobiger et al. (2013) investigated which parts of the Rayleigh wave dispersion curves and ellipticity curves should be jointly inverted for a successful inversion.
In this work, we show the results of a site characterization analysis performed in Switzerland within the Swiss Molasse basin. The seismic data were collected deploying three arrays of increasing size in and around the village of Schafisheim. The arrays, with minimum and maximum inter-station distances of 7.9 m and 26.8 km, were planned to investigate a wide frequency range (0.08-15 Hz) and reconstruct the velocity profiles down to several kilometres.
After a short geological and geographical introduction of the investigated area, we present the configurations of the three deployed seismic arrays and the analysis performed using two single-station methods and three array processing techniques. An overview of the results shows the picked dispersion curves and the ellipticity results together with the given mode interpretation. Using the dinver, toolbox of the Geopsy package (Wathelet et al. 2004) the velocity profiles were reconstructed down to 5 km by inverting the fundamental mode of Rayleigh and Love waves. Due to problems fitting some higher 4248 D. Chieppa Pure Appl. Geophys. modes, a detailed Horizontal-to-Vertical spectral ratio curve analysis was performed for the central part of the investigated area. A second inversion, subdividing the region in two subareas, was performed and the inverted velocity profiles (V P , V S ) were compared to the available geological and geophysical information.

Geographical and Geological Settings
The study site is located in and around the village of Schafisheim in the Canton of Aargau, in northern Switzerland (Fig. 1a). The area is located at the centre of the Swiss Molasse basin, a sedimentary basin between the Alps to the south and the Jura mountains in the north (Sommaruga et al. 2012). The basin was extensively studied in the past decades for hydrocarbon exploration, storage of nuclear waste and CO 2 sequestration (see Sommaruga et al. 2012, and references therein).
The oldest rocks in the area are Precambrian and represent the crystalline basement. In graben structures with north-east to south-west orientation, Carboniferous and Permian sediments were drilled by several boreholes. According to the available deep boreholes and interpreted seismic lines, the graben structures are mainly concentrated in the northern sectors. Above the Precambrian rocks, the Mesozoic units consist of shallow marine carbonates and marls (Triassic) and an alternation of limestones and marls (Jurassic and Cretaceous). The thickness of these units increases from 800 m in the north-east to 3 km in the south-west (3 km). The Tertiary units consist of Molasse deposits in sensu stricto subdivided in four lithostratigraphic groups: Lower Marine Molasse (UMM), Lower Freshwater Molasse (USM), Upper Marine Molasse (OMM) and Upper Freshwater Molasse (OSM). The deposition of these units occurred during a rapid deepening of the basin; it was filled first with turbiditic sediments (UMM), then with gravel fans (USM) and finally with marine sediments (OMM). Above, in the late-middle Miocene, the sedimentation became continental. The investigated area in Schafisheim was chosen for the availability of geological and geophysical information (drilling reports, seismic lines and 3D geological models) and the lack of a detailed site characterization measurement. At the beginning of the eighties, NAGRA, the Swiss technical competence centre in the field of deep geological disposal of radioactive waste, drilled the Schafisheim 1 borehole south-east of Schafisheim (Fig. 1b-Nagra Technischer Bericht 88-11, 1992); it drilled a geological sequence consisting of 220 m of Quaternary sediments, 340 m of Tertiary rocks and 930 m of Mesozoic rocks (Malm, Dogger, Lias, Keuper, Muschelkalk and Bundsandstein). The crystalline basement was found at 1490 m and consisted of granite (upper portion) and syenitemonzonite-diorite (lower portion) rocks (Fig. 1c).
Recently, a 3D model of the Swiss Molasse basin was developed (Allenbach et al. 2017). It was obtained by the re-interpretation of numerous deep boreholes and thousands of kilometres of seismic lines acquired inside the sedimentary basin and extending from Lake Geneva to Lake Constance. The geological interpretation given in GeoMol for the location of the borehole Schafisheim 1 interpreted the first 560 m as Lower Marine and Lower Freshwater Molasse sediments (USM) without taking into account the Quaternary sediments above (Fig. 1d, e). The Mesozoic units below, with the exception of the Bundsandstein unit, were identified in the same way as in the borehole.

Measurements
Three seismic arrays of increasing size were performed in 2016 and 2017. The centre of each array was located as close as possible to the Schafisheim 1 borehole location. Depending on the array size and on the frequency of interest, different recording times and sensor types were selected. We used both Lennartz 5 s and Trillium Compact 120 s sensors. The first were mainly used to investigate the shallow subsurface (hundreds of meters), while the second allowed us to measure at lower frequencies and to reach deeper structures. Both sensors were installed on the free surface; in soft soils, to have a better coupling, a small hole was dug in the ground (10-20 cm) and the sensor placed inside. The parameters of the three arrays are given in Table 1. Vol. 177, (2020) Ambient Vibration Analysis on Large Scale Arrays The small array was deployed in February 2017 and consisted of sixteen Lennartz 5 s sensors distributed as follows: one central station close to the Schafisheim 1 borehole and five rings of increasing diameter with three stations per ring. Each triangle was rotated clockwise with respect to the inner one, thus creating a spiral configuration (Fig. 2a). The minimum and maximum inter-station distances were 7.9 and 346 m, respectively, and the recording time was 120 min.
The intermediate array was deployed in October 2016. Due to the local topography, consisting of a hilly area with valleys oriented in north-south direction, the array configuration presents an oval shape (NNW-SSE), parallel to the extension of the valley. The minimum and maximum inter-station distances are 331 m and 2.4 kms, respectively (Fig. 2b). Twenty-four stations were simultaneously deployed at twelve selected sites: at each location, a Lennartz 5 s sensor and a Trillium Compact 120 s sensor were installed few meters apart one from the other. This choice was performed to test the operation of both sensors at low frequency, mainly as single stations, and to compare the resulting dispersion curves. Removing the first two hours of recordings because of the stabilization of the Trillium sensors, the total recordings were ten hours long; the recording was performed during the night.
The big array measurement was performed in July 2017; it consisted of one central station, 250 m southwest of the Schafisheim 1 borehole location, and 15 stations distributed in five triangles of increasing size. For the small array, each triangle was rotated with respect to the inner one, giving the array a spiral configuration (Fig. 2c); the minimum and maximum inter-station distances were 885 m and 22.8 km, respectively. The 16 installed sensors were all Trillium Compact 120 s and were deployed for 7 days. Due to technical field problems and GPS signal interruption, only 14 stations continuously recorded for 3 days; the other two worked properly for less than 24 h.
A spiral array configuration, as the one designed for the small and big arrays, presents a good compromise between good azimuthal coverage of the study area at a large number of investigated distance ranges and frequencies using a limited number of sensors (Kennett et al. 2015).
The seismic noise data were processed using single-station and array processing techniques. The data recorded using the big array were processed together with the recordings of the two stations DAGMA and EMMET of the Swiss Seismic Network. This increased the maximum inter-station distance to 26.8 km together with the array resolution limits.
The three array measurements were designed to cover a wide frequency range taking into account the available locations and the resulting array responses (bottom row of Fig. 2). These were computed using k min and k max values of the recording arrays obtained; k min is the diameter of the central lobe, while k max corresponds to the largest radius for which the theoretical array response (A) is lower than 0.5 ). The three deployed arrays have wide array responses (k max /k min ): 4.83 for the small, 5.92 for the intermediate (using Trillium sensors) and 3.26 for the big.

Single-Station Analysis
The recorded seismic data were analysed using two single-station techniques: Horizontal-to-Vertical spectral ratio (H/V-Nakamura 1989) and RayDec (Hobiger et al. 2009). The first method evaluates the resonance frequency of the subsoil, while the second directly estimates the Rayleigh wave ellipticity.

Horizontal-to-Vertical Spectral Ratio: H/V
The H/V curves were computed for all three arrays using 500 s window length, 5% window overlap and Konno-Ohmachi smoothing window with a b value of 40 (Konno and Ohmachi 1998). The H/V curves computed for the small, intermediate (Lennartz 5 s and Trillium Compact 120 s) and big arrays are manually picked and are shown in Fig. 3.
The H/V curves computed for the small array were analysed in the frequency range 0.2-20 Hz. They overlap perfectly over a wide frequency range and show small variabilities only above 7 Hz (Fig. 3a). All curves show a peak at 0.86 Hz and a second one for the inner stations above 15 Hz. The picking was performed analysing (1) the shape of the H/V curves (blue stars) and (2) the amplitude spectra of the three components (red stars) as shown in Castellaro (2016). The low-frequency peak was identified at the same frequency by both methods, while the high frequency one was identified by the blue stars at lower frequencies than for the red stars. The H/V curves for the intermediate array were computed between 0.01 and 20 Hz for both the Lennartz 5 s (Fig. 3b) and the Trillium Compact 120 s (Fig. 3c) sensors. In the first case, two peaks were identified at 0.79-0.99 Hz and between 9.9 and 18 Hz: the first was recognized at all sites by both picking methods, while the second was picked at three sites based on the H/V curve shape and at one taking into account the spectra. For the Trillium sensors, three peaks were identified using the H/V curve shape (0.09-0.1 Hz, 0.79-0.99 Hz and 7.1-9.7 Hz) and three using the spectra information (0.04-0.05, 0.79-0.96 and 9-18.2 Hz). Of the two peaks identified below 0.2 Hz, one was recognized using the H/V curve shape and the other using the amplitude spectra method. This is recognizable also in Fig. 3b even though with different amplitude values, below the eigenfrequency threshold of the instrument (red dashed line).
The H/V curves computed for the big array are shown in Fig. 3d. For almost all sites, taking into account the shape of the spectra, two peaks were identified at 0.1 and 0.9-1 Hz (red stars); at four sites a third peak was recognized (4.34, 5.45, 8.01 and 8.44 Hz). The second technique allowed the identification of one peak with different amplitude values at high frequency over a wide frequency range (0.67-15.3 Hz) and one peak at low frequency (0.02, 0.068 Hz) at just two sites. As shown in Fig. 3, the shape of the H/V curves provides information about the investigated subsurface and the main existing contrasts. At small scale, for example, the H/V curves overlap, meaning that the investigated subsurface is homogeneous; at larger scale, increasing the size of the array, the shallow subsurface becomes heterogeneous and the homogeneous portion, where all the H/V curves overlap, shifts to lower frequencies.
Below 0.15 Hz, the intermediate and big arrays$'$ H/V curves show large variability in terms of shapes and amplitudes related to the more heterogeneous subsurface or to the low energy contribution of the sources.
Comparing the H/V curves obtained using the Trillium Compact sensors (Fig. 3c, d) two differences and two common features were recognized: (1) the peak between 0.86 and 1 Hz disappears or is shifted to higher frequency moving away from Schafisheim, (2) high variability in terms of shape and amplitudes for the H/V curves below 0.2 Hz, (3) the red stars at 0.1 Hz identify an unclear peak and (4) a welldefined trough is always visible between 0.1 and 0.15 Hz.

Rayleigh Wave Ellipticity: RayDec
The second single-station method, RayDec, estimates the Rayleigh wave ellipticity curve using a random decrement approach (Hobiger et al. 2009). The Rayleigh wave ellipticity curves for the small, intermediate and big arrays are reported in Fig. 4. These were calculated between 0.2 and 20 Hz for the stations of the small array and between 0.01 and 20 Hz for the other two.
The ellipticity curves computed for the small array overlap over a wide frequency range (0.2-7 Hz); a wide peak (0.8-0.9 Hz) and a wide trough (1.5-3 Hz) were recognized. At 6-7 Hz all the ellipticity curves show a box-shape trend with amplitudes between 2 and 4; this peak (not visible in the H/V curves), generated by the increase of the power spectral density of the two horizontal components, could be linked to a factory nearby the investigated area. Above 7 Hz, the variability of the ellipticity curves increases and a second peak is identified at 11 sites between 10 and 20 Hz.
For the intermediate array, we also computed the ellipticity curves using data from both sensors. The shapes of the ellipticity curves computed using the two sensors are comparable down to 0.1 Hz, similar to the H/V curves. Also in this case, as previously noticed for the small array, between 6 and 7 Hz the ellipticity curves show a box-shape peak.
The ellipticity curves computed for the big array show good agreement between 0.1 and 0.9 Hz and 4254 D. Chieppa Pure Appl. Geophys. variability at lower and higher frequencies. While at high frequency, no specific trend can be recognized, a common trend made of descending flanks and plateaus can be identified at low frequency. Between 6 and 9 Hz, the three stations installed inside the village of Schafisheim itself show one peak or two narrow peaks similar to the observations for the other two arrays. This peak widens and shifts to higher frequencies with increasing distance from Schafisheim.

Array Analysis
The array analysis was performed using the three-  component alone and calculates the Rayleigh wave dispersion curves. For each array, a tuning phase was necessary in order to identify the best frequency and velocity ranges and parameters (e.g. time window length, window overlap, smoothing, etc.). The results reported below show the picked dispersion curves using the three above-mentioned techniques and the related standard deviation for the three arrays. For the intermediate array, the dispersion curves were computed using both sensors separately, but for the sake of simplicity only the curves calculated using the broadband sensors (Trillium Compact 120 s) are shown, as both sensors gave very similar results.

Three-Component High Resolution fk (3CFK)
The 3CFK method analyses the seismic recordings of the three components and calculates the dispersion curves for the transverse, vertical and radial components independently. The vertical component is used alone, while the two horizontal components are projected into possible azimuths with a 5°step in order to separate transverse and radial components. The possible azimuth with the highest energy of the respective component is identified. The analysis, taking into account the energy of the signal, allows the identification of Rayleigh waves on the vertical and radial components, and of Love waves on the transverse. For each wave arrival on the vertical and radial components, an ellipticity value is also derived (Poggi and Fäh 2010). Figure 5 shows the results for the small (top row), intermediate (central row) and big (low row) arrays for the transverse (left column), vertical (central column) and radial (right column) components.
For the small array, two modes were picked on the transverse component (Fig. 5a): the left one starts at 1.7 Hz close to the low-frequency array resolution limit and ends at 2.6 Hz; the right one ranges from 3.9 to 16.4 Hz. On the vertical and radial components (Fig. 5b, c), two similar curves were picked between 1.7 and 14.6 Hz.
Two dispersion curve branches were picked for the transverse component of the intermediate array between 0.5 and 3.6 Hz (Fig. 5d): the first mode extends from 0.5 to 3.0 Hz; the second starts at 1.55 and ends at 3.6 Hz, close to the upper-frequency array resolution limit. Two modes were picked for the vertical component (Fig. 5e): one from 0.7 to 2.4 Hz; the other, at higher velocities, covers the frequency range 1.5-2.7 Hz. An unclear and partially discontinuous mode was picked for the radial component ( Fig. 5f) between 0.63 and 1.7 Hz.
For the big array, two modes were picked for all three components. The first mode in the transverse dispersion panel extends from 0.076 to 1.01 Hz; the second one stretches between 0.35 and 0.74 Hz (Fig. 5g). For the vertical component, two modes were picked at 0.070-1.33 Hz and 0.28-0.91 Hz (Fig. 5h). The radial component presents two modes at 0.080-0.31 Hz and between 0.29 and 1.08 Hz (Fig. 5i).
The ellipticity curves derived for the vertical and radial components of the three arrays are reported in Fig. A1 of the Electronic Supplement. A similar trend is found for the ellipticity curves of vertical and radial components of the small array but not for the first modes of intermediate and big arrays. A moderate ellipticity peak is seen in plots C (small array) and E (intermediate array). The method does not allow us to retrieve the sense of rotation of the particle motion.

Wavefield Decomposition (WaveDec)
The second employed array technique is WaveDec (Maranò et al. 2012); it extracts the Rayleigh and Love wave information from the seismic recordings using a maximum likelihood (ML) approach. Together with the dispersion curves, WaveDec provides information about the Rayleigh wave ellipticity angle and the sense of rotation. The analysis can be performed using a pure Maximum likelihood approach (MIL), a Bayesian Information Criterion (BIC) or a mix of these two approaches. The chosen approach and the model complexity are controlled by a value called c, ranging from 0 (MIL) to 1 (BIC). For all arrays, an almost pure Maximum Likelihood approach (c = 0.1) was used and a maximum number of 3 modelled waves. In addition to c, the frequency range, the window length and the maximum number of simultaneously modelled waves control the dispersion curve estimation. The results are shown in 4256 D. Chieppa Pure Appl. Geophys. For the small array, one dispersion curve was picked for the Love waves between 3.75 and 14.8 Hz (Fig. 6a) and one for the Rayleigh waves from 1.6 to 13.9 Hz (Fig. 6b). The ellipticity angle derived from the Rayleigh wave dispersion curve is shown in Fig. 6g and presents negative values, denoting a retrograde particle motion.
For the intermediate array, in order to investigate lower frequencies, but higher velocities, a longer window length and a minimum velocity of 500 m/s were selected. The results for the intermediate array are shown in Fig. 6c, d, h, i. No dispersion curves were picked for the Love waves, while two Rayleigh wave modes were picked between 1.1-1.6 Hz and 1.8-1.92 Hz (Fig. 6d). Both ellipticity angles have ellipticity values around 0, making their particle motion definition unclear (Fig. 6h, i). For the big array, one curve was picked for the Love waves from 0.185 to 0.575 Hz (Fig. 6e) and two for the Rayleigh waves between 0.21 and 0.72 (Fig. 6f). The ellipticity angle associated with the first Rayleigh wave dispersion curve has negative values (j), while the second has positive ones (k), indicating retrograde and prograde particle motions, respectively.

Modified SPatial AutoCorrelation (MSPAC)
The Modified SPatial AutoCorrelation (Bettig et al. 2001) computes the Rayleigh wave dispersion curve using the vertical component of the seismic recordings. This technique calculates the cross-correlations between pairs of stations and groups them in rings of at least three station pairs, according to their interstation distance. The averaged cross-correlation functions of the rings are interpreted as Bessel functions which are inverted for the Rayleigh wave dispersion curve (Aki 1957). Differently from the original formulation of (Aki 1957), this technique does not request perfect circular array geometries, making the deployment and analysis easier. The analysis was performed using a frequencydependent window of 50 cycles (small and intermediate arrays) and 100 cycles (big array) and different frequency ranges depending on the array size. The results are shown in Fig. A2 of the in Electronic Supplementary material.
Using the vertical component of the seismic recordings, three Rayleigh wave dispersion curves were picked, one for each array. The curve picked for the small array extends from 1.02 to 10.88 Hz, the one for the intermediate array stretches from 0.52 and 1.63 Hz and the one for the big array covers the frequency range between 0.15 and 0.78 Hz.

Overview of the Results
The results of the different methods are shown in Fig. 7 for the Love and Rayleigh waves without showing the respective uncertainties. The Love wave dispersion curves were picked using 3CFK for all three arrays and WaveDec for the big and small arrays only (Fig. 7a). The curves picked using the two methods overlap over the entire frequency range (0.0075-15 Hz) with the exception of a small segment picked using WaveDec between 2 and 3 Hz. The Rayleigh wave dispersion curves computed for the three arrays are shown in Fig. 7b. They show perfect overlap for the small and intermediate arrays, while for the big one some discrepancies are seen. The results of WaveDec (solid light blue curves) have similar shape to the ones of 3CFK (solid black curves), but lower velocities. The curve picked using MSPAC (solid grey curve), instead, is located between two modes and can be interpreted as an apparent curve (Tokimatsu et al. 1992).
Based on this comparison, we can affirm that 3CFK and WaveDec methods performed in a similar way for all arrays and that the scatter of the WaveDec picked curves is generally higher than the one for 3CFK. The results of MSPAC are in agreement for the small and intermediate arrays but not for the big one where the inter-station distance increases and the signal correlation decreases. The ellipticity results are shown in Fig. 8 for the ellipticity curves (a) and the ellipticity angles (b). We compared the results obtained using the 3CFK and WaveDec array techniques and the RayDec singlestation method. Given the different ellipticity outputs obtained using the different processing techniques, we converted the WaveDec ellipticity angles into ellipticity curves using the equation where e is the Rayleigh wave ellipticity and n is the ellipticity angle (Maranò et al. 2012). In a similar way, using Eq. (1), the ellipticity computed using 3CFK (vertical and radial components) and RayDec were converted to ellipticity angles. Due to the lack of information concerning the sense of rotation, the 3CFK and RayDec ellipticity curves were calculated for positive and negative ellipticity angles. As seen for the dispersion curves, the different ellipticity curves and ellipticity angles show similar shapes over the entire frequency range. The ellipticity curves picked for the vertical components of the big and intermediate arrays match the corresponding RayDec ellipticity curves in the frequency range 0.2-0.95 Hz, while the results for the radial component of the same arrays overlap the RayDec curve of the small array. Above the H/V peak (0.8-0.9 Hz), the ellipticities computed using 3CFK are upward-shifted and the trough at 2.25 Hz is not well represented. The WaveDec information shows overlap with the Ray-Dec curve for the big array. The advantages of WaveDec are the information about the ellipticity angle and the particle sense of rotation. The dispersion curves picked using WaveDec (small and big arrays) show negative ellipticity angles and a retrograde sense of rotation, while the modes for the intermediate array present values around zero and a sense of rotation changing from prograde to retrograde.
In addition to the good agreement between the ellipticity curves and the ellipticity angles, the results in Fig. 8 show the advantage of using WaveDec compared to RayDec. Using WaveDec and the ambient seismic data recorded using the big array, we were able to pick two Rayleigh wave dispersion curves (Fig. 6f) with prograde and retrograde particle motions (Fig. 6j, k). Transforming these angles in ellipticity curves using Eq. (1), we see a similar trend between these two curves and the RayDec ellipticity curve for the big array. Based on our analysis, this similarity shows that the ellipticity curve computed using RayDec probably consists in a mix of several modes and not in a single mode.

Interpretation and Discussion
We inverted the Rayleigh and Love wave fundamental modes using the curves of Fig. 7. The fundamental mode of Love waves (L0) was defined taking into account the curve stretching from 0.077 to 3.06 Hz and computed using the three-component high-resolution fk; this was obtained by the match of the curves at low velocity picked using the three arrays (Fig. 7a). The corresponding mode for the Rayleigh waves (R0) consists of the curve picked for the big array (0.07-1.33 Hz), the mode picked for the intermediate array at lower frequency and the curve picked for the small array up to 14.6 Hz (Fig. 7b). Also in this case, the 3CFK curves were used for the inversion. The ellipticity angle derived from Wave-Dec for the fundamental mode shows a retrograde particle motion for the dispersion curves picked using the big and the small arrays, in agreement with our interpretation of fundamental mode. The particle motion for the left mode of the intermediate array is unclear and not considered in our interpretation.
The last step before the inversion concerns the maximum investigated depth. Taking into account the geometry of the big array, the minimum frequency of Rayleigh wave fundamental mode (0.070 Hz) and its maximum velocity (3690 m/s) and assuming to be able to investigate between 1/3 and 1/4 of the maximum wavelength, we estimate a theoretical penetration depth of around 14 km.
The inversion was performed using the fundamental modes of Love and Rayleigh waves as defined above, a velocity model consisting of 14 layers over the half-space and a maximum investigated depth of 5000 m. Based on the shape of the Rayleigh wave dispersion curve between 2 and 4 Hz (Maraschini and Foti 2010), low-velocity layers were allowed. The inversion was conducted using dinver, an open access package of Geopsy based on the Conditional Neighbourhood Algorithm (Wathelet 2008). The inversion algorithm creates an initial set of random velocity models and calculates the misfit of these models with the measured data using the formula where N is the total number of data points, D and M are the measured and modelled data points, respectively, and r i is the uncertainty of the measurement. i indicates the data points used for the inversion. In the neighbourhood of the best models, a new generation of random models is created. The misfit value is calculated using the square root formulation between the measured data and the modelled one.
The top row of Fig. 9 shows the results of the inversion. The velocity model with the lowest misfit (best model) and the corresponding dispersion curves are shown in grey. The bottom row reports the comparison between the generated dispersion curves (orange lines) using the velocity model with the lowest misfit and the picked dispersion curves (bottom row-left and central panels) and the ellipticity angle curves (bottom right panel).
The modelled Rayleigh and Love wave fundamental mode curves are in agreement with the measured dispersion curves. Some disagreements are found with the picked curves between 1.5 and 3.7 Hz for the Love waves ( Fig. 9-bottom row left) and between 1.5 and 2.6 Hz for Rayleigh waves (Fig. 9bottom row centre).
To better understand the investigated area and to improve the fit of the dispersion curves, the H/V curves were analysed. For the small array (Fig. 3a), the curves overlap pretty well from 0.2 to 7.5 Hz. Above, as seen by the H/V peak variability, some local heterogeneities occur in the shallow subsurface. The shape and the amplitude of the H/V curves calculated for the intermediate array are comparable down to 0.1 Hz (Fig. 3b, c) and good overlap is seen in the frequency range 0.1-0.5 Hz. At higher and lower frequencies, the curves show higher variability; the high frequency variation is related to the subsurface variation from site to site, while the low frequency one, as mentioned in 3.1.1, is influenced by the strong noise increase at low frequency, by the installation and the short recording times. Above the H/V peak at 0.86 Hz, two different trends were recognized: a steep one and a gentler one. The H/V curves picked for the big array show high variability at low (left of 0.1 Hz) and high (right of 0.7 Hz) Vol. 177, (2020) Ambient Vibration Analysis on Large Scale Arrays 4261 frequencies and pretty good overlap between 0.1 and 0.7 Hz, as seen for the intermediate array. The highfrequency scatter explains the local variability from site to site; at low frequency the H/V curves have plateaus with much lower values compared to the curves of Fig. 3c.
With the identification of two trends above the 0.86 Hz peak for the H/V curves computed using the intermediate array and to understand if any spatial distribution exists, we performed a detailed H/V analysis in the centre of the investigated area. 40 additional points were measured following a regular grid from north to south. The deployed sensors, recording the seismic ambient vibrations from 45 min to two hours, were Lennartz 5 s sensors installed on the free surface using metallic tripods. The computed H/V curves were analysed calculating the H/V slope angle between the peak and the respective trough in the frequency range 0.8-3 Hz. The calculated H/V slope angles are reported in Fig. 10: red points indicate steep H/V flanks; yellow points gentle flanks.
The information provided by the H/V curve slope angles (Fig. 10a) and the multilevel b-spline interpolation (Fig. 10b) show the co-existence of a central area with a main north-south orientation, where the H/V curves present a steep flank and a narrow trough (red points), and two surrounding zones to the east and west with gentler flanks and wider troughs (yellow points). The hilly topography of the area and the geological map support this observation (Fig. 10c). The mapped geological units show the transition from coarse gravels and glacioalluvial sediments in the central and eastern areas to a thin cover of unconsolidated deposit on the western side. The transition to areas with different morphological and geophysical properties can explain the impossibility to simultaneously fit Love and Rayleigh wave dispersion curves over the entire frequency range and the three arrays.
Based on this idea, we assume that the properties of the Rayleigh and Love waves differ in the central and surrounding areas and that we therefore measured at the same time two different fundamental modes for Rayleigh and two for Love waves. To verify our hypothesis, we performed two independent inversion runs using similar starting models and different dispersion curves. The new mode interpretation is reported in Fig. 11: L0/L1 and R0/R1 labels report the dispersion curves common to both zones and attributed to the Love and Rayleigh waves, respectively, and the subscripts inner and outer are used to refer to the central and to the surrounding areas, respectively. The inversion was performed using the fundamental and first higher modes of Rayleigh and Love waves and the negative ellipticity angle derived using the RayDec ellipticity curves for the intermediate array, assuming a retrograde particle motion. The resulting ellipticity angles with retrograde particle motion show a steep flank for the inner area and a less inclined trend in the outer two between 0.86 and 3 Hz. The locations of the two seismic stations used in the inversion are marked by black circles in Fig. 10b.
The inversion results are shown for the inner (Fig. 12) and the outer areas (Fig. 13). While the number of layers (14) and the maximum investigated depth (15 km) were equal for both inversions, a lowvelocity zone was allowed only in the shallow layers of the inner zone.
A good fit was achieved for both inversions. For the inner area (Fig. 12), the best velocity model, represented in dark grey, explains the inverted dispersion curves and the ellipticity angle. Small disagreements exist with the Rayleigh wave fundamental mode at low frequency and with the Rayleigh wave first higher mode between 0.7 and 0.9 Hz; this has similar shape as the theoretical dark grey curve, but it is shifted towards lower velocities. The use of low-velocity layers allowed a good fit of the Rayleigh wave fundamental mode between 1.5 and 4 Hz. Also for the outer zones (Fig. 13), the velocity model with the lowest misfit (dark grey) is in good agreement with the dispersion curves and the ellipticity angle. As for the inner area, small disagreements can be seen with the Rayleigh wave fundamental mode at low frequency. The best P-wave velocity profile and the ten best S-wave velocity profiles for each investigated area are plotted in the top row of Fig. 14; good agreement can be seen down to 15 km for the S-wave velocity profiles of the two zones. The plot in Fig. 14b shows a zoom of the first 3500 m and the comparison with the measured P-wave velocity profile (thin light grey curve-Nagra Technischer Bericht 88-11, 1992) and its interpretation as reported in the Seismic Atlas of the Swiss Molasse basin (black curve- Sommaruga et al. 2012). Fig. 14c shows the interpreted geological sequence of the Schafisheim 1 borehole (Nagra Technischer Bericht 88-11, 1992).
The inverted P-and S-wave velocity profiles allowed the identification of the bottom of the Quaternary sediments, the transition between Tertiary and Malm rocks and a deep contrast at 2000 m (Fig. 14b, c). The first contrast at around 215 m, responsible of the H/V peak at 0.86 Hz, is well marked in all velocity profiles and confirmed by the geological information reported in the Schafisheim 1 borehole report (Nagra Technischer Bericht 88-11 1992). The second contrast is located at around 590 m and is well visible in the P-wave velocity profile of the outer area as a smooth increase in the S-wave velocity; the third one, at around 2000 m, is recognizable by a velocity increase in all S-wave profiles. At higher depths, using the available data, it is not possible to clearly see the top of the crystalline basement in the velocity profiles at the same depth as the borehole interpretation. In the first 3.5 km, some variability in terms of thickness and velocity are observable for the P-and S-wave velocity profiles. While for the outer area a velocity profile with a gradient was found, for the inner zone a low-velocity zone was derived in the first three layers. It was identified between 50 and 150 m where, based on the geological and geophysical data of the Schafisheim 1 report, an alternation of sands and gravels of Quaternary age exists.

Figure 11
New mode attribution for the Love (top) and Rayleigh (bottom) wave dispersion curves. The L0/L1 and R0/R1 labels refer to the dispersion curves portions common to both investigated areas, while the words inner and outer explain to which area each portion belongs to. The high-resolution fk results are in black (thick), the WaveDec results in blue, the MSPAC results in grey and the array resolution limits in black (light dashed lines) c Figure 12 Inversion results for the inner area. Top row: fundamental and first higher modes of the Rayleigh wave dispersion curves and Rayleigh wave ellipticity angle. Central row: fundamental and first higher modes of Love wave dispersion curves. Bottom row: P-and S-wave velocity profiles and zoom to the first thousands of meters for the P-and S-wave velocity profiles. The model with the lowest misfit is shown in grey

Conclusion
In this work, we show the results of an array analysis performed in the Swiss Molasse Basin. The recorded seismic data were processed using singlestation and array processing techniques and inverted for the P-and S-wave velocity profiles estimation. The first inversion produced a good fit among the picked fundamental modes of Rayleigh and Love waves and the theoretical ones together with a remarkable disagreement with all the other modes. A closer look at the H/V curves computed for the small and intermediate arrays highlighted two different flanks above the 0.8-0.9 Hz H/V peak. These curves and additional single-station measurements confirmed the existence of a north-south elongated central area, surrounded by two zones with different properties (Fig. 10). Based on this information and the geological map, we hypothesized the occurrence of lateral variations at the centre of the investigated area and the simultaneous presence, in the dispersion curve plots, of two different fundamental modes for the Rayleigh and two for the Love waves. Using the new observations, a new mode interpretation was provided and two inversions were performed: one for the inner area and one for the surroundings. In both cases, the inversions were performed inverting fundamental and higher modes of Rayleigh and Love waves together with the ellipticity angles. The inverted P-and S-wave velocity profiles allowed the reconstruction of the subsurface down to 15 km. The investigated site shows the importance of a joint H/V and dispersion curve interpretation, as well as the importance of a reasoned modal attribution for the extracted dispersion curves.
The tested array processing techniques show comparable results at low and high frequency for the Love and Rayleigh waves. Based on the results of our analysis, 3CFK shows continuous dispersion curves and computes the dispersion curves for all arrays, even if lateral variations occur. For these reasons, we used the dispersion curves derived from 3CFK for the inversion; the other curves, even if providing similar information, were not chosen for their scatter and less continuity over the analysed frequency range (e.g. WaveDec) but were used to confirm the results of 3CFK. For the ellipticity information, RayDec and WaveDec results have to be considered together. WaveDec provides the information about the particle motion, while RayDec gives an estimation of the ellipticity over a wide frequency range. For these reasons, the ellipticity interpretation should be performed using the results of both methods together.
As mentioned at the beginning of this work, two seismic stations were installed at each of the twelve sites of the intermediate array, close to each other. This choice was performed to test the sensor b Figure 13 Inversion results for the outer area. Top row: fundamental and first higher modes of the Rayleigh wave dispersion curves and Rayleigh wave ellipticity angle. Central row: fundamental and first higher modes of Love wave dispersion curves. Bottom row: P-and S-wave velocity profiles and zoom to the first thousands of meters for the P-and S-wave velocity profiles. The model with the lowest misfit is shown in grey performances at low frequency and to check the advantages in using broadband sensors. Based on the deployed array configurations, we compared the sensor performances only in terms of single-station methods and conclude that the Lennartz 5 s sensors perform properly down to 0.1 Hz. In order to compare the performances of the short-period and broadband sensors in terms of dispersion curves, a larger array should be deployed. This work confirms the potential of surface waves and proves the possibility to apply the same methodologies developed for microzonation studies to crustal investigations.
Inverting the available dispersion curves and ellipticity angles, the maximum investigated depth is 15 km in agreement with our initial estimation. Due to the increase of the uncertainty with depth, the use of additional data (e.g. geological or geophysical measurements) and models is strongly recommended.