Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings: Application to Seismic Stations in Switzerland

We performed a study to test the performances of the Hybrid Heuristic Search (HHS) algorithm (Nagashima et al. 2014) using earthquake recordings at 23 instrumented sites in Switzerland. The HHS algorithm is based on the diffuse field theory and estimates the P- and S-wave velocity profiles by inverting horizontal-to vertical spectral ratio (H/V) curves of earthquake recordings. In this study we inverted the H/V curves computed using the arrivals of local and regional earthquakes with a maximum local magnitude of 4.6. We compared the results of the inversion technique to the information available for each investigated site (e.g., site characterization analysis). Our results emphasize the importance of site characterization data for the definition of the parameter space in the near surface, and the potential to extend the investigated depth to much higher depths than the standard site characterization measurements. The additional constraints in the inversion are important to drive the inversion towards a unique solution. Our analysis confirms the potential of the HHS algorithm to invert the full earthquake H/V curve for subsurface investigation and its applicability to areas with low magnitude earthquakes. In addition, the HHS algorithm can be a useful technique to investigate deeper structures and extend the velocity profile to depths that are not resolved by classical site characterization techniques.


Introduction
Passive seismic methods are a group of geophysical techniques developed to investigate the subsurface structure and estimate the P-and S-wave velocity profiles without actively generating signals. The data used by these methods consists of ambient seismic vibrations, which are also known as seismic noise or microtremors and are a mix of mainly Rayleigh and Love waves propagating along the topographic surface. Ambient seismic vibrations have the advantage of being continuously recorded by temporary and permanent installations from less than one hour to several days and are used in site characterization analysis to estimate the elastic properties of the subsurface (e.g. Hobiger et al., 2021;Michel et al., 2014;Picozzi et al., 2009). Earthquake recordings, instead, consist of body (P-and S-waves) and surface waves. The body waves propagate from the hypocenter in all directions following a spherical pattern; these waves, at the transition with the topographic surface, are partly converted to surface waves and start to propagate parallel to the surface in a cylindrical way and are partly reflected downwards. While the energy of earthquake recordings is several orders higher than that of ambient seismic vibrations, the disadvantages of earthquake data are the restricted use in areas with low or no seismicity and their infrequent occurrence so that the detection becomes only possible using permanent networks. Applications of earthquake recordings are in tomographic studies (e.g. Diehl et al., 2009;Haslinger et al., 1999;Husen et al., 2000), for the calculation of receiver functions (e.g. Bertrand & Deschamps, 2000;Dugda et al., 2005;Kummerow et al., 2004;Lombardi et al., 2008;Ozacar & Zandt, 2004;Zor et al., 2003), and Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/s00024-022-03226-2. the computation of the amplification functions (e.g. Edwards et al., 2013;Michel et al., 2014).
In this work, we show results obtained at many sites in Switzerland for the Hybrid Heuristic Search (HHS) inversion algorithm . The algorithm inverts earthquake H/V curves using an initial parametrization and estimates the P-and S-wave velocity models for each site. The inverted H/V curves are computed combining the Fourier spectra of horizontal and vertical components of an earthquake recording starting at the S-wave arrival.
For the examples shown in this work, we propose a solution that combines the horizontal-to-vertical spectral ratio (H/V) curves from earthquakes with the results from site characterization analyses performed near each site using ambient seismic vibrations. To proceed with the inversion and estimate the P-and S-wave velocity profiles, the H/V curves from earthquakes are calculated. Different authors compute the H/V curves from earthquakes by selecting the whole signal duration or shorter windows. Depending on the length of the selected signal, P-, S-, or coda waves are included in the analyzed signal and used for the H/V computation (e.g. Horike et al., 2001;Lachet et al., 1996;Lermo & Chávez-García, 1993). The H/V curves computed in this work are obtained for the signals after the S-wave arrival in order to select the S-wave and the coda information. The H/V curves computed for each individual event are then averaged and inverted for the velocity profiles of the site. The HHS method, previously tested by Nagashima et al. (2014) at 5 sites in Japan, allowed the estimation of S-wave velocity profiles down to the bedrock. The data used by these authors consisted of 77 regional earthquakes with magnitudes (M JMA ) between 4.0 and 7.3 and epicenters localized in a restricted area, east of the study area. The assumption on which the inversion algorithm is based is a diffusive wavefield after the arrival of the S-wave; this condition, based on the studies of Kawase et al. (2011), is achieved by averaging the S-wave arrival of many events over various incidence and azimuthal angles.
In this study, (1) we test the applicability of the HHS method in Switzerland, a country where the seismicity is much lower than in Japan. The local magnitude of the analyzed earthquakes in our study ranges between 2.5 and 4.6.(2) The epicenters of the analyzed earthquakes are not one-directional as in the work of Nagashima et al. (2014) but distributed around each permanent installation to investigate the properties of the subsurface from all directions. In a preliminary phase, we tested the HHS algorithm at station SBERN for the same input data (average H/V curve) and two different parameterizations consisting of completely free layers and of a mix of partially constrained and unconstrained layers. The motivation behind this preliminary phase relates to the results obtained by Nagashima et al. (2014) that set the Pand S-waves of the half-space as the only constraint. Our results show the importance of constraints for the definition of the parameter space, the possibility to improve the resolution of the subsurface by adding some variability to the parameters of the shallow layers and the potential to extend the investigated depth at much higher depths than the standard site characterization measurements. Based on the findings for the station SBERN, we extend the analysis to 23 sites in Switzerland using the parameterization consisting of a mix of partially constrained and free deep layers. The inversion results for three of the 23 stations are shown in detail in the section Inversion results, while the remaining results can be found in the Electronic Supplementary Information. In the section Discussion, we illustrate in detail the features and the recurrent trends identified for the measured curves at different sites and group the analyzed stations in three classes. Finally, we summarize our findings, from the choice of the initial parameterization to the analysis performed at all sites.
Our choice to constrain the shallow layers of the initial model to the velocity profiles estimated by the array measurements is driven by the possibility to use geophysical information collected for the study area, to reduce the space for the parameter search, and to decrease the computation time. At each of the 23 selected sites, passive and active array measurements were performed by the Swiss Seismological Service (Michel et al., 2014) or by a private company (Fäh et al., 2009), respectively. The array measurements performed at each site provide a different number of Rayleigh-wave ellipticity curves and/or dispersion curves that are jointly inverted for the P-and S-wave velocity profiles. In order to compare the final results of the HHS inversion with the available data and to calculate the misfit, synthetic dispersion curves (Wathelet et al., 2020) and the theoretical SH-wave transfer functions (Knopoff, 1964) were calculated. While the first were compared with the measured dispersion curves, the others were corrected with respect to the Swiss reference rock velocity profile (Poggi et al., 2011) and then compared with the empirical amplification function for the elastic case computed by the Empirical Spectral Modelling method (ESM-Edwards et al., 2013).
In the following, we present the developed and applied processing methods and the inversion techniques for site characterization analysis used by the Engineering Seismology group of the Swiss Seismological Service (http://www.seismo.ethz.ch/en/ home/). The analyzed data consists of Rayleigh and Love waves with negligible contributions of body waves recorded by single sensors or by arrays of sensors deployed on the topographic surface following ad-hoc configurations. Three-component sensors, usually Lennartz 5 s sensors, were used for seismic passive acquisition, while vertical 4.5 Hz geophones were used for active measurements. Combining the Fourier spectra of the two horizontal components and the spectrum of the vertical component, the fundamental frequency of resonance (f 0 ) of the sediments is retrieved; it is inversely proportional to the depth of the bedrock and gives an indication of the sediment thickness. The techniques used for this purpose are the horizontal-to-vertical spectral ratio curves (H/V- Nakamura, 1989) and the Rayleigh-wave ellipticity decrement (RayDec-Hobiger et al., 2009). The array methods, instead, exploit the information collected by several sensors deployed on the topographic surface and simultaneously recording the seismic wavefield propagating through the array. Depending on the type of the deployed sensors, the number of components and the chosen array configuration, different processing techniques were developed to extract the dispersion curves (e.g. Aki, 1957;Bettig et al., 2001;Capon, 1969;Lacoss et al., 1969). The first techniques focused on the analysis of the vertical component alone and estimated only the Rayleighwave dispersion curves. With the development of three-component sensors and their diffusion, new techniques were introduced to extract the dispersion curves of Love and Rayleigh waves (e.g. Köhler et al., 2007). The techniques used for site characterization analysis employed by the Swiss Seismological Service are 3-component high-resolution f-k (3CFKFäh et al., 2008;Poggi & Fäh, 2010), Wavefield Decomposition (WaveDec- Maranò et al., 2012Maranò et al., , 2017 and Modified Spatial Auto-Correlation (MSPAC-Bettig et al., 2001). The 3CFK method computes the dispersion curves for the vertical, radial, and transverse components individually. Poggi and Fäh (2010) extracted the dispersion curves for the vertical and radial components and the corresponding ellipticity curves. The second technique, also known as WaveDec, uses the information of all components together and estimates the Love-and Rayleigh-wave dispersion curves by a maximum likelihood approach. For each frequency of the measured Rayleigh-wave dispersion curve, the ellipticity angle curve is computed. It describes the shape of the ellipse along which the particles of Rayleigh wave move together with the rotation sense. The absolute value of the ellipticity angle curve, defined as the ratio of horizontals and vertical axis, is equal to the tangent of the ellipticity (Maranò et al., 2017). The last method is an implementation of the SPAC technique (Aki, 1957) and is based on the use of irregular rings defined by minimum and maximum radius (Bettig et al., 2001). While this technique is implemented in Geopsy for three components (Wathelet et al., 2020), we generally use the method for the vertical component only.
At certain sites, where the very shallow portion of the underground is of major interest, i.e. sites with relatively high resonance frequencies, active seismic techniques are employed and the data are analyzed by the Multi-channel Analysis of Surface Waves technique (MASW-Park et al., 1999;Xia et al., 1999). Seismic waves are here generated by an active source (e.g. sledgehammer or a weight-drop) located at the ends of the deployed geophone line, which is constituted by equally spaced vertical or threecomponent geophones recording the propagation of the generated Rayleigh and/or Love waves through the entire deployment. The seismic recordings of the passive measurements are analyzed for the computation of the dispersion curves for the high-resolution f-k method implemented in the Geopsy software. Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings 757 Finally, the results of single station techniques and array methods are interpreted and inverted using the Conditional Neighbourhood Algorithm implemented in the dinver toolbox (Wathelet et al., 2008) of the Geopsy package (Wathelet et al., 2020). This technique allows the estimation of P-and S-wave velocity profiles. Whenever possible, a joint inversion combining the dispersion and the ellipticity curves is performed to obtain the model with the best fit. Several authors performed joint inversions combining the Rayleigh and/or Love wave dispersion curves alone (Arai & Tokimatsu, 2004;Boxberger et al., 2011;Horike, 1985;Tokimatsu et al., 1992), the dispersion curves and the H/V curves (Arai & Tokimatsu, 2005;Parolai et al., 2005;Picozzi et al., 2005), or a combination of Rayleigh-wave dispersion curves and Rayleigh-wave ellipticity curves (Hobiger et al., 2013), or the ellipticity angle curves (Maranò et al., 2012(Maranò et al., , 2017. Ellipticity and H/V curves better constrain the layer thickness, while the dispersion curves determine the absolute velocities of the layers (Scherbaum et al., 2003).
The array analyses, as performed by the Swiss Seismological Service, end with the computation of the theoretical SH-wave transfer function (Roesset et al., 1970) of the estimated velocity profiles. The SH transfer function shows the changes of amplification with frequency and corresponds to the amplitude ratio between the SH-wave at the surface and the amplitude of the same wave at the bedrock. The theoretical SH-wave transfer function is computed for the S-wave half-space velocity estimated by the site characterization analysis and compared with the empirical amplification function (ESM, Edwards et al., 2013) computed by the Swiss Seismological Service. The ESM curve provides the site amplification at instrumented sites using the earthquake recordings over a wide range of magnitude. The signal used for the Fourier spectra computation consists of the S-wave and coda information from single events. The noise recording is extracted before the P-wave arrival and is artificially increased to intersect the Fourier spectra of the signal and avoid the spectral modelling of noise in case of small and moderate earthquakes (Edwards et al., 2013). To compare the SH-wave transfer functions with the empirical amplification function of each site, both curves are referenced to the Swiss rock velocity profile (Poggi et al., 2011), a gradient velocity profile with shearwave velocities increasing from 1000 m/s at the surface, to 3200 m/s at 4 km depth. The Swiss rock velocity profile was obtained combining the measured velocity profiles and the amplification functions at several instrumented sites in Switzerland. The ESM technique overcomes the installation need of having two seismic stations recording at the same time as for the standard spectral ratio (SSR- Borcherdt, 1970) and the surface-to-borehole spectral ratio (Knopoff, 1964) techniques. These techniques compute the amplification of the seismic waves in sedimentary layers by deploying one seismic sensor on sediments and one on rocks. The SSR uses two sensors which are located on the surface within the same region, while the surface-to-borehole spectral ratio technique has one sensor located at the surface and the other buried underneath, at the transition between the half-space and the sedimentary basin.

Data
The Swiss Seismic Network consists of more than 200 seismic stations distributed over the entire country with higher density in areas of higher seismicity (e.g. Canton of Valais, area around Basel, and Canton of Grisons). The network includes shortperiod, broadband, and strong-motion stations recording local and global earthquakes over a wide range of magnitudes.
In this work, we collect and process the earthquake information recorded at 23 sites of the Swiss Seismic Network from the respective installation date, different for each site, to September 2018. The stations are all located on soil. Three are equipped with short-period sensors, six with broadband sensors and 14 are strong-motion stations. The analyzed earthquakes consist of local events with epicenters in and around Switzerland and local magnitudes ranging between 2.5 and 4.6. The information on stations and sensor type, installation date, number of analyzed earthquakes, minimum and maximum epicentral distances, and frequency range of the average H/V curve are reported in Table 1. The seventh column of Table 1 shows the thickness of the Molasse Basin deposits as reported in Geomol website (https:// viewer.geomol.ch, Allenbach et al., 2017). The Geomol model is a 3D representation of the geology north of the Alps, in the sedimentary basin known as the Swiss Molasse Basin. The model was built combining many datasets including seismic lines, borehole logs, and cross-sections. The last two columns of Table 1 provide information about the Vs 30 and the type of site characterization measurement performed.
The 23 sites were selected taking into account the availability of site characterization measurements at the installation site, the type of installed seismic sensor, and the availability of published reports (Fäh et al., 2009;Michel et al., 2014;Poggi et al., 2017).
The Engineering Seismology group of the Swiss Seismological Service (SED) performed measurements at 18 of these sites using passive (17) and active (1) methods. The passive acquisitions were performed with Lennartz 5 s sensors arranged in an array configuration consisting of rings of increasing size around a central station (type 1 in Table 1). The recording time varied between 2 and 3 hours (Hobiger et al., 2021). In Lausanne at the EPFL campus (station SEPFL), an active measurement was performed using the MASW technique. The installed sensors were three-component geophones with a 4.5 Hz corner frequency deployed in a straight line of 46 m length (type 2 in Table 1). The chosen source was a 120 kg dropping mass device installed at both extremities of the deployment. The depth estimated by these measurements ranges between 38 and Table 1 Station name, sensor type (HG: strong-motion sensor; HH broadband sensor; EH short-period), installation date, number of picked S-wave arrivals, epicentral distances, H/V frequency range, thickness of the Molasse sediments as estimated by the GeoMol project (Allenbach et al., 2017) 265 m. At all sites, Rayleigh-wave dispersion curves and H/V curves were measured; at 16 of the sites, Love-wave dispersion curves were also measured. The description of the analysis performed at these sites can be found in Fäh et al. (2009), Michel et al. (2014 and Poggi et al. (2017). At the remaining 5 sites, a private company performed MASW, shear wave refraction tomography, and/or P-wave hybrid seismic measurements (refraction and reflection seismic data); these sites are classified as type 3 in Table 1. The source used for the active measurements is an 8 kg sledgehammer and the receivers are horizontal and/or vertical geophones. The active measurements performed at each of these sites are described in Fäh et al. (2009) and summarized in Poggi et al. (2011). The maximum investigated depth at type 3 sites is 50 m, which limit the exploration to the shallowest layers. In some cases the measured dispersion curves were not retrieved from the MASW acquisition, so that the velocity profiles were derived from shear-wave refraction tomography (e.g. sites WEIN and WILA). The MASW-derived Rayleigh-wave dispersion curves at FLACH, STEIN and TORNY were digitized from the station report provided by the private company.
Three types of sensors (short-period, broadband, and strong motion sensors) are installed at the 23 sites of interest. Based on the sensor type and on the noise level, the number of recorded seismic events differs: broadband and short-period sensors have a high number of recorded earthquakes both at low and high frequencies, while the strong-motion sensors have many recordings only at high frequency due to instrument noise. The sampling frequencies of the selected sensors are 200 or 250 Hz, whereas the antialias filter is set by the hardware.
The results of site characterization measurements performed by the Swiss Seismological Service and by the private company are stored in the site-characterization database of the Swiss Seismological Service (Swiss Seismological Service (SED) at ETH Zurich, 2015). These include microtremor H/V curves, Rayleigh/Love wave dispersion curves, Rayleighwave ellipticity curves, P-and S-wave velocity profiles, and in some cases also geotechnical acquisitions (e.g. CPT and sCPT campaigns). For the estimated velocity profiles, the theoretical SH-wave transfer function is also computed (Roesset et al., 1970) for all analyzed sites. The site characterization data plays a fundamental role in this work: the estimated velocity profiles from site characterization analysis define the initial parameters for the shallowest layers in the initial velocity profiles which will adjust during the inversion. The dispersion curves and the transfer functions, not directly used in the HHS inversion, are compared with the synthetic curves from the inverted velocity profiles.
The following paragraphs present the geological and geophysical information available for the seismic stations SBERN, SYVP, SARK and STIEG in detail. These sites will be used to highlight the results of the HHS inversion algorithm and the final classification. The remaining sites are shown in the Electronic Supplementary Information.
Station SBERN is located in the city center of Bern (Canton of Bern); it has an Episensor ES-T sensor operating since June 2013. From a geological point of view, the site is located on fluvioglacial and glaciolacustrine gravels. Around the seismic station, a passive array measurement was performed (Michel et al., 2014). The seismic recordings of two arrays, with minimum inter-station distance of 10 m and maximum inter-station distances of 100 and 200 m, respectively, consisting of 16 sensors each were processed for single-station and array processing techniques and allowed the identification of the fundamental mode dispersion curves for Rayleigh and Love waves. The inversion results shown by Michel et al. (2014) allowed the estimation of the Pand S-wave velocity profiles down to 145 m using the two dispersion curves (2.5-14.5 Hz) and the right flank of the Rayleigh-wave ellipticity curve (2.5-3.5 Hz) computed with the time-frequency analysis for the central seismic station of the array. The estimated velocity profiles present two interfaces with mild velocity contrasts at 40 and 120 m. The first one marks the transition between the fluvial sediments above, and the Molasse deposits below, and shows an S-wave velocity increase from 400-500 to 800 m/s; the second interface corresponds to the half-space where the S-wave velocity increases to 1000-1500 m/s at a depth of 160 m. SYVP is located on lacustrine deposits in the town of Yverdon (Canton of Vaud). The installed sensor is an accelerometer (Episensor ES-T) operating since September 2011. The subsurface around the installed seismic station was investigated using passive seismic ambient vibrations by means of two seismic arrays with equal minimum inter-station distances of 15 m and apertures of 120 and 240 m, respectively (Michel et al., 2014). For each array, 16 sensors were installed on metallic tripods. The array analysis allowed the picking of three dispersion curves for the Rayleigh waves and three for the Love waves in the frequency range between 1.14 and 8.2 Hz. The observed curves, whether for Rayleigh or Love waves, were interpreted as fundamental mode, first and second higher modes with the fundamental mode having the lowest respective velocity and the second higher mode the highest. The velocity profiles were estimated down to 140-150 m by jointly inverting the three Rayleigh-wave dispersion curves and the ellipticity curves, picked over the same frequency range as the Rayleigh-wave dispersion curves. The velocity profiles for station SYVP do not present strong impedance contrasts. Above 6 m, the velocity profiles show two thick layers with similar velocities and no big velocity contrasts. In contrast, below 60 m, the velocity profiles present a linear gradient with seismic velocities increasing up to 2440 m/s at about 140 m depth.
SARK is installed in Sarnen (Canton of Obwalden), on a lacustrine plain. The sensor is an accelerometer (Episensor ES-T) and has been in operation since May 2012. Around the permanent seismic station, two passive seismic array measurements of increasing size were performed (Michel et al., 2014) using 16 Lennartz 5 s sensors. The minimum inter-station distance for both arrays is 15 m, while the apertures are 100 and 200 m, respectively. The array processing analysis produced four dispersion curves: two for the Rayleigh waves and two for the Love waves. The modes measured for the Rayleigh waves (2.6-17.6 Hz) were interpreted as fundamental mode and first higher mode. The two modes for the Love waves (2.7-19.3 Hz) were attributed to the fundamental and first higher modes. The velocity profiles were estimated down to 150-170 m inverting the four dispersion curves, the Rayleigh-wave ellipticity curve (2.4-3.8 Hz) was computed with the three-component high-resolution f-k technique (Poggi et al., 2011) and the ellipticity peak at 2.3 Hz. The estimated velocity profiles show two interfaces: one at about 30 m with Vs of 1450 m/ s and one at the transition to the half-space where the S-wave velocities increase up to about 3000 m/s. STIEG is a borehole installation located in Stiegenhof (Canton of Zurich). The station, operating since December 2012, is equipped with a Lennartz 3D-BH sensor 123 m below the surface and an Episensor ES-T sensor at the surface. In January 2013, 50 m south-west of the seismic station, two passive site characterization measurements consisting of two circular arrays (16 sensors) with radii of 100 and 200 m around a central station were performed (Poggi et al., 2017). The array data analysis allowed the measurement of Rayleigh-and Love-wave fundamental modes between 2.1 and 17.7 Hz. The measured dispersion curves and the resonance peak (f 0 ) at 4.4 Hz obtained from the H/V spectral ratio curves allowed the estimation of the P-and S-wave velocity profiles down to 160 m. The estimated S-wave velocity profile shows a linear gradient with depth down to 60 m. At higher depths two velocity contrasts (60 and 160 m) can be observed with S-waves of about 1500 m/s and 2500 m/s, respectively.

Geology
From a geological point of view, the seismic stations used in our study are located north of the Alpine Chain, the collisional orogen resulting from the convergence of the European and African plates and the closure of the Alpine Tethys (e.g. Frisch, 1979;Haas et al., 1995;Stampfli et al., 2001;Tricart, 1984;Trümpy, 1960). As a result of a complex deformation process, several tectonic units with different deformation styles formed. The units of interest in this work are represented, from south to north, by the Helvetic Zone, the Molasse Basin and the Jura Mountains. The Helvetic Zone extends from France to Austria and overlies, towards north, the Oligo-Miocene clastic sequence of the Molasse Basin (Pfiffner, 1993). Internally, the Helvetic Zone is divided in two sub-areas by major thrusts: the Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings Infrahelvetic Complex (below) and the Helvetic Nappes (above). The geological sequence of the Helvetic Nappes consists of limestones, marls and shales of Mesozoic and Triassic age. These rocks overlie the southern edge of the European continental platform with a tectonic contact (Infrahelvetic Complex). The units below consist of a basement of crystalline-metamorphic rocks and sediments of Paleozoic age (Ramsay, 1981) with sediments of Mesozoic and Triassic age on top. The basement, strongly deformed during the Alpine evolution, presents elongated domes. Towards south, the Helvetic Zone is overlain by the Penninic and Austroalpine Nappes representing the European continental margin and the Adriatic margin of Tethys, respectively. North of the Helvetic Zone, the Molasse Basin and the Jura Mountains represent the foreland and the fold-and-thrust belt of the Alpine Chain, respectively. The Swiss Molasse Basin is the northern foreland basin extending from Lake Constance in the north-east to Lake Geneva in the south-west. The thickness of the sediments increases from about 800 m in the north to more than 5 km in the south. The geological sequence of the Swiss Molasse Basin consists of a Precambrian to Paleozoic crystalline basement with Mesozoic and Tertiary units on top.
The crystalline basement is dissected by NW-SE faults generating graben structures parallel to the Alpine Chain (McCann et al., 2006;Naef & Madritisch, 2014). These grabens were filled during the Late Carboniferous and the Permian ages with clastic sediments of fluviolacustrine and alluvial fan sedimentation (Diebold, 1989;Matter, 1987). As shown in Sommaruga et al. (2012), several deep boreholes drilled Permian and Carboniferous sediments (e.g. borehole Weiach 1), while others drilled directly the Carboniferous sediments and no Permian sediments (e.g. borehole Entlebuch 1). Above the crystalline basement, the Mesozoic units with thickness up to 2 km are represented by Lower Triassic to Upper Cretaceous shallow marine sediments (Mock & Herwegh, 2017). The Mesozoic sequence, completely represented in the western sectors, lacks the Cretaceous sediments in eastern Switzerland (Sommaruga et al., 2012).  (Berger et al., 2005;Homewood et al., 1986;Schlunegger et al., 2007;Sinclair & Allen, 1992;Sinclair et al., 1991). North of the Molasse Basin, the Jura fold-and-thrust belt forms the frontal portion of the western Alps. The geological sequence consists of a basement with granite, gneiss and Paleozoic sedimentary rocks above which continental deposits, limestones and evaporites of Triassic age were deposited (Pierce, 1966). The evaporites, still present in the upper Triassic, alternate with shales, reaching a total thickness of about 1 km. The Jurassic and Cretaceous rocks are limestones, shales, and marls. Above these, Eocene to Pliocene strata were deposited. Based on the tectonic style, the Jura is divided into an external portion towards north called Plateaux and in an internal one, characterized by major folds, referred to as Haute Chaîne. The 23 analyzed sites are located in three different tectonic units: the seismic stations SARK and SBUH are located in the Helvetic Nappe of the Helvetic Zone, SBEG is in the eastern sectors of the Jura, while the remaining stations are all distributed in the Molasse Basin (Fig. 1). Station WALHA is located in Germany close to the border with Switzerland and the geology corresponds to sites in the Swiss Molasse Basin.

H/V Calculation
The horizontal-to-vertical spectral ratio curves inverted in this work are calculated for the earthquake recordings from the arrival of S-waves. In order to pick the shear wave, time windows of four minutes in length were selected from the continuous signal, starting one minute before the P-wave arrivals picked by the routine location team of the Swiss Seismological Service as reference. To determine the quality of the signal, the signal-to-noise ratio was extracted between 5 and 1 min before the P-wave arrival of each event. To identify and pick the S-wave arrivals, the data were filtered between 0.2 and 30 Hz using a 4 th order causal Butterworth bandpass filter. The picking was manually performed for each event without attributing an uncertainty or quality value to the picked value. Time windows of 40 s starting at the arrival of the S-waves were used for the computation of the H/V curves. Several tests were performed to choose the appropriate length of the time window and include both the S-wave and the coda; as already shown by Kawase et al. (2011). The S-waves and the coda give similar H/V curves and are both of interest in this work. For this reason, and in order to reach the condition of a diffuse wavefield, a single time window of 40 s was selected starting at the arrival of the S-wave. The seismic noise window of four minutes length, extracted before the P-wave arrival, was divided into 6 sub-windows of 40 s each to have windows of equal length and equal Fourier spectrum frequency range with the noise before the earthquakes. The earthquake signal and the noise windows of 40 s length were tapered at both extremities using a 5% Tukey window, zero-padded Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings at the end, and converted into the frequency domain by Fast Fourier Transform (FFT). The six noise spectra were then averaged with the geometric mean method to remove spurious signals. The final earthquake and noise spectra were then smoothed using a Konno-Ohmachi window with a b value of 60 (Konno & Ohmachi, 1998). The signal (earthquake) to noise (ambient seismic noise) ratio was calculated for each frequency and component. Only for frequencies with a signal-to-noise ratio higher than 3 for all components simultaneously, the signal was used for the horizontal-to-vertical spectral ratio computation. The H/V ratio was calculated by dividing the squared average of the two horizontal components by the absolute value of the spectrum for the vertical component (Nakamura, 1989), as shown in the formula: where |N(f)|, |E(f)| and |Z(f)| are the absolute spectra of the north, east and vertical components, respectively. The same procedure was repeated for all recorded earthquakes. Once the H/V curves were computed for each selected event, an average H/V curve was calculated using the geometric mean formulation. The mean was calculated for frequencies where the signal-to-noise ratio was higher than 3 for all three components simultaneously using a minimum of 5 events. The frequency range of each computed H/V curve depends on the sensor type, the number of the picked events, the signal-to-noise ratio and the magnitude and distance of the recorded earthquakes. The values for all stations are given in Table 1.

Inversion
The theoretical expression of the H/V ratio for the S-waves in the diffuse-field assumption (Sánchez-Sesma et al., 2011) is: where Im G Eq horizontal ð0; 0; xÞ Â Ã and Im G Eq vertical ð0; 0; xÞ Â Ã are the imaginary parts of the Green's function for the horizontal and the vertical seismic motion, respectively. This formulation works for a specific circular frequency (x ¼ 2pf ) when the source and the receiver are co-located. Following Kawase et al. (2011), Eq.
(2) can be transformed into: here, a H and b H are the P-and S-wave velocities of the seismological bedrock, respectively, and TF 1 x ð Þ and TF 3 x ð Þ are the transfer functions for the horizontal motion due to the S-wave and for the vertical motion due to the P-wave, respectively. In Eq. (3), the incident angles of single P-and S-waves are not taken into account because the average of multiple events with different azimuths and incident angles generates a wavefield with vertically incident waves. Readers are referred to Kawase et al. (2011) for further explanations and numerical verifications.
To reconstruct the velocity profiles below the investigated site, the averaged H/V curves were inverted with the inversion algorithm proposed by Nagashima et al. (2014) and based on the HHS method, a combination of a genetic algorithm, and simulated annealing. The genetic algorithm starts by generating a population of individuals with a specific set of properties. After each iteration, a set of individuals (parent models) is stochastically selected from the initial population by a fitness function, which defines how close a given solution is to the optimum solution for the desired problem. The properties of the selected individuals are later combined with each other to create a new population of individuals (child models) using the crossover and mutation operators. In the HHS technique, the choice of new models is similar to that in the annealing technique, where the probability to create new child models is linked to the increase of the generation number (k) by where T is the temperature and k increases of one unit every ten iterations.
To evaluate the fit between the results of the inversion (d inv ) and the observed H/V curve (d obs ), the HHS algorithm computes the misfit (m) as: where f min and f max define the frequency range of the observed H/V curve. After each iteration, the misfit of the child model is compared to the misfit computed for the parent model and the difference between the misfits of the two generations (D m ) is calculated. If D m is zero or negative, the child model passes to the next generation. Otherwise, if the value is positive, the child is carried over to the next generation with a probability of e ÀD m =T k ð Þ : If D m is positive but the child model is not carried to the next generation with a probability of 1-e ÀD m =T k ð Þ , the parent model moves to the next generation. The HHS inversion is a mix of the simulated annealing technique and genetic algorithm routines. During the inversion, 200 iterations are performed using the initial population of 200 randomly generated models within the parameter space defined by the user. The initial temperature (T) is an operator of the simulated annealing technique defining the starting point of the inversion; its initial value is set to 100 and it decreases after each iteration following Eq. (4). Mutation and crossover are two operators in the genetic algorithm, which increase the diversity in a population and introduce recombination between the individuals of the same population, respectively. These operators are set to 0.1 and 0.7, respectively, and act from the second generation onwards.
The inversion algorithm explores the parameter space and looks for the model with the lowest misfit, able to fit the average H/V curve given as input. The initial profile consists in layers where the initial parameters (seismic velocities, thickness, Poisson's ratio and density) are defined according to the available information. A range of variability is then defined to allow each parameter to change and adjust according to the input curve: i.e. a variability of ± 10% for an initial S-wave velocity of 1000 m/ s will corresponds to a range for the shear-wave velocity that varies from 900 to 1100 m/s. The variability values tested for each inversion, same input H/V curve and parametrization, are ± 10%, ± 20% and ± 50%. A small value (i.e. 10%) will give an inverted velocity profile which will be close to the initial profile, while a value of 20 or 50% will allow the algorithm to look for parameters which are quite different from the initial values.
In order to show the performances of the HHS inversion algorithm and the importance of defining the parameters for the initial profile, we inverted the average H/V curve from earthquakes for the seismic station SBERN with two different initial models. In the first inversion (Fig. 2), the initial profile consists of 20 layers with an equal thickness of 100 m. The Pand S-wave velocities and densities, instead, are set with linear gradients. To define the half-space seismic and a unique mandatory constraint for the HHS inversion algorithm, Nagashima et al. (2014) defined the velocity structure and the bedrock velocity based on the studies performed by the National research Institute for Earth science and Disaster prevention (NIED) and of Kawase and Matsuo (2004), respectively. This type of information is not available for the Swiss stations and we decided to use literature data. Campus and Fäh (1997) estimated the velocity profiles for several areas of Switzerland among which the Molasse Basin (Mittelland) is of interest in this work. Taking into account the values found by these authors for the P-and S-waves and the S-wave results of Chieppa et al., (2020aChieppa et al., ( , 2020b for two sites in the Swiss Molasse Basin, we defined minimum values of 100 m/s for the P-wave velocity and of 50 m/s for the S-wave velocity at the surface and fixed values of 6200 and 3500 m/s at the half-space, respectively. Their initial values increase with depth, following a regular step. The minimum value for the Poisson's ratio allowed by the inversion algorithm is 0.25. During the inversion, the seismic velocities and the thicknesses of the layers are allowed to explore the entire parameter space without any constraints, in order to find the model best fitting the H/V curve. The values for the density are calculated taking into account the estimated S-wave velocity Vs (Nagashima et al., 2017): Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings 765 The velocity model for the second inversion consists of shallow and deep layers. The number of shallow layers and the initial parameters of the initial profile are site-dependent and are linked to the velocity profiles estimated by the site characterization measurements performed nearby (Table 1). The number of deep layers is fixed to 10; these are connected to the shallow layers and explore the deep subsurface using the initial values set by the user. Based on the information contained in the H/V curve, the values will change and adjust without restrictions during the inversion (Fig. 3). The thicknesses and the 766 D. Chieppa et al. Pure Appl. Geophys. shear wave velocities of the shallow layers are allowed to adjust during the inversion within userdefined ranges, while the same parameters for the deep layers are defined investigating the entire parameter space without restrictions. The density and the P-wave velocity are linked to the estimated S-wave velocity by Eq. (6) and the formulation reported by Nagashima and Kawase (2021), respectively. The P-and S-wave velocities are in general forced to increase with depth. At certain sites, in agreement with the information provided by the site characterization analysis, a low-velocity zone (LVZ) is allowed in the shallow layers. This is kept by the HHS algorithm until the fit of the earthquake H/V curve is better than using a similar velocity model with a gradient: in that case the LVZ disappears. The P-and S-wave half-space velocities are set as for the first inversion to 6200 and 3500 m/s. In the inversion, the maximum investigated depth is not known in advance and is directly linked to the information contained in the H/V curve given as input. For this reason, the inversions performed in this work, with the exception of the examples shown for station SBERN (Figs. 2, 3), are repeated several times changing the shear wave velocity of the halfspace to 2500, 3000, 3500 and 4200 m/s. In a similar way, the P-wave velocities of the half-space were adjusted, while the initial thickness and the density did not change. The seismic velocities were defined taking into account the distribution of the velocity profiles estimated for different tectonic regions in Switzerland by Campus and Fäh (1997), as well as the results obtained by Chieppa et al., (2020aChieppa et al., ( , 2020b at two sites in the Swiss Molasse Basin. After assuming that 2500 and 4200 m/s are probably two extreme values for the 23 analyzed sites, we discarded these values and kept only the results where the S-waves for the half-space is equal to 3000 or 3500 m/s. The inversions were performed ten times for the same H/V curve using the same initial model and the same ranges of variability for the shallow and deep layers. In each inversion a different number of seeds is created to check the algorithm performances in terms of stability and similarity; the resulting models are then plotted together to show the standard deviation ranges around the average H/V curve. While the deep layers were always allowed to change without constraints, three different variability values were used for the shallow layers: 10%, 20% and 50%. The higher the variability value chosen, the higher the size of the parameter space that will be investigated, meaning that the computation time will be longer and the final velocity profile might be quite different from the initial profile.
The difference between the observed H/V curve is represented in terms of h H=V using the formula where N is the total number of frequencies in the input curve, d obs i is the input H/V curve, d inv i the inverted H/V curve and r i is the standard deviation of the logarithmic H/V curves for each seismic station. The h H=V defines the color of the curves in all subplots (Figs. 2,3,4,5,6,7,8,9). The color bar in each figure ranges between 0 (dark brown) and 0.3 (dark yellow) according to the h H=V value calculated for each inversion.
In a similar way, we compared the measured Rayleigh-and/or Love-wave dispersion curves from site characterization analysis, whenever available, with the synthetic dispersion curves computed for the estimated velocity models from each inversion. While it might sound important and reasonable to compare the measured dispersion curves only with the synthetic results for the HHS inversion using a b Figure 2 Inversion results for the strong-motion station SBERN, performed using 20 completely unconstrained layers. a H/V curves, b S-wave velocity profiles: zoom to the shallow structure around the site characterization velocity profiles (left) and entire profiles down to the half-space (right), c Rayleigh-wave dispersion curves, d Lovewave dispersion curves, e empirical amplification function and SHwave transfer functions, f map showing the epicenters of the recorded earthquakes. The input H/V curve and the measured amplification function curve are represented in black, the H/V curve and the estimated velocity profile from site characterization measurement and the measured Rayleigh-and Love-wave dispersion curves are in dark green, and the initial profiles in grey. The Swave half-space velocity is 3500 m/s. The color of the inverted results is referred to the h H=V calculated for each inverted H/V curve and ranges from dark brown (low h H=V ) to yellow (high h H=V ).
Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings velocity profile with unconstrained layers, we will show that the same comparison is needed also when constraints are set since the initial profile might be improved. In contrast to the H/V curves, the dispersion curves were not used in the inversion and are here used only for comparison purposes as they had been used to invert for the velocity profiles during the site characterization analysis. Also in this case, in order to objectively evaluate the quality of our results, we computed the difference between the measured dispersion curves and the theoretic curves for the inverted model taking into account the dispersion curve type (e.g. Rayleigh or Love) and the corresponding mode attribution, e.g. the measured fundamental mode for Love waves with the theoretical fundamental mode for Love waves. Due to the different number of measured dispersion curves and the different number of points (frequencies) in each curve, we calculated a total h disp value taking into account the number of frequencies in each curve as shown in the following formula: where N R and N L are the total number of frequencies of the Rayleigh-and Love-wave dispersion curves of all modes, respectively, d obs i is the measured dispersion curve and d inv i the synthetic dispersion curve computed for the velocity profile estimated by the Geopsy package (Wathelet et al., 2020). To illustrate the attribution to a class rather than another, the h disp is used as described in the section Inversion results for all presented sites.
Besides the comparison of the H/V and dispersion curves, the 1D SH-transfer functions and the empirical amplification function (ESM) curves are shown (Figs. 2,3,4,5,6,7,8,9, plot e). The 1D SH-transfer function is calculated for each estimated velocity profile for vertically propagating SH-waves (e.g. Roesset et al. 1970) taking into account the ratio between the amplitude of the SH-waves at the surface and the amplitude of the same wave at the bedrock interface. The obtained transfer function will depend on the S-wave velocity of the half-space, which is fixed during the entire inversion and different at different sites. To compare the obtained results with the ESM curve both curves are corrected for the Swiss reference rock profile (Poggi et al., 2011), as described in Edwards et al. (2013). The empirical amplification function is computed for the elastic and the anelastic cases of the recorded earthquakes with a signal-to-noise ratio higher than 3. Once the amplification spectrum is calculated for each event individually, the averaged ESM curve is computed (Edwards et al., 2013). For the calculation of the anelastic ESM, a site-dependent attenuation operator (j) is used; it depends on the epicentral distance and the path of the seismic rays during their propagation. The value of j is continuously updated with new recordings.
While the ESM measures the propagation of all seismic waves, the SH-wave transfer function considers the contributions of SH-waves only. Based on this information, we do not expect one-to-one correspondence as the source composition is different, but we are interested in similarities between these two curves in terms of shape and/or amplitude. As for the dispersion curves, the color used for the SH-wave transfer functions is related to the h H=V calculated using Eq. (7).

Station SBERN and the Benefit of Constraining the Shallow Structure
Two inversions were performed for station SBERN using two initial profiles. The first velocity profile consists of 20 unconstrained layers with a thickness of 100 m each and seismic velocities and densities increasing with depth (Fig. 2). The second model is a mix of twelve shallow layers partially constrained to the a-priori information and ten deep and completely unconstrained layers (Fig. 3). The initial parameters of the shallow layers are attributed using the velocity profiles from the previous ambient noise vibration measurements. 10% variability was allowed for the initial parameters of the shallow layers to adjust and fit the earthquake H/V curve. Below the shallow layers, ten layers were added to explore the deep structure. These 10 layers investigate the deeper portion of the subsurface and adjust the maximum investigated depth with respect to the Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings 769 initial model. The seismic velocities and densities of these layers follow a gradient, while the thickness is fixed to 100 m. The P-and S-wave velocities of the half-space are fixed to 6200 and 3500 m/s, respectively, and represent the only constraint requested by the HHS algorithm. Both tests show agreement between the H/V curves generated by the inversion and the measured curve over the entire frequency range (Figs. 2a, 3a). The results in Fig. 2a show a wide distribution of the inverted H/V curves around the measured one; in Fig. 3a, the inverted H/V curves are closer to the measured curve and less scattered. With the halfspace velocities as the only constraint (Fig. 2), the S-wave velocity profiles show a gentle gradient down to about 800 m with larger differences at greater depths. The velocity profiles obtained when constraints are set in the shallow subsurface show a steep gradient down to 200 m (Fig. 3), almost constant velocity between 200 and 800 m, and several impedance contrasts below (e.g. at about 800 m, 1200 m, and 1900-2100 m).
The wide distribution of the inverted H/V curves seen in Fig. 2 is confirmed by the slightly higher misfit values for the H/V curves (h H=V : 0.13-0.22) when compared with the misfit values in Fig. 3 ranging from 0.08 to 0.10. For each test, the measured Rayleigh-and Love-wave dispersion curves (green color) are shown together with the ones derived from the velocity profiles of the H/V inversion. These were not directly used to constrain the HHS inversions but, for comparison purposes, to discriminate which of the obtained results fits the input H/V curve and the available data (e.g. dispersion curves) at the same time. The measured dispersion curves in Fig. 2c, d, with h disp between 0.52 and 0.85, do not fit the synthetic ones. The fit of the Rayleigh-and Love-wave dispersion curves in Fig. 3c, d is much better (h disp : 0.05-0.08); small and negligible disagreements can be observed for the Rayleigh-wave fundamental mode below 3 Hz and between 5 and 8 Hz. In addition to the dispersion curves, we plotted the SH-wave transfer function for each inverted velocity profile and the observed empirical amplification function for station SBERN (Figs. 2e, 3e). The SH-transfer functions from the inversion (colored curves) are similar in terms of shape with the empirical amplification function (black solid line) for the seismic station SBERN, but present amplification values which are higher than the ESM upper standard deviation curve (black dashed line). In Fig. 3e, the SH-amplification curves overlap over the entire frequency range. Up to 0.8 Hz the SH-transfer curves and the ESM curve have similar amplification values with a gentle drift starting at higher frequencies up to 20 Hz. Between 0.5 and 9 Hz, the SH-amplification functions are still located within the boundaries defined by the standard deviation of the ESM function and fit perfectly the ESM curve between 1 and 3 Hz. At higher frequencies, similar peaks with higher amplification values can be observed up to 10 Hz.
The good agreement found for the H/V and dispersion curves (Fig. 3a, c, d) confirms the importance to constrain the layers of the initial model for two reasons: to drive the inversion to a more realistic solution thanks to the a-priori information and to reduce the number of unknowns in the inversion process -and consequently, the computation timeleaving only the 10 deeper layers completely unconstrained. The disagreement found at high frequency between the SH-wave transfer functions and the empirical amplification function might be explained by taking into account the shorter frequency range of the inverted H/V curves and/or the measured dispersion curves, the Swiss rock reference used for the correction (Poggi et al., 2011), or the construction site where the station is installed.
The example in Fig. 3 shows a way to combine the site characterization results with the earthquake recording information and the possibility to extend the velocity profiles to larger depths. Based on the findings in Fig. 3, the velocity profiles used from now on will consist in a variable number of partly constrained layers, linked to the site characterization measurement, and ten unconstrained layers at depth. For the following examples, the S-wave velocities of the half-space were set to 3000 m/s or 3500 m/s, while the variability allowed to the shallow layers was set to 10%, 20% or 50%. The quality of the velocity profile estimated by the site characterization analysis is not known in advance as well as the maximum depth investigated. For this reason, we tested increasing values of variability and half-space 770 D. Chieppa et al. Pure Appl. Geophys. shear-wave velocity and analyzed the obtained results. When the estimated velocity profile from site characterization is good, a small value (10%) is enough and no improvement is needed for the shallow layers; on the contrary when the velocity profile has no or low resolution in the shallow layers (e.g. few and thick layers), the shallowest layers need more space to adjust. For station SBERN, an additional example is shown in Fig. A2 (Electronic Supplementary Information) setting variabilities for the shallow layers equal to 10% and 50% and a half-space velocity for the S-wave equal to 3000 m/s. The inverted velocity profiles with increasing S-wave velocities for the half-space overlap and follow a similar gradient down to the bedrock ( Fig. A3-Electronic Supplementary Information).

Inversion Results
The results shown in this section summarize an intense analysis performed using the HHS algorithm at 23 sites. For each analyzed site, we defined an initial model consisting of shallow layers, with parameters linked to the results of the site characterization measurements, and deep layers, with initial values free to change during the inversion. Starting from the same initial profile with a shear-wave velocity for the half-space equal to 3500 m/s, several independent inversions were performed, changing the parameter variability attributed to the shallowest layers (± 10%, ± 20% and ± 50%). The entire analysis was then repeated with a half-space velocity to 3000 m/s. Not all tested combinations produced noteworthy results and therefore, only one of the many configurations is displayed for each investigated station in this work. The results for the 23 sites were divided in three classes taking into account some recurring features (e.g. simultaneous fit of H/V and measured dispersion curves) and not the computed h H=V values, as the result ranges of the h H=V are pretty narrow for all tested inversions and do not show significant variations. The three classes have the following features: • Class 1 groups the sites where the results of the site characterization are good and the estimated velocity profiles explain all picked curves in a satisfactory way. The parameters of these velocity profiles require small variability to adjust during the HHS inversion and a good fit of the H/V curve and the measured dispersion curves is achieved. • Class 2 contains the stations whose inverted velocity profiles share similar characteristics irrespective of the allowed variability in the shallow layers, providing relatively stable results in the deeper structures. • Class 3 includes the remaining sites where the inverted H/V curves do not fit the average H/V curve from earthquakes, the synthetic dispersion curves from the inverted velocity profiles are not in agreement with the measured dispersion curves, and/or the SH-wave transfer functions are off when compared with the empirical amplification function. This fit is evaluated taking into account all the above-mentioned curves at the same time. This class includes sites where the site characterization measurements were generally difficult to perform (e.g. rock sites, low quality site characterization analysis).
In Table 2, we give an overview of the group assignment for the different seismic stations together with the parameters chosen, including the maximum S-wave velocity in the half-space and the allowed variabilities shown for each site here in the main document and in the Electronic Supplementary information. We were looking for solutions in Class 1 or 2 explaining all the observations. We started with a half-space velocity of 3500 m/s and in case of a unsatisfactory or no solution, the S-wave velocity for the half-space was decreased to 3000 m/s. In case even this value did not provide a good fit, the station was attributed to Class 3. In the following, we will present an example for each class; the results for the 20 remaining sites can be found in the Electronic Supplementary information. The colors used to plot the H/V curves, dispersion curves and SH-wave transfer functions are the same for the curves of the same model in all plots and are based on the h H=V value calculated between the H/V curve from earthquakes and the inverted H/V curves as defined in Eq. (7). Based on the good performances of the HHS inversion algorithm and the low values obtained for Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings 771 the analyzed sites, the h H=V misfit ranges are displayed between 0 and 0.3.

Class 1: Good Fit of All Observations with Small Variabilities (Example: SYVP)
The average H/V curve computed for site SYVP was inverted between 0.19 and 28.7 Hz. The initial profile consisted of 26 layers over the half-space: the 16 shallowest layers were constrained by the parameters estimated by the site characterization (Michel et al., 2014), while the remaining 10 layers were set following a linear gradient. The results shown in Figs. 4, 5 were obtained using an S-wave half-space velocity of 3500 m/s and a variability value for the parameters of the 16 shallow layers of 10% in Fig. 4 and 50% in Fig. 5. In both inversions, the inverted H/ V curves fit the input data over the entire frequency range. The peak (f 0 ) at around 1 Hz is partially represented in Fig. 4a with amplitudes around 10; the fit becomes even better in Fig. 5a in terms of amplitude and shape. Also the trough, above the f 0 peak frequency, is well represented in both cases. The Love-wave dispersion curves calculated for the inverted velocity profiles in Fig. 4d show perfect agreement with the measured curve for the fundamental mode and the first higher mode. Disagreement exists for the second higher mode (4.3-5.6 Hz), where the inverted dispersion curves are shifted towards higher velocities. The Rayleigh-wave dispersion curves in Fig. 4c show perfect agreement for all three modes. Increasing the variability of the shallow layers to 50% (Fig. 5c, d), a systematic shift towards higher frequencies is seen for the synthetic Rayleigh-and Love-wave dispersion curves from the inverted velocity models when compared to the measured dispersion curves. The similar fit of the H/V curves in Figs. 4a and 5a is confirmed by the h H=V values computed with Eq. (7): 0.20-0.21 and 0.12-0.16, respectively. The difference calculated between the inverted and measured dispersion curves ranges from 0.15 to 0.35 in Fig. 4c, d and between 0.43 and 0.72 in Fig. 5c, d, confirming the better fit when a low variability value is used. The depth of the half-space for the velocity profiles in Fig. 4b spreads between 600 and 1700 m with a big velocity contrast between 500 and 700 m. The S-wave velocity profiles in Fig. 5b (right column) reach the values of the half-space at slightly shallower depths. The earthquake empirical amplification function measured at SYVP, shows a curve with several sharp peaks. Similar peaks can be observed for the SHwave transfer functions of the inverted velocity profiles in Figs. 4e and 5e. The main peak at low frequency in Fig. 4e is comparable in terms of amplitude and shape to the peak of the inverted curves, but it is shifted towards higher frequencies (around 1.2 Hz). In Fig. 5e, the observed peak is centered at the same frequency but shifted towards higher amplitudes. Several peaks with high amplitudes are also seen at high frequencies in the SHwave amplification functions when small variability values are set (Fig. 4e); the same peaks have smaller amplitudes in Fig. 5e and fit slightly better the ESM curve. The initial model from site characterization explains the H/V curve calculated from earthquakes Table 2 Classification of the analyzed seismic stations in three classes. The values with the percentage in each cell refer to the allowed variabilities for the shallow layers shown in the main document or in the Electronic supplementary information; the velocity corresponds to the maximum S-wave velocity for the half-space. The example shown in Figs. 2 and 3 for station SBERN was obtained using variability values of 10 and 100% with an S-wave half-space velocity of 3500 m/s and it is not included in this table (h H=V is 0.50). Small variability values, e.g. 10%, generate a good fit of the H/V and dispersion curves at the same time; higher values (50%) slightly improves the fit of the H/V curve, while the fit of the dispersion curves gets worse. The velocity structure in the shallow subsurface is well constrained  by the site characterization measurements and we can effectively determine the velocity structures in the deeper part by using the H/V curves.

Class 2: Good Fit of all Observations Irrespective of the Variabilities (Example: SARK)
At SARK, we inverted the averaged H/V curve computed for 84 earthquakes. Figures 6, 7 show the results obtained for a S-wave half-space velocity of 3500 m/s and a shallow layer variability of 10% (Fig. 6) and 50% (Fig. 7). A good fit with the input H/V curve is obtained for both tests over the entire frequency range. The picked fundamental and the first higher modes of Rayleigh-and Love-wave dispersion curves fit the generated dispersion curves. While the inverted dispersion curves overlap quite well over the entire frequency range in Fig. 6c, d, the dispersion curves in Fig. 7c, d present similar shapes but are scattered around the measured curves. The h H=V values of the H/V curves are similar in both inversions and range between 0.07 and 0.09 in Fig. 6 and between 0.08 and 0.10 in Fig. 7. Similar results can be seen for the h disp : 0.27-0.42 in Fig. 6 and 0.32-0.64 in Fig. 7. The inverted S-wave velocity profiles in Fig. 6b and Fig. 7b show a steep gradient starting at about 200 m, reaching high velocities similar to the bedrock velocity at a depth of around 400 m. The SH-amplification functions computed by the inverted velocity profiles of both tests agree with the ESM curve (black solid curve) in terms of shape, but not amplification. This shift, not visible at other sites (see Electronic Supplementary information), might be explained in two ways: taking into account the tectonic evolution of the area or in terms of choice for the Swiss rock reference. From a tectonic point of view, the seismic station SARK is located in the Helvetic Nappes, overthrusting the Swiss Molasse Basin towards north. The depth of the thrust dividing the two units is located at around 4-4.5 km, based on the geological interpretation in Sommaruga et al. (2012). The thrust, if interpreted as the bedrock halfspace, might influence the calculation of the empirical amplification function since part of the energy generated by the earthquakes might be reflected downward and, therefore, not reach the surface. Taking into account the chosen reference rock (Poggi et al., 2011) instead, we might say that the rocks beneath the SARK seismic station are more rigid than the rocks at the sites used for the definition of the Swiss reference rock, where the S-wave half-space velocity is set to 3200 m/s. The site characterization results provided P-and S-waves velocity profiles with a good resolution so that the inverted velocity profiles do not move away from the starting model, irrespective of the chosen variability. These are the sites where the velocity structures both in the shallower and deeper parts are well constrained by the H/V curve alone. Moreover, the stations of this class show that the HHS algorithm can be a useful technique to investigate deeper structures and reach depths that are not resolved by a classical site characterization measurement. Accessing these depths would require the deployment of arrays with inter-station distances of several kilometers, which is not everywhere feasible, for example, in valleys or cities where the space is limited.

Class 3: No Simultaneous Good Fit of All Observations Possible (Example: STIEG)
The results shown in Figs. 8, 9 were obtained inverting the averaged H/V curve computed for 122 seismic events in the frequency range 0.25-30.11 Hz. The initial profile consists of 8 shallow layers, identified using passive ambient vibrations, and 10 deeper unconstrained layers. The S-wave velocity of the half-space was set to 3500 m/s. The first test was performed with a variability of 10% for the shallow layers (Fig. 8). The inverted H/V curves are in general agreement with the input data between 0.43 and 1.85 Hz; at higher frequencies, the inverted H/V curves present shapes which are different from the measured H/V curve and from the initial profile. While the fit of the measured and inverted H/V curves is achieved in a narrow frequency range (0.5-1 Hz), the fundamental modes of Rayleigh and Love waves match the inverted dispersion curves over the entire frequency range. The SH-amplification function curves have shapes which are comparable with the empirical amplification function and values which are within the standard deviation curves (Fig. 8e). Below 2.5 Hz, we see a rapid decrease of the amplification values towards 1. The second inversion Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings  was performed using a variability for the shallow layers equal to 50%. The results in Fig. 9 show a very good fit of the H/V curves over the entire frequency range between 0.4 and 20 Hz. The dispersion curves of the inverted velocity profiles do not fit the measured dispersion curves but are shifted towards higher velocities. The SH-amplification functions reported in Fig. 9e have the same trend as in Fig. 8e for frequencies below 2.5 Hz. Above 2.5 Hz, the amplification curves in Fig. 9e fit the ESM curve reproducing the wide peak between 10 and 16 Hz and the troughs. The h H=V values range between 0.19 and 0.20 in Fig. 8a and 0.08 and 0.11 in Fig. 9a; the h disp changes from 0.11 to 0.25 in Fig. 8c, d to 0.32-0.98 in Fig. 9c, d. The velocity profiles of both examples do not move away from the input profile (gray curve) and investigate the subsurface without any strong impedance contrasts, with a moderate gradient down to 700 m. The 50% variability improves the fit of the measured H/V curve and of the empirical amplification function, but reduces the fit of the measured dispersion curves. Vice versa, the results in Fig. 8a are in disagreement with the measured H/V curve. All the stations of this class have Vs 30 ranging from 422 to 787 m/s corresponding to classes B and C in the Swiss building code SIA261. The reasons why the inverted curves do not fit the measured curves are numerous: the array deployment was not performed exactly at the location of the permanent installation, lateral variations occur underneath the study area, sites with more compact rocks, or the uncertainty of the analysis/interpretation of the measured curves that result in a slightly biased estimation of the velocity profiles. In both tests, the inverted H/V curves are flat at low frequency and do not fit the narrow peak at about 0.3 Hz in the observed H/V curve. This could be linked to the ease of the chosen velocity profile that does not allow a proper investigation of the deeper structures.

Discussion
The analysis performed for station SBERN showed the importance of a-priori data (e.g. geophysical information) to constrain the parameter space. Two velocity profiles were used for the H/V inversion: the first consisting in 20 completely unconstrained layers and the second with 12 shallow layers constrained to the available information (e.g. previous site characterization) and 10 deep unconstrained layers. Using both initial models, a good fit was found between the H/V curve from earthquakes and the inverted velocity profiles. In order to discriminate which of the inverted velocity profile explains the subsurface better, the measured dispersion curves from site characterization were used for comparison. For most of the runs of Fig. 2, i.e. without constraints by the site characterization measurements, we noticed disagreement between the synthetic (generated by the estimated velocity profiles) and the measured dispersion curves. A good agreement was instead found (shown in Fig. 3), where the superficial layers were constrained to the velocity profiles estimated by the site characterization measurements. The example of station SBERN shows the importance of constraints in the inversion and the possibility to restrict the parametrization space to drive the inversion towards the solution with the lowest misfit corresponding to the most reasonable minimum of the inversion. In contrast to our expectations, the h H=V values computed for Eq. (7) are smaller when constraints are set (Fig. 3) than when the layers are completely unconstrained (Fig. 2). The higher values of Fig. 2 are linked to the larger accessible parameter space of the unconstrained inversion and therefore to the slower conversion towards the solution. This shows a twofold advantage of the constrained inversions: (1) the inversion converges faster to the solution and (2) the improvement of the subsurface velocity profiles without moving away from the initial profile. Another feature we recognize in Figs. 2, 3 is the distribution of the velocity profiles in relation with the distribution of the inverted H/V curves. In Fig. 2b, the scatter of the inverted H/V curves reflects in S-wave profiles spread over a wide range of depths (800-2800 m); in Fig. 3b, the better fit of the inverted H/V curves corresponds in more similar velocity profiles with comparable impedance contrasts and maximum investigated depths.
Based on the main findings for station SBERN, further analyses of the other sites were performed using a velocity profile with a mix of shallow partly 780 D. Chieppa et al. Pure Appl. Geophys. constrained layers and deep unconstrained layers. The parameters of the shallow layers are constrained by the thickness, velocity and density information collected using site characterization measurements around the investigated site. A certain degree of variability is then allowed to the seismic velocities and thicknesses of the shallow layers. If the site characterization velocity profiles explain the measured dispersion/ellipticity curves computed for the site characterization analysis in a satisfactory way, it means that the subsurface properties at the permanent installation and in the area where the array measurement was performed are comparable. In this case, a variability value not exceeding 10% should be allowed to the shallow layers in the H/V inversion. Stable velocity profiles not moving from the initial model as in Class 2 were observed in about one third of the cases (8). A too high variability given to the parameters of the velocity profiles with good initial velocity profiles might generate a shift of the results towards less realistic solutions in another third of the cases (Class 1). In all performed tests, the layers added at the bottom of the velocity profile from site characterization are set without restrictions and expected to investigate the deeper structures by fitting the H/V curves. The ten deeper unconstrained layers added at the bottom of each velocity profile allowed an increase in investigated depth at all sites. While in most of the cases these layers added information at depths and increased the depth of the estimated bedrock, at sites like BOBI and SGEV (see Electronic Supplementary information), the half-space and the unconstrained layers were pushed upwards, close to the location of the half-space found in site characterization. At these two sites, probably due to the low sensibility of the H/V curves, the inverted velocity profiles were not able to resolve larger depth ranges, and located the half-space at 250 m in BOBI and 275 m in SGEV. Taking into account the depth of the half-space and the shape of the inverted velocity profiles, we identified P-and S-wave velocity profiles with three different gradients: steep (DAGMA, FLACH, SBEG, SBERN, SBUH, SEPFL, SOLB, STEIN and WALHA), moderate (HAMIK, SZUZ, BOBI, EMMET and ZUR), and low (SGEV, STGK, TORNY, WEIN, WILA and SOLZ). The steep gradient is obtained when the half-space is deep and the velocity profiles investigate deep portions of the subsurface; the low gradient instead presents a shallow half-space and high shear-wave velocities in the first 1-1.5 km. The combination between the type of gradient and the areal distribution of seismic stations does not show any specific relation which can be summarized here. This is probably due to the lateral variations of thickness for the Swiss Molasse Basin sediments and the far distance between the analyzed sites. We identified recurring features for the analyzed sites and classified them into three classes. The first class shows the most recurring trend and is represented by the seismic station SYVP; the a-priori information retrieved using ambient seismic noise measurements constrains the shallow layers of the initial velocity model. When a small variability value (10%) is set, the initial parameters of the shallow layers are allowed to change by a small percentage and adjust to fit the earthquake H/V curve. In Fig. 4, the allowed variability allows the fit of the H/V curve, dispersion curves, and amplification function at the same time. Increasing the variability to 50% (Fig. 5), the fit of the H/V curve and the empirical amplification function do not change, while the fit of the dispersion curves degrades. The inversion of the H/V curve using the HHS algorithm generates velocity profiles with lower investigation depth not able to explain the earthquake H/V curve and the measured dispersion curves at the same time. For the seismic stations of this class, we suggest using a variability not higher than 10% to allow the parameters of the shallowest layers to improve the resolution of the subsurface without moving away from the initial profile.
The features of the second class are shown for the strong motion station SARK. Similar to Class 1, the good quality of the a-priori information combined with the deep unconstrained layers and a small variability for the parameters of the shallow layers generates models that remain close to the initial parametrization. This occurs independently from the allowed layer variabilities. The shallow layers, as in Fig. 4, drive the inversion towards the presumably best result without moving away from the initial information; the deep layers below reconstruct the deep interfaces in the subsurface. For a site with a Vol. 180, (2023) Identification of Subsurface Structures Using H/V Curves from Earthquake Recordings good initial profile such as SARK, we can either allow a low or a high variability value and the inversion will not diverge from the initial model as confirmed by the h H=V values ranging between 0.08 and 0.10 independent of the chosen S-wave velocities of the half-space. These are the sites where the velocity structures both in the shallower and deeper parts are well constrained by the H/V curves alone as long as the initial structure lies within an adequate search range. The last group is presented by the results obtained for site STIEG. This example was selected to show the complexity in fitting the H/V curve from earthquakes, the measured dispersion curves from site characterization and the empirical amplification function at the same time. With a small variability value (e.g. 10%), the velocities and the thicknesses of the shallow layers do not move away from the initial parametrization and the inverted dispersion curves still fit the measured curves over the entire frequency range, while the inverted H/V curves show some disagreement with the measured H/V curve. Using a 50% variability, the fit of the calculated H/V curve and of the SH-amplification function improves, while the disagreement with the measured dispersion curves increases. This example represents a group of stations where site characterization was difficult due to the lack of sufficient space for the acquisition, or lateral variability in the structure, so that the site characterization might not completely be representative of the subsurface underneath the permanent seismic station and the calculated dispersion curves are not representative of the subsurface of the seismic station. Another possible explanation for the stations of this class is a possible directionality in the earthquake sources. Similar to the sites of Class 1, higher variability values (i.e. Figure 5) move the inverted velocity profile away without providing a simultaneous fit of the H/V curves, of the measured dispersion curves and of the amplification functions. Due to the difficulty in explaining the observations, the estimated velocity profiles for the stations of this class need to be treated with more uncertainty.
The a-priori information used in this work to set the initial parameters of the shallow layers was derived by the ambient seismic vibrations, MASW, or a mix of active and passive techniques, as reported in Table 1. In several cases, the inversion of dispersion and ellipticity curves generated P-and S-wave velocity profiles with strong velocity contrasts at the transition to the half-space; these were generally replaced by more gentle contrasts and velocity gradients after the H/V inversion.
We should also note that the H/V curves computed for the velocity models obtained for site characterization (green lines) tend to show disagreements in the high frequency range when compared with the observed H/V curves. This behavior is seen at almost all investigated sites at frequencies above 10 Hz and suggests a very shallow structure often not captured by the site characterization surveys and related to the instrumentation spacing or to the lateral variability of the structure close to the seismic station. The inverted H/V curves have a shape much closer to the observed earthquake H/V ratio (black line) suggesting that the inversion might improve the resolution of the very shallow layers. The details of the housing of the seismic station might also affect the H/V ratios at higher frequency which are not captured by site characterization or H/V inversion. For sites like SEPFL, SGEV in Class 1, HAMIK, STGK, WEIN in Class 2 or SOLZ in Class 3 of the Electronic Supplementary information, the theoretical H/V curves from site characterization (green curves) show shapes similar to the observed H/V curves, a lower peak amplitude at the frequency of the fundamental peak and a general shift towards higher frequencies. This is the direct consequence of the insufficient modeling of the velocity structures down to the seismological bedrock. The problem associated with too small S-wave velocity in the halfspace for the H/V modelling is discussed in Kawase et al. (2019). For site SBERN, we see an excellent fit of the measured H/V curve over the entire frequency range down to 0.33 Hz with a particular fit of the trough and the peak at 0.6 and 0.45 Hz, respectively. The good fit of these two elements allows a detailed reconstruction of the velocity profiles down to 2500-2600 m.
In order to confirm or reject the results of the H/V inversion, the inverted velocity profiles could be verified using other available information, such as measurements in deep boreholes. Unfortunately, most of the boreholes drilled in the Swiss Molasse Basin 782 D. Chieppa et al. Pure Appl. Geophys. have P-wave velocity profiles and only in rare cases the S-wave velocities were measured. For the H/V inversion performed in this work, the S-wave velocity profiles are much better constrained than the P-wave profiles and only these should be used for comparison. To conclude our analysis, we plotted the velocity profile with the lowest misfit for all shown cases with the results of passive site characterizations performed in the Swiss Molasse Basin (Chieppa et al., 2020a(Chieppa et al., , 2020b and with three velocity profiles (Mitteland, Innerschweiz, and Innerschweiz1) from Campus and Fäh (1997). The results in Fig. 10 show a general agreement in terms of seismic velocities and depths for all investigated sites.

Conclusions
In this work, we tested the Hybrid Heuristic Search (HHS) algorithm developed by Nagashima et al. (2014) and estimated for each investigated site, P-and S-wave velocity profiles. The analysis was performed inverting the average Horizontal-to-Vertical spectral ratio curves computed for the earthquake arrivals of local and regional events.
Several motivations drove this work. The first and most important one is the possibility to test the HHS technique to sites outside Japan (e.g. Switzerland), where the number of earthquakes and the local magnitude is lower. The remaining reasons are linked to the work of Nagashima et al. (2014) which inverted the H/V curves using earthquakes coming from a unique direction (east) and set the half-space velocities as the only constraint for the initial profile. In this work, the average H/V curves were computed for several seismic events distributed around the seismic station to support the diffusive field theory and show the different results obtained when additional constraints (e.g. near surface data) are set in the shallow layers of the initial profile. Finally, in order to test the stability of the HHS algorithm, each inversion was repeated ten times using the same calculated H/V curve and the same parametrization but a different number of seeds. The results of the ten HHS inversions, not providing an uncertainty value, were then plotted together for comparison to show the range of possible solutions.
The definition of the initial profile proved to be a key factor for the inversion. We used an initial profile that consists of a mix of shallow layers, with parameters defined by the available geophysical information, and deeper unconstrained layers investigating the interfaces in the subsurface down to the upper crust. Moreover, depending on the reliability of the a-priori data and on the goal of the analysis, a certain degree of variability should be allowed to the shallow layers. This value is intended to improve the information of the velocity profiles and fit the input H/V curve and the geophysical data. The allowed variability values are 10%, 20% and 50%. A small value should be selected when the a-priori information explains the recorded data in a satisfactory way and the HHS inverted velocity profiles should not move away from the velocity profiles estimated by other methodologies. A large value, instead, should be used when the disagreement between the inverted and the measured data (e.g. dispersion curves) is big and a different model should fit the available data. In all tests, the deep layers are free to move and expected to investigate the parameter space without restrictions. After identifying the initial profile for each site using the existing site characterization Figure 10 Shear-wave velocity profiles comparison. The models with the lowest misfit for all shown inversions discussed in this work and in the Electronic Supplementary information are shown in gray. The S-wave velocity profiles obtained by Chieppa et al., (2020aChieppa et al., ( , 2020b at two sites in Switzerland are shown in red. The results for Mitteland, Innerschweiz, and Innerschweiz1 areas from Campus and Fäh (1997)  results, we performed several inversions changing the seismic velocity of the half-space and the variability allowed to the shallow layers. This phase allowed the identification of main and recurrent features and the classification of all analyzed sites into three classes. For the first and most common class, a good parametrization of the velocity profile and a small variability value drive the HHS algorithm towards the fit of the H/V curve and dispersion curves at the same time (e.g. SYVP). The fit of the dispersion curves proves that the inverted velocity profiles remain close to the site characterization results, while the H/V curve provides information both for the shallow and deep layers. The second class shows sites that, independently of the allowed variability and thanks to the good initial models used, do not move the velocity profiles away from the initial information (e.g. SARK). The third group includes the sites, where the initial profile does not provide information to fit the H/V, dispersion curves and SH-amplification functions at the same time, both for low and high variabilities (e.g. STIEG). The discrepancy between the measured and inverted data in Class 3 can be explained taking into account the occurrence of lateral variations between the instrumented site and the area of the array deployment, by the quality of site characterization measurements and/or the stiffness of the site. The 23 analyzed sites are classified as follows: 8 in Class 1, 8 in Class 2 and 7 in Class 3. For all sites, we extracted the depth values reached for three shear-wave velocities, namely 1000, 1500 and 3000 m/s. The depth for Vs = 1000 m/s gives a depth range between 10 and 86 m, for 1500 m/s between 54 and 370 m and for 3000 m/s the depth ranges between 267 and 2315 m. The depth values were then used to compute the Vs based thickness maps for the study area. Of the 23 analyzed stations, 20 are located in the Swiss Molasse Basin and are inhomogeneously distributed, with higher concentration in the north-east of Switzerland than in the south-west. The final maps were then compared with the work of Sommaruga et al. (2012) but due to the uneven distribution of the stations and the unsignificant match between the two sets of data, the maps are not shown in this work.
Our analysis confirms the potential of HHS algorithm to invert the full earthquake H/V curve for subsurface investigation and its applicability to low magnitude earthquakes. We suggest exploring the parameter space and repeat the H/V inversion several times changing the allowed variabilities for the shallow layers and the half-space velocities. This step is needed to have a range of possible solutions and recognize recurrent features which can be ascribable to one of the three identified classes.
A joint inversion would be the next step towards an improved estimation of the velocity profiles: microtremor/earthquake H/V curve, Rayleigh/Lovewave dispersion curves, Rayleigh-wave ellipticity curves computed using single station methods or array techniques, and SH-wave transfer function/ empirical amplification functions.
Open access funding provided by Swiss Federal Institute of Technology Zurich. This work was funded by the SNF projects 200021_153633/1 and 200020_175475/1 and by Strong Motion Renewal Project funded by the Federal Office of the Environment (BAFU).

Data availability
The earthquake recordings analyzed in this work are open and can be downloaded using the ORFEUS website (http://orfeus-eu.org/). The inverted P-and S-wave velocity profiles can be requested to the authors.