# Profiling the Quito basin (Ecuador) using seismic ambient noise

D Pacheco, E Mercerat, Françoise Courboulex, Luis Fabian Bonilla, A Laurendeau, A Alvarado

► To cite this version:

D Pacheco, E Mercerat, Françoise Courboulex, Luis Fabian Bonilla, A Laurendeau, et al.. Profiling the Quito basin (Ecuador) using seismic ambient noise. Geophysical Journal International, 2021, 228 (2), pp.1419-1437. 10.1093/gji/ggab408. hal-03442293

HAL Id: hal-03442293

https://hal.science/hal-03442293v1

Submitted on 28 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Profiling the Quito basin (Ecuador) using seismic ambient noise

D. Pacheco<sup>1,2</sup> E. D. Mercerat,<sup>3</sup> F. Courboulex,<sup>1</sup> L. F. Bonilla<sup>1,4</sup> A. Laurendeau<sup>1,5</sup> and A. Alvarado<sup>2</sup>

Received 2021 July 26; in original form 2021 January 30

## SUMMARY

Quito, the capital of Ecuador, with more than 2.5 M inhabitants, is exposed to a high seismic hazard due to its proximity to the Pacific subduction zone and active crustal faults, both capable of generating significant earthquakes. Furthermore, the city is located in an intermontane piggy-back basin prone to seismic wave amplification. To understand the basin's seismic response and characterize its geological structure, 20 broad and medium frequency band seismic stations were deployed in Quito's urban area between May 2016 and July 2018 that continuously recorded ambient seismic noise. We first compute horizontal-to-vertical spectral ratios to determine the resonant frequency distribution in the entire basin. Secondly, we cross-correlate seismic stations operating simultaneously to retrieve interstations surface-wave Green's functions in the frequency range of 0.1–2 Hz. We find that Love waves travelling in the basin's longitudinal direction (NNE–SSW) show much clearer correlograms than those from Rayleigh waves. We then compute Love wave phase-velocity dispersion curves and invert them in conjunction with the HVSR curves to obtain shear-wave velocity profiles throughout the city. The inversions highlight a clear difference in the basin's structure between its northern and southern parts. In the centre and northern areas, the estimated basin depth and mean shear-wave velocity are about 200 m and 1800 $m^{-1}$, respectively, showing resonance frequency values between 0.6 and 0.7 Hz. On the contrary, the basement's depth and shear-wave velocity in the southern part are about 900 m and 2500 $m^{-1}$, having a low resonance frequency value of around 0.3 Hz. This difference in structure between the centre-north and the south of the basin explains the spatial distribution of low-frequency seismic amplifications observed during the $M_w$ 7.8 Pedernales earthquake in April 2016 in Quito.

**Key words:** Seismic interferometry; Seismic noise; Side effects; Surface waves and free oscillations.

## 1 INTRODUCTION

Quito, the capital city of Ecuador (Fig. 1a), is situated in a zone of high seismic and volcanic activities. With a population of 2.5 million people living in a significant percentage of dwellings that do not follow any seismic code, the city is susceptible to high seismic risk. Quito is located inland at about 180 km from the Pacific subduction zone, at 2800 m above sea level in the Inter-Andean Valley. Quito's shape is elongated in the NNE–SSW direction and has a length of approximately 40 km, while in its central area, the width (east–west) is only 5 km (Fig. 1b).

The city has been built on a piggy-back basin developed on the hanging wall of an active reverse fault system named the Quito

Fault System (QFS, Fig. 1b). The QFS is a blind N–S striking, west-dipping thrust, whose geomorphological expression at the eastern part of the city is a succession of moderate hills (Alvarado et al. 2014). In recent studies, Alvarado et al. (2014) and Mariniere et al. (2019) report that the QFS accommodates an average of 3–5 mm $yr^{-1}$ of horizontal shortening. Based on the analysis and modelling of GPS data and InSAR, a weak rate of elastic deformation accumulation and a shallow locking depth (less than 3 km) have been identified. Moreover, an aseismic slip is taking place at shallow depths in the central segment of the QFS and its associated shallow creep decreases from the centre towards the south and north segments of the fault system, suggesting a heterogeneous strain accumulation in the fault system (Mariniere et al. 2019). Beauval

<footnote n="1"><sup>1</sup>Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, Valbonne, France. E-mail: daniel.pacheco.l@icloud.com</footnote>

<footnote n="2"><sup>2</sup>Instituto Geofísico, Escuela Politécnica Nacional, Ladrón de Guevara E11-253, Quito, Ecuador</footnote>

<footnote n="3"><sup>3</sup>CEREMA, Agence de Sophia Antipolis, Valbonne, France</footnote>

<footnote n="4"><sup>4</sup>GERS Department, Université Gustave Eiffel, 77447 Marne-la-Vallée Cedex 2, France</footnote>

<footnote n="5"><sup>5</sup>Institut de Radioprotection et Sûreté Nucléaire (IRSN), PSE-ENV, SCAN, BERSSIN, Fontenay-aux-Roses, 92262, France</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Map of Ecuador with the location of the city of Quito represented by the letter Q. (b) Map of the 4B-array temporal seismic stations installed in Quito between May 2016 and July 2018. Colours in triangles indicate the different types of sensor-datalogger combinations. The lines on the map represent the surface traces of reverse faults (solid red lines), strike-slip faults (dotted-dashed red lines) and folds (dotted lines) described in Alvarado *et al.* (2014) as the Quito Fault System (QFS). The yellow stars mark the location of the two largest earthquakes instrumentally recorded in the Quito basin, according to the IG-EPN.</figcaption></figure>

*et al.* (2010) report that five earthquakes were felt in the city with MSK intensities in the VII–VIII range in the last 500 yr, most of them related to the Nazca Plate subduction zone. Nevertheless, the authors mention a crustal event in 1587 probably related to the QFS with an estimated magnitude between 6.3 and 6.5. In general, the QFS produces moderate-size earthquakes with predominantly reverse focal mechanisms (Alvarado *et al.* 2014; Vaca *et al.* 2019), and since the development of the Ecuadorian seismological network

in the early 1990s, seven earthquakes of magnitude larger than 4.0 have been recorded and located in the Quito area (Alvarado *et al.* 2018). Two of them had magnitudes greater than 5.0, one in 1990 (*M*<sub>w</sub> 5.3) and the other in 2014 (*M*<sub>w</sub> 5.1, Fig. 1b), both associated with the QFS (Beauval *et al.* 2014).

The basin is filled with volcano-sediments, lacustrine-type sediments, alluvial, colluvial and lahar deposits (Alvarado *et al.* 2014). The nature and thickness of these infilling materials are poorly
known, and studies such as Jaya (2009) and Alvarado et al. (2014) suggest that the spatial extent (horizontal and vertical) of these deposits is highly heterogeneous. Characterizing the seismic response of the Quito basin is important because basin effects can have a dramatic influence on the ground motions generated by earthquakes (e.g. Bindi et al. 2009). Guéguen et al. (2000) performed ambient seismic vibration measurements at 673 sites along the basin using three-component seismic sensors of 1 Hz. Applying the so-called horizontal-to-vertical spectral ratio (HVSR) technique, introduced by Nakamura (1989), on signals recorded at a sampling frequency of 100 Hz during 2 min, they observed one frequency peak around 1–2 Hz at most sites. Ten per cent of sites located in the basin’s central axis also presented a second peak between 5 and 8 Hz. Nevertheless, the authors could not study frequencies lower than 1 Hz because of the limited frequency band of the short-period sensors, the lack of instrumental correction in data processing, and the quite short recording time.

Recently, Laurendeau et al. (2017) analysed earthquake recordings of the National Strong Motion Network (RENAC) operated since 2010 in Quito by the Instituto Geofísico-Escuela Politécnica Nacional (IG-EPN). They identified a basin response that varies between the northern and southern parts of the city. Computing the eHVSR (earthquake horizontal-to-vertical spectral ratio) of 179 earthquakes and the Standard Spectral Ratio (SSR) of 28 earthquakes, the authors identified the main frequencies of potential site amplification at each of the 18 stations. For all the southern stations, a low-frequency peak (around 0.3–0.4 Hz) is present in both spectral ratios (eHVSR and SSR). In contrast, in the stations located in the city’s central and northern parts, the spectral ratios remain relatively flat at frequencies below 2 Hz. This different seismic response in the northern and southern part of the basin was confirmed by the largest earthquake recorded by the network (Pedernales earthquake, 16 April 2016, M<sub>w</sub> 7.8), with larger amplitudes and much longer durations in the southern part of the city. These results suggest that the bedrock interface and overlying sediments are far from being continuous along the basin. Apart from these observations, the fundamental parameters controlling the seismic amplification, such as the thickness and the shear-wave velocities of the infilling material and the seismic bedrock’s geometry, remain unknown.

The high population density and building concentration make geotechnical (boreholes) and geophysical (active seismic) studies rather difficult and expensive to carry out. When performed, they only explore the first few tens of metres of soil. In such cases, passive seismic methods based on continuous recordings of ambient seismic noise by temporary networks installed in urban areas are attractive because of the field logistic’s simplicity. Therefore, they represent an alternative and powerful tool to investigate both shallow and deeper geological structures, from tens to thousands of metres (e.g. Ma & Clayton 2016; Manea et al. 2016; Vassallo et al. 2019).

It has been theoretically shown that the cross-correlation of ambient noise recorded at pairs of seismic stations can provide an estimate of the Green's functions between receivers (e.g. Larose et al. 2004; Wapenaar 2004). This technique has been proven by several applications in different areas of the world and different geological environments, from global to regional and local scales (e.g. Shapiro et al. 2005; Bensen et al. 2007; Zheng et al. 2008; Kästle et al. 2018). Since not long ago, ambient noise cross-correlation studies in urban environments were used to obtain the geometry of the bedrock-sediment interface and the shear-wave velocity profiles in densely populated cities located on sedimentary basins. In Jakarta-Indonesia, for example, using Green's functions recovered

from seismic noise, Saygin et al. (2015) performed a Rayleigh wave tomography and developed a high-resolution shear-wave velocity model, in the frequency range of 0.2–2 Hz. Ma & Clayton (2016), applying ambient noise and receivers function analysis, derived a shear-wave velocity and structural model for Los Angeles basin. Using a dense network of sensors, they found that the depth of the sediments in Los Angeles varies from 4 to 8 km deep. On a smaller scale, Pastén et al. (2016) using surface-wave group and phase velocities obtained from ambient noise cross-correlations, characterized the geological structure under three zones in the Santiago basin (Chile). Similarly, in the urban area of Benevento city (southern Italy), Vassallo et al. (2019) cross-correlated 1 month of ambient seismic noise to reconstruct Rayleigh waves travelling between pairs of stations in the basin. They inverted phase-velocity dispersion curves to estimate 1-D representative shear-wave velocity profiles, finding two main seismic interfaces at a depth of about 250 and 400 m in the basin.

Because seismic noise in cities is predominantly driven by numerous processes such as cultural life, traffic and the production and transportation of all kinds of goods, much of the recurring and time-dependent seismic signals are generated and emitted by human-made physical processes. These anthropic ground movements overlap with natural ground vibrations (wind-induced tremors, ocean waves or earthquakes) to form highly variable spatial and temporal seismic environmental noise. These characteristics make the analysis of recorded ambient noise within cities particularly challenging.

Inspired by these previous studies and with the aim to image the depth of the seismic bedrock and determine the seismic wave's velocities, we deployed a temporal seismic network throughout the Quito basin from May 2016 to July 2018 (Table 1). This network installation was part of a large research project involving several French and Ecuadorian institutions (Cerema, Géoazur and IG-EPN) to study the hazard and seismic risk.

In this paper, we present the ambient noise experiment results in the city of Quito, which focuses on the location of the sediment-bedrock interface that plays a critical role in seismic amplification and its spatial variability. In the coming sections, we briefly introduce the Quito seismic experiment and the characteristics of recorded seismic noise in the urban environment. We calculate HVSR curves and determine the fundamental frequency ($f_0$) per seismic station and its spatial variation throughout the basin. The cross-correlation computation of continuous ambient noise for all seismic station pairs is carried out, and we retrieve Love wave phase-velocity dispersion curves along the basin. Finally, we jointly invert HVSR and phase-velocity dispersion curves to obtain 1-D shear-wave velocity models for ten seismic stations located throughout the city. The results introduced by this research are the first step towards the detailed structure of the Quito basin we attempt to obtain in the future.

## 2 QUITO SEISMIC EXPERIMENT

Between May 2016 and July 2018, 20 broad and medium frequency band three-component seismological stations were deployed progressively throughout the city to record ambient seismic noise for a year or more (referred from here as the 4B-array, Mercerat et al. 2016). The interstation distances of the 4B-array varied from 1.8 to 32 km (see Fig. 1b and Table 1). The altitude difference between the highest (ROQE) and lowest (PUEM) stations was 780 m. The stations located in the south are, on average, 100 m higher than the stations located in the northern part of the city (Table 1).
<table><caption>Table 1. Description of the 4B-array.</caption><thead><tr><th>Station</th><th>Lon.</th><th>Lat.</th><th>Alt.</th><th>Start date</th><th>End date</th><th>Digitizer</th><th>Samp. freq.</th><th>Sensor</th><th>Corner freq.</th></tr></thead><tbody><tr><td>CRON</td><td>-78.42</td><td>-0.06</td><td>2801</td><td>Dec-2016</td><td>Jul-2018</td><td>Agecodagis Kephren</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>MAST</td><td>-78.46</td><td>-0.07</td><td>2648</td><td>Jul-2017</td><td>May-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>HCAL</td><td>-78.43</td><td>-0.09</td><td>2721</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>DELA</td><td>-78.49</td><td>-0.11</td><td>2785</td><td>May-2016</td><td>Jul-2018</td><td>Agecodagis Kephren</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>PUEM</td><td>-78.36</td><td>-0.13</td><td>2404</td><td>Jul-2017</td><td>Still working</td><td>Agecodagis Kephren</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>SOCA</td><td>-78.46</td><td>-0.13</td><td>2888</td><td>Dec-2016</td><td>Nov-2018</td><td>Reftek RT 130</td><td>200</td><td>Trillium Compact</td><td>0.008</td></tr><tr><td>INCA</td><td>-78.47</td><td>-0.15</td><td>2837</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>PARU</td><td>-78.50</td><td>-0.17</td><td>2919</td><td>Dec-2016</td><td>Still working</td><td>Reftek RT 130</td><td>200</td><td>Trillium Compact</td><td>0.008</td></tr><tr><td>JBPC</td><td>-78.48</td><td>-0.18</td><td>2788</td><td>May-2016</td><td>Aug-2018</td><td>Agecodagis Kephren</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>INAQ</td><td>-78.46</td><td>-0.19</td><td>2908</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>ROQE</td><td>-78.53</td><td>-0.21</td><td>3184</td><td>May-2016</td><td>Jul-2018</td><td>Agecodagis Kephren</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>OBSA</td><td>-78.50</td><td>-0.21</td><td>2824</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>CSEK</td><td>-78.44</td><td>-0.20</td><td>2876</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>HSUR</td><td>-78.54</td><td>-0.24</td><td>2908</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>QUEM</td><td>-78.49</td><td>-0.23</td><td>3038</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>ARGE</td><td>-78.52</td><td>-0.27</td><td>2911</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>QUIB</td><td>-78.54</td><td>-0.29</td><td>2918</td><td>May-2016</td><td>Jul-2018</td><td>Agecodagis Kephren</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>HLUZ</td><td>-78.55</td><td>-0.31</td><td>2980</td><td>Jul-2017</td><td>Jul-2018</td><td>Guralp CMG DM24</td><td>200</td><td>Guralp CMG40T</td><td>0.2</td></tr><tr><td>SVAC</td><td>-78.49</td><td>-0.32</td><td>2573</td><td>Jul-2017</td><td>Still working</td><td>Reftek RT 130</td><td>200</td><td>Trillium Compact</td><td>0.008</td></tr><tr><td>GARO</td><td>-78.53</td><td>-0.33</td><td>3011</td><td>Oct-2016</td><td>Jul-2018</td><td>Reftek RT 130</td><td>200</td><td>Trillium Compact</td><td>0.008</td></tr></tbody></table>

The array was equipped with different three-component seismic instruments: 16 stations consisted of medium frequency band velocity sensors (Guralp CMG-40T 5s), 10 of them had a 24-bit ADC Guralp acquisition and digitization system (CMG DM24) and 6 of them had a 24-bit ADC Agecodagis acquisition and digitization system (Kephren). Four stations were equipped with a broadband velocity sensor (Nanometrics Trillium Compact 120s) and a 24-bit ADC Reftek acquisition and digitization system (RT 130). The velocity sensors were buried 50–80 cm below the ground surface, placed over a 3-cm-thick granite plate, and entirely covered by the ground soil. All digitizers were continuously synchronized using GPS antennas, and data were recorded continuously at a sampling frequency of 200 Hz and collected from the field every 3 months maximum. The acquired raw data was then converted and organized into a 1-day-long mseed database, composed of 9497 three-component (north, east and vertical) traces.

## 3 SEISMIC NOISE ANALYSIS

We compute hourly power spectral densities (PSD) and their probability density functions (PDF), to characterize the background seismic noise level in the city of Quito. Continuous three-component (north, east and vertical) seismic data acquired in the 4B-array is analysed following the methodology detailed in McNamara & Buland (2004) and Nakata et al. (2019). We parse 1-day-long files, for each station component, into 1-hour time segments. To further reduce the final PSD estimate variance, each 1-hr segment is divided into 13 subsegments, overlapped by 75 per cent.

We compute the spectrum in each of the 13 subsegments; individual spectrums are then averaged together to estimate the hourly PSD. We remove the nominal instrumental response to attain hourly PSD in ground motion units, expressed in decibels (dB) with respect to acceleration (0 dB = 1 m²/s⁴/Hz). Acceleration PSD is a standard format for referencing noise at a seismic station to both global background noise models (Peterson 1993) and self-noise models of seismic instrumentation (e.g. Ringler & Hutt 2010). Once the hourly PSDs are calculated, to reduce the number of samples, we

use one octave of data around each extracted frequency at one octave intervals. The daily PDFs are an aggregate of 47 PSDs, each representing 60 min of ground motion and overlapping with the next hourly PSD by 30 min. We obtain the daily representative spectrum (statistical mode) by extracting the most frequent noise level for each frequency from daily PDFs.

To fully characterize the seismic noise levels, we generate long-term PDFs from hundreds of hourly PSDs processed using the method described before for each station component. The long-term representative spectrums are computed from long-term PDFs. New Low Noise Model (NLNM) and New High Noise Model (NHNM, Peterson 1993) are used as a reference to assess the quality of the chosen sites. For the computations, we use an open-source bundle of three highly configurable Python scripts (Noise Toolkit NTK from Hutko et al. 2017).

Fig. 2 shows a comparison of the long-term PDFs for four seismic stations of the 4B-array. For frequencies between 0.1 and 1 Hz, the representative spectra of long-term PDFs are within the zone defined by NLNM and NHNM levels for all the stations. On seismic stations composed of medium frequency band sensors (T = 5 s) such as ARGE, DELA or PUEM (Fig. 2), the instrumental noise increases at frequencies lower than 0.1 Hz, and the horizontal components overcome the NHNM level for frequencies less than 0.05 Hz. Conversely, for the stations with broad-band sensors like PARU (Fig. 2), the power level in the primary microseism band (0.05–0.07 Hz) is well under the NHNM level.

For frequencies greater than 1 Hz, all stations exhibit noise levels higher than the NHNM level (Figs 2 and 3), which have been related to the anthropic activity in highly urbanized areas (McNamara & Buland 2004; Boese et al. 2015; Díaz et al. 2017). To analyse seismic noise subdaily variations, such as anthropogenic influences due to night-day differences, we accumulate the PSDs in hourly bins and compute PDFs for each hour of the day over a year (July 2017–July 2018) in all station components. We also plot the PDF of representative spectra as a function of the hour [in Coordinated Universal Time (UTC) and local time (LT)] at stations DELA and PUEM (Fig. 3). At lower frequencies (<1 Hz), the noise levels in the horizontal components (north and east) are higher than in the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Long-term PDF of the continuous signal recorded on four three-component stations. The stations are located at the south (ARGE), at the north (DELA), at the west (PARU) and outside (PUEM) of Quito city. Long-term representative spectrums (statistical mode, dotted lines) of PDF are indicated and dashed bounding curves show low- and high-noise models (NLNM and NHNM; Peterson 1993).</figcaption></figure>

vertical (Figs 2 and 3). We observe this behaviour in all 4B-array stations in this frequency band.

Between 1 and 10 Hz, the night-day difference is in the order of 15–20 dB (Fig. 3). This pattern of increased noise during the daylight working hours is observed at every station in the array. The level of seismic noise and the duration of the quiet period vary between stations. For example, in Fig. 3, we observe a decrease in the noise level in PUEM between 01- and 10-hr UTC (20- and 05-hr LT) for frequencies higher than 5 Hz, while in DELA, the seismic noise decreases between 04- and 08-hr UTC (23- and 03-hr in LT) for frequencies greater than 1 Hz. We compute monthly and

weekly PDF for all the stations, but no evident seasonal variations are observed.

The frequency range between 0.01 and 0.1 Hz is not of interest to our study because of the very high instrumental noise, mainly observed in the seismic stations with Guralp CMG-40T 5 s sensors (Figs 2 and 3). Ambient noise from oceanic sources dominates the frequency band between 0.1 and 1 Hz, even with the 4B-array seismic stations installed in a highly urbanized area. In this frequency band, representative PDF spectra is within the NLNM and NHNM curves calculated by Peterson (1993). The lack of seasonal variations indicates that oceanic noise sources were stable during the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Night-day variations for DELA and PUEM. The dotted black line indicates the frequency of the second microseism.</figcaption></figure>

study period. The cultural noise (e.g. Boese *et al.* 2015; Díaz *et al.* 2017) dominates the frequency spectrum at high frequencies (approximately >1 Hz). In this range, the spectrograms clearly show diurnal and weekly variations. Due to the issues discussed herein, we decide to use both oceanic and coherent cultural noise in the frequency band between 0.1 and 10 Hz for this study.

## 4 HORIZONTAL-TO-VERTICAL SPECTRAL RATIOS (HVSR)

To obtain the resonant frequency for each seismic station of the 4B-array, we use the well-known HVSR technique. HVSR provides a measure of the fundamental resonance frequency at a given site (e.g. Nakamura 1989; Castellaro 2016; Molnar *et al.* 2018; Nakamura 2019) and works particularly well in the case of 1-D plane-parallel stratigraphy.

Similar to the approach described in McNamara *et al.* (2015), we use available PDFs to compute HVSRs using daily and long-term representative spectrums. We then calculate the spectral ratio of the horizontal by vertical representative spectra derived from daily and long-term PDFs at a single three-component seismic station. Because two perpendicular (north and east) ground motion components are measured at the seismic stations, to combine the horizontal Fourier spectra, we use the Diffuse Field Assumption Method (DFA) proposed by Sánchez-Sesma (2011) (M1), and also the methods referenced by Albarello & Lunedei (2013): arithmetic mean (M2), geometric mean (M3), vector summation (M4), quadratic mean (M5) and maximum horizontal value (M6). All calculations were done using the IRIS Toolbox proposed by Bahavar *et al.* (2020). By using six HVSR calculation methods (Fig. 4b), we expect robust results

in determining the resonance frequency at each site. To determine the resonance frequency, we identify all significant peaks ($f_s$) in the spectral ratio curves (larger than 1.5), and then choose the lowest frequency peak ($f_0$) as the site resonance frequency. We compute HVSR in the 0.1–10 Hz frequency range because Guéguen *et al.* (2000) and Laurendeau *et al.* (2017) identified the main resonance frequency for the Quito basin within this frequency band.

Fig. 4(a) shows the long-term PDFs representative spectra calculate for the north, east and vertical components for ARGE (south of the city), DELA (north of the city) and PARU (west of the city) stations. The long-term HVSR curves calculated using the six methods described above are shown in Fig. 4(b). In each example, the shape of the spectral ratio curves is similar (Fig. 4b), and differences in the location of HVSR peaks in the frequency axis are negligible; however, M1 and M4 methods show higher amplitudes in the HVSR curves than the other methods (Fig. 4b). The larger HVSR values for M4 are due to the $\sqrt{2}$ factor in the technique (magnitude of the sum of two-unit orthogonal vectors $\sqrt{2}$, Albarello & Lunedei 2013). These larger M4 values are comparable to M1 values because the PSD of the Fourier transform of the *i* component of displacement used in the DFA approach behaves as vectors of averages (Sánchez-Sesma 2011; Bahavar *et al.* 2020). In most of the seismic stations of the 4B-array, we observe that the shape of the daily HVSR curves does not change substantially over the year (as an example, see ARGE and PARU stations in Fig. 5), and the value of $f_0$ is stable over time (Fig. 5).

To investigate the azimuthal variability in the HVSR curves, we rotate the two horizontal components (north and east) to obtain a representative horizontal time-series by 10° steps, starting from 0°N to 180°N. We first compute hourly rotated HVSR curves at each
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Long-term representative spectrums computed for the three-component stations ARGE, DELA and PARU (a), together with the long-term HVSR curves calculated with six methods (b). For the location of seismic stations, see Fig. 1. M1 = Diffuse Field Assumption Method, M2 = Arithmetic Mean, M3 = Geometric Mean, M4 = Vector Summation, M5 = Quadratic Mean and M6 = Maximum Horizontal Value. Shaded areas behind HVSR curves represent the standard deviation of each one.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Daily HVSR curves computed for ARGE and PARU stations for 1-yr period. Difference panels show the computed methods of HVSR.</figcaption></figure>

station to observe variations related to the anthropic activity. The polar diagrams present in Fig. 6(b) show the azimuthal variation of the HVSR curves for ARGE at 07:00 (UTC) and 17:00 (UTC), hours in which there is the lowest and highest anthropic activity in Quito, respectively (Fig. 3). Although a slight difference is observed in curves' amplitudes, no significant changes in shape are present. Moreover, when we compare the long-term HVSR curve calculated in ARGE (M3 in Figs 4b and 6c) with the azimuthal average HVSR curves of 07:00 and 17:00, the variations observed are minor (Fig. 6c). In a more general view, when we calculate and average the rotated daily HVSR curves for each station, directional effects are not evident in the polar diagrams (Fig. 6d). The *f*₀ pick recorded in the stations is almost independent of the azimuth, and

at the same time, the amplitude of the *f*₀ peak does not vary much, for example, in SOCA (north), ROQE (centre) and HLUZ (south); while HCAL (north), INAQ (centre) and ARGE (south) *f*₀ amplitude variations with respect to azimuth are observed. Nor are directional effects observed in stations located outside the Quito basin, such as PUEM or SVAC (Fig. 6).

Considering that *f*₀ is independent of the time window, azimuth, or the method used for the HVSR computations, we map *f*₀ for the seismic stations throughout Quito to study the resonant frequency's spatial variability (Fig. 7). In the southern part of the city (ARGE, QUIB and HLUZ), we estimate *f*₀ near 0.3 Hz, while in the basin's centre (DELA, INCA, OBSA, QUEM and HSUR), *f*₀ is about 0.6 Hz. SOCA and PARU have *f*₀ near 0.7 Hz. GARO and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Map with the location of the seismic stations present in sections b to d. (b) Rotated HVSR polar diagrams of ARGE at two different hours, where 07:00 (UTC) shows the time with low anthropic activity and 17:00 (UTC) the time with high anthropic activity. (c) Comparison between the long-term HVSR curve computed in ARGE with the mean 07:00 and 17:00 HVSR curves. Shaded areas behind HVSR curves represent the standard deviation of each one. (d) Rotated HVSR polar plots diagrams of HCAL, SOCA and INAQ (north), with OBSA and ROQE (centre), ARGE and HLUZ (south) and PUEM and SVAC (outside Quito basin).</figcaption></figure>

SVAC, located in city's southeastern part, shows $f_0$ at 0.8 Hz. In the northern part, MAST has $f_0$ at 0.4 Hz. In CRON and HCAL, we estimate $f_0$ around 1.2 Hz, and in ROQE (eastern flank of the Pichincha volcano, Fig. 1b) at 2.6 Hz. PUEM, located outside the Quito basin, has a complex shape in its HVSR curve, without a clearly identifiable peak. At frequencies higher than 2 Hz, we also find high-frequency peaks ($f_p$) around 3 Hz in the JBPC, ARGE, HLUZ and SVAC. Finally, a west-east profile formed by JBPC, INAQ and CSEK shows rather complex shapes and the presence of several peaks. In JBPC, no clear peak is identified, while INAQ presents $f_0$ around 0.5 Hz, and in CSEK, the $f_0$ is near 0.4 Hz (Fig. 7a).

To summarize all these results, we present in Fig. 7(a) an image of the long-term HVSR curves obtained using the M4 method for all the stations ordered from north to south, and we mark $f_0$ at each

station (black vertical lines in Fig. 7a). The seismic stations located in the western part of the basin (the western side in Fig. 7a) and those located at the east (the eastern side in Fig. 7a) do not show any evident variation of $f_0$ as a function of latitude. In the central part of the basin (central profile in Fig. 7a), the value of $f_0$ clearly decreases southward (ARGE, QUIB and HLUZ). Applying a linear interpolation, we use the measured $f_0$ values to draw up an iso-resonance frequency map. To improve the results of the linear interpolation, we utilize the $f_0$ values from the strong-motion (accelerometer) stations VILF (0.41 Hz), ZALD (0.34 Hz) and LILI (0.36 Hz) published in Laurendeau *et al.* (2017). These $f_0$ values were obtained using eHVSR and SSR methods. These three strong-motion stations are located in Quito's southern zone and complement and densify the values of $f_0$ in this sector (Fig. 7b).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) HVSR curves calculated for each station. The curves are organized in function of the station's latitude and distributed in three station profiles. Black solid lines mark the value of $f_0$, and the colour scale highlights the different amplitudes of HVSR curves. (b) map of the linear interpolation of the $f_0$ values shown in the left-hand panel plus the $f_0$ values for the VILF, ZALD, LILI accelerometer stations taken from Laurendeau *et al.* (2017).</figcaption></figure>

# 5 MULTICHANNEL CROSS-CORRELATION ANALYSIS

In addition to the previous *single-station* analysis of ambient noise data recorded by the 4B-array, we apply a *multistation* analysis to fully exploit the data's information. Cross-correlation based analysis is a powerful tool to extract information about wave propagation inside the medium and, if possible, invert dispersion data to obtain shear wave velocities. The spatial and temporal coverage of the 4B-array represents the first opportunity to carry out this type of study in the Quito basin.

## 5.1 Data processing

We analyse continuous data from the 20 stations of the 4B-array recorded between May 2016 and July 2018 to compute cross-correlations between all possible pairs of station-components.

To estimate the interstation Green's functions, we follow the procedure introduced by Bensen *et al.* (2007) with some modifications. We divide continuous 1-day-long data into time segments with a duration of 60 min (3600 s). We then remove the mean and the linear trend before applying an instrumental response correction. We down sample each trace to 50 Hz. We then apply a spectral whitening to the 1-hr traces (1) to equalize their spectral amplitude content, (2) to minimize the contribution of the most energetic signals present in the data and the possible persistent monochromatic sources and (3) to broaden the frequency band of the surface waves dispersion measurements (Lin *et al.* 2008). We attempt to suppress the effects of earthquakes and local non-stationary noise sources around the stations. Instead of applying the traditional 1-bit

normalization, we used an automatic gain control (AGC) operator with a window length of 5 s proposed by Behm *et al.* (2013). We perform hourly cross-correlation in the time domain between all possible station-components, using time lags between -200 and 200 s, because we expect to observe the highest dispersed surface-wave train at the longest inter-station distance (32 km) in this time window. To enhance the coherent signals, we stack the hourly cross-correlation functions for every station pair using two algorithms: (1) a linear time stack (e.g. Bensen *et al.* 2007; Lin *et al.* 2008) and (2) a non-linear timescale phase-weighted stack (ts-PWS, Ventosa *et al.* 2017), that uses phase coherence to improve the coherent signals of the cross-correlation functions.

Figs 8(a) and (e) show a comparison of the two algorithms (linear stack and ts-PWS) for the cross-correlation functions calculated between DELA and INCA (inter-station distance = 5.75 km). We see how the causal (positive lag time) and acausal (negative lag time) signals emerge as the time-series' length increases from a 1-day to 1-year stack. The emergence of coherent signals as a function of stack time is particularly well exemplified in the stacked cross-correlations' amplitude spectra. With increasing stack time, the spectra's energy increases in the frequency bands between 0.1-0.5 Hz and 0.5-2 Hz (Figs 8b and f). The spectrum of the 1-yr signal corresponding to the ts-PWS (Fig. 8f) clearly shows these two frequency band peaks. The non-coherent high-frequency amplitudes decrease drastically comparing 1-yr ts-PWS with 1-day or 1-week ts-PWS (Figs 8e and g).

To quantify signal emergence with increasing time-series length in cross-correlations, we define the 'signal' as the maximum amplitude peak measured between lag times 0–25 s and the 'noise' as the root-mean-square (RMS) of the amplitude values between lag times
<figure><img src="image_10.png" /><figcaption>**Figure 8.** DELA-INCA north-component example of the signal's emergence for increasingly long time-series using the traditional linear stack scheme (a) and ts-PWS (e). In the amplitude spectra of the stacked cross-correlations, using linear stack (b) and ts-PWS (f), it is possible to observe how the coherent signals' energy increases as a stacking time function. The filtered stacked cross-correlations in the 0.1–0.5 and 0.5–2 Hz bands are also shown (c and g). To calculate the signal-to-noise ratio, the cross-correlation functions (circles) were used, as well as the filtered stacked cross-correlations between 0.1–0.5 Hz (triangles) and 0.5–2 Hz (squares) for the linear stack (d) and ts-PWS (h).</figcaption></figure>

30–50 s. The resulting ratio (SNR) per frequency band is shown in Figs 8(d) and (h). The SNR values of the stacked signals, obtained from cross-correlations filtered between 0.1 and 0.5 Hz (triangles in Figs 8d and h), are about an order of magnitude greater using the ts-PWS compared with the linear-stack. In this frequency band, the SNR value observed using the 1-month time-series is almost the same as that observed for time-series of 6-months or 1-yr. For cross-correlations filtered between 0.5 and 2 Hz (squares in Figs 8d and h), the SNR values using the linear stack algorithm are low, even using 1-yr of data.

Since ts-PWS is a non-linear stacking method, it can extract a signal, discarding the incoherent noise between 0.5 and 2 Hz, mainly present in an urban context (traffic noise), and clearly exhibits higher SNR values.

We obtain the best results using the ts-PWS algorithm (Figs 8e to h). In general, for the cross-correlation functions of the stations of the 4B-array in the frequency range between 0.1 and 0.5 Hz a coherent signal emerges from the stack of 1 week, and it is stable after 1 month of data. Between 0.5 and 2 Hz, a coherent signal emerges from 1 month, and it is stable after stacking 6 months of data. As exemplified in Fig. 8(f), coherent signals only were observed for frequencies between 0.1 and 2 Hz, even when using 1-yr data. As a final processing step, the nine components of the Green's tensor

were obtained by rotating the original nine components (EE, EN, EZ, NE, NN, NZ, ZE, ZN, ZZ; E = east, N = north, V = vertical) using the procedure described in Lin *et al.* (2008). The first and second characters in the components of the Green Tensor (RR, RT, RZ, TR, TT, TZ, ZR, ZT and ZZ; R = radial, T = transverse, Z = vertical, Fig. 9) represent the direction of the source-side and the receiver-side stations, respectively.

## 5.2 Cross-correlations features in the urban area of Quito

Spectral analysis for the nine components of the Green's tensor shows that coherent signals that emerged from ts-PWS have energy mainly between 0.1 and 2 Hz (Fig. 9), for this reason, we filter the Green's tensor in this frequency band. Clear wave propagation patterns are identified in this plot at both positive and negative correlation lags. As expected, the RR, RZ, TT, ZR and ZZ components show much more consistent wave energy than the others. For most station pairs, the predominant signal visible in the TT component is outstanding (Figs 9 and 10).

In the ZZ component, a large amplitude signal is observed around $t = 0$ up to an inter-station distance of 25 km (Fig. 9). Clear wave-forms near $t = 0$ are often observed in previous ambient noise correlation studies (e.g. Villaseñor *et al.* 2007; Poli *et al.* 2011;
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Interstation distance (km) as a function of the correlation time (s) for the nine components of the correlation tensor (the components are indicated above the panels). Correlation functions were filtered between 0.1 and 2.0 Hz, and they were stacked for each 1 km distance bin. Colour scale indicates positive and negative normalized amplitudes, with the same scale in all panels. Record sections for every combination of three-component motion are labelled as follows: R = radial, T = transverse, Z = vertical. RR, RZ, ZR and ZZ components show Rayleigh waves, and TT shows clear Love waves.</figcaption></figure>

Taylor *et al.* 2019) though the physical reason is not completely clear. In our case, they are possibly associated with the ocean generated ambient seismic noise coming from the Pacific or the Atlantic oceans nearly with an EW azimuth and therefore arriving almost at the same time to all stations in the basin. We think this is a viable explanation because the dominant frequency for this zero-lag peak is around 0.3 Hz. However, teleseismic body wave energy that arrives at the stations at a near-vertical incidence angle could also explain the zero-lag peaks in the vertical component's correlation functions (e.g. Landès *et al.* 2010; Hillers *et al.* 2013). Rayleigh waves are observed in the cross-correlation functions between the vertical and radial components (RR, RZ, ZR and ZZ), whereas Love waves are observed in TT components (Fig. 9).

We use the RR, TT and ZZ cross-correlation functions obtained by cross-correlating JBPC with all other stations to analyse the seismic wave propagation along the Quito basin. JBPC is located approximately in the centre of the 4B-array (Figs 1 and 10), with distances to the northernmost (CRON) and southernmost (GARO) stations around 15 and 17 km, respectively. In Fig. 10, we take the cross-correlation functions presented for JBPC in Fig. 9 (Fig. 10a) and filter them between 0.1–0.5 Hz (Fig. 10c) and 0.5–2 Hz (Fig. 10d). The acausal part of the cross-correlation functions corresponds to energy propagating from JBPC, and the causal parts represent the energy propagating towards JBPC. In general, we observe two groups of waves travelling across the basin, one at low frequencies (0.1–0.5 Hz, sections b and c in Fig. 10), and another at high frequencies (0.5–2 Hz, sections b and d in Fig. 10). In the case of low-frequency waves, the acausal and causal parts are symmetrical and without significant variations in amplitude or phase (panel c in Fig. 10). On the other hand, significant asymmetry is observed in the amplitude and spectral content of acausal and causal

signals with frequencies above 0.5 Hz (section d in Fig. 10). This behaviour suggests that the distribution of noise sources is frequency-dependent. Waveforms in the TT component are much clearer than waveforms in the RR or ZZ components, especially at frequencies higher than 0.5 Hz (sections b and d in Fig. 10). The strongest signals in the TT component's causal parts between QUEM and JBPC and TT component's acausal parts north of JBPC may indicate that the noise propagation is more efficient in the south-to-north direction for Love waves at frequencies higher than 0.5 Hz (section d in Fig. 10). This behaviour is also dimly observed in the RR and ZZ components. In this frequency band, we do not observe clear, coherent signals south of QUEM in any component (section d in Fig. 10).

To invert the dispersion curves and further suppress the incoherent signals, we average the positive and negative sides of the cross-correlation functions between 0.1 and 2 Hz to obtain what is called the symmetric signal. This operation has been commonly used in various applications of the cross-correlation technique (Bensen *et al.* 2007; Lin *et al.* 2008; Asano *et al.* 2017; Vassallo *et al.* 2019) in order to better approximate the Green's function between stations pairs. To perform the time-frequency analysis, we focused on the symmetric cross-correlation functions. In Fig. 11, we present amplitude cross-correlation spectrograms corresponding to the 10-station profile that is approximately aligned with basin's elongation and take JBPC as a virtual source. The term virtual source refers to the principle of reconstructing seismic responses by cross-correlating seismic wavefield records obtained at different receiver locations. The function or waveform obtained by cross-correlating the two receiver recordings can be interpreted as the response that would be measured at one of the receiver locations as if there was a source at the other location (e.g. Wapenaar *et al.* 2010). The symmetric
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Acausal and causal parts of the obtained inter-station cross-correlation functions for JBPC and all others in the 4B-array in three frequency bands. The map at the top (a) indicates the location of the seismic stations and in red JBPC. Red traces represent the autocorrelation of JBPC. The RR, TT and ZZ were filtered between 0.1–2 Hz (b), 0.1–0.5 Hz (c) and 0.5–2 Hz (d).</figcaption></figure>

cross-correlation signal represents, as best as we can, the interstation Green's function mainly composed of surface waves travelling from one seismic station to the others (e.g. Wapenaar *et al.* 2010).

In general, clear dispersion of surface waves is observed in the TT component north of QUEM (central part of the basin). This dispersion effect suggests a satisfactory recovery of Love waves dispersion. On the other hand, waveforms with dispersive characteristics are relatively rare in RR and ZZ components, as we can

readily see in Figs 9, 10 and 11. In the RR component, we observe the dispersion effect at distances less than 4 km to the north of JBPC and at distances less than 6 km to the south of JBPC (Fig. 11). A signal of approximately constant frequency content around 0.3 Hz dominates the spectrograms of the ZZ component, and the energy at higher frequencies is almost absent (Fig. 11). In the spectrograms, frequencies higher than 0.5 Hz are not observed for stations south of QUEM (Fig. 11). This affects the possibility to perform satisfactory
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Time frequency analysis of the symmetric Green's function for the 10-station which are approximately aligned with basin's elongation and that take JBPC as a virtual source. The graphical scale shows the amplitude normalized with respect to the maximum amplitude of each symmetric signal.</figcaption></figure>

surface waves group velocity measurements. Signals with frequencies lower than 0.5 Hz are present in all spectrograms. The duration they reach (around 15–20 s) in the symmetric cross-correlation functions between JBPC and ARGE, QUIB, HLUZ and GARO stations is outstanding (Fig. 11). This may indicate the presence of low frequency waves trapped in the southern part of the basin, which deserves more attention and research in the future. The attenuation of frequencies higher than 0.5 Hz in the cross-correlation functions with ARGE, QUIB and HLUZ stations may also be related to structural changes (depth of the basement) or changes in the anelastic properties of the sedimentary infill. Love waves on the TT components protrude, and they exhibit clear phase velocity dispersion and present much higher SNR than Rayleigh waves (RR or ZZ components). Due to the evidence described herein, we decided to use the symmetric signals of the TT component in the next section when inverting the velocity structure of the basin.

## 6 1-D SHEAR-WAVE VELOCITY PROFILES FROM JOINT INVERSION OF PHASE VELOCITY DISPERSION AND HVSR CURVES

Due to the elongated shape of the Quito basin and the seismic stations' distribution in the 4B-array, we use the symmetric cross-correlation signals of the eleven stations analysed in Fig. 11, together with the profile formed by the PARU, JBPC, INAQ and CSEK stations, to obtain, as a first step, phase velocity dispersion curves. This 11-station profile has the advantage of having seismic stations distributed from north to south and is more or less aligned with the direction of maximum elongation of the Quito basin. We divide it into three virtual source gathers that traverse the city (Fig. 12).

A section of symmetric cross-correlation traces in the north part of the basin, with a virtual source at CRON (Line-1 on map in
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Map with the location of the four virtual source gathers, together with the symmetric signals of the virtual source gathers, and the dispersion images calculated from the symmetric signals are shown. The bottom panel shows the best solutions for inverting the picked dispersion curves in the dispersion images for each of the virtual sources gathers. The dashed red lines in the *V*<sub>s</sub> profiles represent the best solutions. The dashed red lines on the dispersion images show the best theoretical dispersion curve, and the green dots represent the manually picked points. On the map, the colour circles represent the receivers, and the red ones represent the virtual sources. To facilitate visualization in the map, ellipses were drawn around each seismic source gather in the map. The amplitudes in the Dispersion Images are normalized in relation to the maximum, and to facilitate visualization, only values greater than 0.5 have been plotted.</figcaption></figure>

Fig. 12), is compared with gathers located in the centre (Lines 2 and 4 on the map in Fig. 12), and in the south of the basin (Line-3 on the map in Fig. 12). Lines 1, 2 and 3 are approximately parallel to the basin's NNE-SSW elongation, while Line 4 perpendicularly transects the basin near its middle zone. In Fig. 12, the symmetric cross-correlation sections of the four lines are shown in the waveforms panel. Traces are arranged according to their distance (offset) between the virtual source (red dots on the map in Fig. 12) and the receiver (colour dots on the map in Fig. 12). The azimuths between virtual source and receivers for Lines 1, 2, 4 are relatively homogeneous, with variations less than 15°. Line 3 (covering the southern part of the basin) presents slightly higher differences in interstation azimuths. For the calculation of phase dispersion curves (dispersion images panel in Fig. 12), we compute the Love wave phase velocity spectra applying a frequency-phase velocity

transform. This technique consists of the summation, at each frequency, of wave fields present in all traces travelling with a selected phase velocity. This operation produces a so-called dispersion image or dispersion spectrum (Park *et al.* 1999; Park 2011).

Clear phase dispersion trends are observed between 0.15 and 1.25 Hz for Lines 1, 2 and 4. On the other hand, for Line 3, the coherent energy in the dispersion image fades for frequencies higher than 0.35 Hz. We manually picked the phase dispersion curves for the fundamental mode (green dots in dispersion images panel in Fig. 12), looking for the maximum amplitude around each selected point in Love wave phase velocity dispersion images. Manually picked dispersion curves in Lines 1, 2 and 4 are more or less similar, with smooth curves between 0.15 and 1.25 Hz. The dispersion curve in Line 3 (southern zone) shows a different shape, with a rapid increase in phase velocity between 0.2 and 0.35 Hz. Parametrization
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Shear wave models (c) obtained from the best fit of the dispersion curves (a) and HVSR (b) inversion functions. The green lines show the measured curves, while the dashed red lines are the best models.</figcaption></figure>

and Neighbourhood algorithm (NA) inversion (Wathelet et al. 2004; Wathelet 2008) were performed using Geopsy open-source software package (Wathelet et al. 2020). Assuming a horizontally layered 1-D medium below each profile, we perform 1-D inversions of the phase dispersion curves. Since little baseline information is available for the Quito basin, we use a fairly general parametrization. Taking into account the geological descriptions presented in Jaya (2009) and Alvarado et al. (2014) for the Quito basin, we explore a simple three-layer model: two sedimentary layers and a volcanic bedrock. Love-wave phase velocities are inverted, using a layered model with uniform velocities, and shear-wave velocity ($V_s$) is set to increase with depth. The maximum half-space depth was fixed at 2000 m, and as a rule of thumb, the minimum resolved depth was set to a third of the shortest wavelength ($\lambda_{min}$), which was estimated as $\lambda_{min} = V_{min}/f_{max}$ (Wathelet et al. 2004). For Lines 1, 2 and 4, the minimum resolved thickness is around 215 m, while for Line 3, it is around 500 m. Using Geopsy, we sought the simplest $V_s$ profiles that satisfy the phase dispersion curves picked in the dispersion images in Fig. 12.

Computed models are represented according to their misfit with a grey colour-scale in the lower panel of Fig. 12. We computed 5030 models for each manually picked dispersion curve. The associated

theoretical Love-wave phase velocity dispersion curve of the best profile was plotted with a dashed red line in the dispersion images panel (Fig. 12) to show the fit with the picked points.

Taking advantage of the fact that the HVSR curves can give us much more information in addition to $f_0$ and that in several works, they have been successfully inverted to generate velocity profiles (e.g. Garcia-Jerez et al. 2016; Piña Flores et al. 2021); at a later stage, we invert the HVSR curves computed through DFA (M1 method in Fig. 4) for the 11-station profile, together with the Love-wave dispersion curves previously obtained on lines 1, 2 and 3 (Fig. 12) to generate 1-D shear-wave velocity profiles. In the joint inversion, we use HVSR software (Garcia-Jerez et al. 2016) and we process together HVSR and dispersion curves to reduce the non-uniqueness problem. HVSRs carry information associated with velocity structure, such as strong $V_s$ contrast. Dispersion curves, on another side, are sensitive to the absolute velocity variation with depth due to their frequency dependence. In our case, HVSRs are related to punctual structures along an essentially vertical path, whereas dispersion curves are expressed across segments.

Attending to the stratigraphic description presented in Jaya (2009) and Alvarado et al. (2014), where around three geological
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Simplified shear-velocity model obtained. The symbols used in the figure are detailed in the symbology section.</figcaption></figure>

units are distinguished, the inversion procedure began with establishing parameters ranges for general models of two layers over a half-space. Seismic velocities, density, thickness and Poisson's ratio of each layer obtained from the previous dispersion curves inversion (Fig. 12) were used as initial models. We carried out the joint inversion assigning equal weights to HVSR and dispersion curves. Unlike other seismic stations, since three peaks are visible on the HVSRs of JBPC and HLUZ, a general model of three layers over a half-space has been considered. In the joint inversion of the HVSRs of the CRON, HCAL and SOCA, we use the dispersion curve corresponding to Line 1 (Fig. 12), with the HVSRs of JBPC, OBSA and QUEM, we use the dispersion curve corresponding to Line 2 (Fig. 12), while with the HVSRs of ARGE, QUIB, HLUZ and GARO, we use the dispersion curve of Line 3 (Fig. 12). We decided not to do a joint inversion using the INCA HVSR due to the anomalous $f_0$ amplitude value (close to 10, Fig. 7a). We have no certainty that this value is not influenced by soil saturation or a tilt in the seismometer base, among other things.

Figs 13(a) and (b) show the computed models according to the ratio misfit/min(misfit) calculated for both the dispersion curve and the HVSR. We carry out more than 5000 models for each dispersion curve–HVSR pair, but we only plot those models with a ratio less than 7. Measured curves are represented by green lines, while the red dashed lines highlight the best models.

In $V_s$ models (Fig. 13c), the difference in the half-space depth is significant. On CRON, HCAL and SOCA, this limit is around 200 m deep. On OBSA and QUEM, it is observed slightly deeper, around 250 m. However, the half-space depth on ARGE, QUIB and HLUZ is greater than 700 m. The half-space depth in JBPC is around 500 m; however, there is no good agreement between the best models and the measured curves. On GARO, the depth of the half-space is around 100 m, but in this station, we have little confidence in the HVSR low-frequency part due to the high level of noise (Figs 7a and 13b). The lack of low frequencies resolution in the HVSR curve could explain why we cannot observe a contrast at a depth such as for HLUZ or QUIB.

The $V_s$ of the half-space in the south stations is higher (around 2500 ms⁻¹) than the observed in the north and centre stations (around 1700 ms⁻¹).

## 7 DISCUSSION

The background seismic noise in the stations of 4B-array is characterized by the presence of coherent ocean-generated seismic noise that dominates at low frequencies (less than 1 Hz) and coherent cultural noise at high frequencies (greater than 1 Hz) as it has been also seen by Díaz *et al.* (2017) in other ambient noise measurements in urban zones. From the PDFs representative spectra, we observe a large variability in the ambient noise level along the Quito basin. In the frequency range between 0.1 and 0.5 Hz, stations located in the south-western part (ARGE, QUIB and HLUZ) present systematically higher noise levels (Figs 2 and 4). This may reflect the presence of different infilling sediments as it was recently observed by Smith & Tape (2019) in the Nenana basin (Alaska), where they report that underlying sediments strongly influence the seismic wavefield amplitude levels.

When calculating HVSR and mapping $f_0$ through the Quito basin, it is observed that the value of $f_0$ decreases towards the south-west, which can be interpreted as a deepening of the basin or a change in the superficial velocity of the sediments. The zone comprising ARGE, QUIB and HLUZ stations presents the lowest values of $f_0$, around 0.3 Hz (Fig. 7). This suggests that this zone shows the greatest sediment thicknesses or the presence of low velocities. To the northeast, CRON and HCAL stations with $f_0$ around 1.2, would probably be marking the basin's eastern limit. The value of $f_0$ in the east-west line formed by the PARU-JBPC-INAQ-CSEK stations decreases towards the east (Fig. 7). It could then be interpreted as a deepening of the basin toward the east. This is in accordance with the QFS trace and its overall behaviour reported by Alvarado *et al.* (2014) (Figs 1 and 7).

We also investigated the effect of azimuth on HVSRs by rotating horizontal components through a set of azimuths (Fig. 6). By not
finding significant variability in both $f_0$ and its amplitude, associated with different azimuths, we can assume that the ambient noise wavefield is azimuthally isotropic and that the sites, in their great majority, are 1-D. Therefore, we can interpret that the subsurface in Quito is horizontally stratified with homogeneous layers that extend for great distances. While this 1-D assumption seems to work well at most sites, it is likely to be less convincing in the centre of the basin, at JBPC, given the HVSR curve at this location is complex, without any clear peaks, and with amplitudes larger than 2 for an extended range of frequencies. Morphologically in this area, the Quito basin narrows, measuring about 3.5 km in the EW direction.

Cross-correlation of more than 1 yr of ambient seismic noise was computed to reconstruct Rayleigh and Love waves traveling between pairs of stations. Coherent Love waves are much more energetic and clear than Rayleigh waves (Figs 9, 10 and 11) and exhibit significant energy in the frequency range 0.1–2 Hz. In the Kanto basin (with an area about 17 000 km², Koketsu & Kikuchi 2000), where Rayleigh waves are predominant, the shear-wave velocity inside the basin gradually increases with depth (e.g. Takemura et al. 2015); however, in the smaller Osaka basin (of elliptical shape and with an area about 3600 km², Asano et al. 2016), where Love waves predominate over Rayleigh waves, a clear unconformity at the sediment-bedrock boundary (between 1 and 1.5 km deep) with a strong impedance contrast in shear-wave velocity is observed at the subsurface (Asano et al. 2017). The free surface condition is a sufficient condition for Rayleigh waves to develop, even in homogeneous (or smooth varying) media. In contrast, Love waves are subjected to the presence of a strong velocity contrast in depth. This suggests that the clear predominance of Love waves in the Quito basin would indicate the preference of a layered velocity model presenting a strong impedance contrast in depth, similar to the case of Osaka basin, maybe the best approximation for the Quito subsurface.

When analysing the cross-correlation functions, it is observed that signals between 0.5 and 2 Hz fade to the basin's southwest (Fig. 10). Viens et al. (2016) show the low quality of Green's functions recovered from cross-correlations between pairs of stations in a profile of seismic stations in the Kanto basin can be explained by important changes in the basin's geometry. In the case of Viens et al. (2016), when a bay separated the virtual source from the receiver stations, new scattered wave sources are found that violate the assumption of distant ambient noise sources. In Quito's case, the seismic waves undergo a kind of low-pass filter before entering the southwestern zone (Figs 10 and 11). This behaviour could be associated with an important change in the basin's structure between the north, central and the southern parts.

Based on the results obtained in this work, we can argue that: (1) the measurements of $f_0$ in ARGE, QUIB and HLUZ point to the interface bedrock-sediments (Figs 13 and 14) found to the south of the basin through the joint inversion of Love-wave dispersion and HVSR curves. The frequency signals close to 0.3 Hz observed in the cross-correlation functions in the three components throughout the basin, along with low-frequency $f_0$ picks found at sites like CSEK ($f_0$ = 0.4) and MAST ($f_0$ = 0.4) could indicate that this bedrock is present all along the basin. (2) In the north-central part, the bedrock seems to be overlaid by a geological stratum (yellowish layers in Fig. 14) not found in the basin's southern part. This layer could correspond to a volcanic-origin formation identified by Jaya (2009) based on fieldwork and drilling records. According to the authors this formation is not found to the south of the Quito's basin, and it is only present in the centre and north with an identified maximum thickness of 200 m. A gradual transition of $V_s$ in the north and central

parts (Fig. 14) could cause the joint inversion to be insensitive to the bedrock interface. The limit between shallow sediments (bluish layers in Fig. 14) and the volcanic-origin formation (yellowish layers in Fig. 14) could be the interface that excites the HVSR peak at around 0.6–0.8 Hz observed in the north-central basin's part. The boundary between surface sediments (bluish layers in Fig. 14) and intermediate sediments (greenish layers in Fig. 14) in the south could be the interface that excites a second/third peak in the HVSR curves at around 3 Hz.

## 8 CONCLUSIONS

We present an original study based on the analysis of almost 2 yr of continuous noise recordings of the 4B-array in the Quito basin. The goal of our experiment is to obtain a first model of the deep structure that could explain the seismic amplifications observed on earthquakes recordings (Laurendeau et al. 2017). For this, we use both cross-correlation functions and HVSR curves.

Using four virtual source gathers and inverting Love waves phase velocity dispersion curves jointly with HVSR curves, we find that the sedimentary column in the southern part of the basin is much thicker (around 900 m) than in the north-central part of the basin, where a strong seismic reflector is found at a depth of around 200–300 m. This result is consistent with the rather low fundamental frequency value ($f_0$ = 0.28–0.4 Hz) obtained using HVSR in the southern part of the basin and the higher value in the centre-north part. Based on geological studies and the shallow $V_s$ values found from the inversion, we favour the hypothesis of differential infilling of the basin. Indeed, a large volcanic deposit exists in the north-central part of the Quito basin that is not present in the southern part. This thick volcanic layer that covers the central and northern part of the city may play the role of a natural seismic isolator that protects this part of the city from the strong low-frequency amplification that exists in the south.

The results presented in this paper constitute the first velocity model down to 1 km depth of the Quito basin, a highly populated urban area exposed to high seismic risk. It is a first step towards the detailed structure of the Quito basin we wish to obtain in the future.

## ACKNOWLEDGEMENTS

We are thankful to two anonymous reviewers and the editor, Victor Cruz-Atienza, for their insightful and constructive comments, which helped us improve the manuscript. We also want to thank Christophe Maron for his help in constructing the instrumental responses and Anne Deschamps for her help in the seismic database creation. We thank the IG-EPN staff for the effort to maintain the 4B network. The authors gratefully acknowledge the comments of Guust Nolet and Allie Hutchison for improving the quality of this paper. This work has been founded by the IRD and IG-EPN in the frame of the International Joint Laboratory *Earthquakes and Volcanoes in the Northern Andes* and *Remake* project ANR-15-CE04-0004.

## DATA AVAILABILITY

The 4B-array dataset will be made publicly available. It will either be hosted online or freely sent on external hard disks upon request.
REFERENCES

Albarello, D. & Lunedei, E., 2013. Combining horizontal ambient vibration components for H/V spectral ratio estimates, *Geophys. J. Int.*, **194**(2), 936–951.

Alvarado, A. et al., 2014. Active tectonics in Quito, Ecuador, assessed by geomorphological studies, GPS data, and crustal seismicity, *Tectonics*, **33**(2), 67–83.

Alvarado, A. et al., 2018. Seismic, volcanic, and geodetic networks in ecuador: building capacity for monitoring and research, *Seismol. Res. Lett.*, **89**(2A), 432–439.

Asano, K., Sekiguchi, H., Iwata, T., Yoshimi, M., Hayashida, T., Saomoto, H. & Horikawa, H., 2016. Modelling of wave propagation and attenuation in the Osaka sedimentary basin, western Japan, during the 2013 Awaji Island earthquake, *Geophys. J. Int.*, **204**(3), 1678–1694.

Asano, K., Iwata, T., Sekiguchi, H., Somei, K., Miyakoshi, K., Aoi, S. & Kunugi, T., 2017. Surface wave group velocity in the Osaka sedimentary basin, Japan, estimated using ambient noise cross-correlation functions, *Earth, Planets Space*, **69**(1), 108.

Bahavar, M., Spica, Z.J., Sánchez-Sesma, F.J., Trabant, C., Zandieh, A. & Toro, G., 2020. Horizontal-to-vertical spectral ratio (HVSR) IRIS station toolbox, *Seismol. Res. Lett.*, **91**(6), 3539–3549.

Beauval, C., Yepes, H., Bakun, W.H., Egred, J., Alvarado, A. & Singaucho, J.-C., 2010. Locations and magnitudes of historical earthquakes in the Sierra de Ecuador (1587-1996), *Geophys. J. Int.*, **181**(3), 1613–1633.

Beauval, C., Yepes, H., Audin, L., Alvarado, A., Nocquet, J.M., Monelli, D. & Danciu, L., 2014. Probabilistic seismic-hazard assessment in Quito, estimates and uncertainties, *Seismol. Res. Lett.* **85**(6), 1316–1327.

Behm, M., Leahy, G.M. & Snieder, R., 2013. Retrieval of local surface wave velocities from traffic noise—an example from the La Barge basin (Wyoming), *Geophys. Prospect.*, **62**(2), 223–243.

Bensen, G., Ritzwoller, M., Barmin, M., Levshin, A., Lin, F., Moschetti, M., Shapiro, N. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys. J. Int.*, **169**(3), 1239–1260.

Bindi, D. et al., 2009. Site amplifications observed in the Gubbio Basin, Central Italy: hints for lateral propagation effects, *Bull. seism. Soc. Am.*, **99**(2A), 741–760.

Boese, C., Wotherspoon, L., Alvarez, M. & Malin, P., 2015. Analysis of anthropogenic and natural noise from multilevel borehole seismometers in an urban environment, Auckland, New Zealand, *Bull. seism. Soc. Am.*, **105**(1), 285–299.

Castellaro, S., 2016. The complementarity of H/V and dispersion curves, *Geophysics*, **81**(6), T323–T338.

Díaz, J., Ruiz, M., Sánchez-Pastor, P.S. & Romero, P., 2017. Urban seismology: on the origin of earth vibrations within a city, *Scient. Rep.*, **7**(1), doi:10.1038/s41598-017-15499-y.

García-Jerez, A., Piña-Flores, J., Sánchez-Sesma, F.J., Luzón, F. & Perton, M., 2016. A computer code for forward calculation and inversion of the H/V spectral ratio under the diffuse field assumption, *Comput. Geosci.*, **97**, 67–78.

Guéguen, P., Chatelain, L.J., Guillier, B. & Yepes, H., 2000. An indication of the soil topmost layer response in Quito (Ecuador) using noise H/V spectral ratio, *Soil Dyn. Earthq. Eng.*, **19**(2), 127–133.

Hillers, G., Ben-Zion, Y., Lanès, M. & Campillo, M., 2013. Interaction of microseisms with crustal heterogeneity: a case study from the San Jacinto fault zone area, *Geochem. Geophys. Geosyst.*, **14**(7), 2182–2197.

Hutko, A.R., Bahavar, M., Trabant, C., Weekly, R.T., Fossen, M. & Ahern, T., 2017. Data products at the IRIS-DMC: growth and usage, *Seismol. Res. Lett.*, **88**(3), 892–903.

Jaya, D., 2009. Origen de los depósitos volcánicos de la cuenca de Quito: Implicación de la evolución del complejo Pichincha y del sistema de fallas de Quito, *Master 2 SGT PREFALC*, Université Sophie Antipolis Nice (in Spanish).

Kästle, E.D., El-Sharkawy, A., Boschi, L., Meier, T., Rosenberg, C., Bellahsen, N., Cristiano, L. & Weidle, C., 2018. Surface wave tomography of the alps using ambient-noise and earthquake phase velocity measurements, *J. geophys. Res.*, **123**(2), 1770–1792.

Koketsu, K. & Kikuchi, M., 2000. Propagation of seismic ground motion in the Kanto Basin, Japan, *Science*, **288**(5469), 1237–1239.

Landès, M., Hubans, F., Shapiro, N.M., Paul, A. & Campillo, M., 2010. Origin of deep ocean microseisms by using teleseismic body waves, *J. geophys. Res.*, **115**(B5), doi:10.1029/2009JB006918.

Larose, E., Derode, A., Campillo, M. & Fink, M., 2004. Imaging from one-bit correlations of wideband diffuse wave fields, *J. appl. Phys.*, **95**(12), 8393–8399.

Laurendeau, A. et al., 2017. Low-frequency seismic amplification in the Quito Basin (Ecuador) revealed by accelerometric recordings of the RENAC network, *Bull. seism. Soc. Am.*, **107**(6), 2917–2926.

Lin, F., Moschetti, M.P. & Ritzwoller, M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, *Geophys. J. Int.*, **173**(1), 281–298.

Ma, Y. & Clayton, R.W., 2016. Structure of the Los Angeles Basin from ambient noise and receiver functions, *Geophys. J. Int.*, **206**(3), 1645–1651.

Manea, E., Michel, C., Poggi, V., Fäh, D., Radulian, M. & Balan, F., 2016. Improving the shear wave velocity structure beneath Bucharest (Romania) using ambient vibrations, *Geophys. J. Int.*, **207**(2), 848–861.

Mariniere, J., Nocquet, J.-M., Beauval, C., Champenois, J., Audin, L., Alvarado, A., Baize, S. & Socquet, A., 2019. Geodetic evidence for shallow creep along the Quito fault, Ecuador, *Geophys. J. Int.*, 2039–2055.

McNamara, D.E. & Buland, R.P., 2004. Ambient noise levels in the Continental United States, *Bull. seism. Soc. Am.*, **94**(4), 1517–1527.

McNamara, D.E., Stephenson, W.J., Odum, J.K., Williams, R. & Gee, L., 2015. Site response in the eastern United States: a comparison of Vs30 measurements with estimates from horizontal:vertical spectral ratios, *USGS Publ. Warehouse*, **509**, 67–79.

Mercerat, D., Laurendeau, A., Courboulex, F., Alvarado, A., Bonilla, L. F., Bertrand, E., Gueguén, P., Charvis, P., Langlaude, P., Perrault, M., Barros, J. G. & Reyes, M. F., 2016. ANR-remake Quito basin (Remake Quito), International Federation of Digital Seismograph Networks. Dataset/Seismic Network, doi:10.7914/SN/4B_2016.

Molnar, S., Cassidy, J., Castellaro, S., Cornou, C., Crow, H., Hunter, J., Matsushima, S., Sánchez-Sesma, F. & Yong, A., 2018. Application of microtremor horizontal-to-vertical spectral ratio (MHVSR) analysis for site characterization: state of the art, *Surv. Geophys.*, **39**(4), 613–631.

Nakamura, Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, *Quater. Rep. Railway Tech. Res. Inst.*, **30**(1), 25–33.

Nakamura, Y., 2019. What is the nakamura method?, *Seismol. Res. Lett.*, **90**(4), 1437–1443.

Nakata, N., Gualtieri, L. & Fichtner, A., 2019. *Seismic Ambient Noise*, Cambridge Univ. Press.

Park, C., 2011. Imaging dispersion of MASW data—full vs. selective offset scheme, *J. Environ. Eng. Geophys.*, **16**(1), 13–23.

Park, C.B., Miller, R.D. & Xia, J., 1999. Multichannel analysis of surface waves, *Geophysics*, **64**(3), 800–808.

Pastén, C., Sáez, M., Ruiz, S., Leyton, F., Salomón, J. & Poli, P., 2016. Deep characterization of the Santiago Basin using HVSR and cross-correlation of ambient seismic noise, *Eng. Geol.*, **201**, 57–66.

Peterson, J.R., 1993. Observations and modeling of seismic background noise, Tech. rep.

Piña Flores, J., Cárdenas-Soto, M., Sarabia-González, A., García-Jerez, A., Sierra-Alvarez, A., C., Sáenz-Castillo, A., M., Luzón, F. & Sánchez-Sesma, J.F., 2021. Imaging the structure of the Sun Pyramid (Teotihuacán, Mexico) from passive seismic methods, *Eng. Geol.*, **281**, 105969, doi:10.1016/j.enggeo.2020.105969.

Poli, P., Pedersen, H. & Campillo, M., 2011. Emergence of body waves from cross-correlation of short period seismic noise, *Geophys. J. Int.*, **188**(2), 549–558.

Ringler, A.T. & Hutt, C.R., 2010. Self-noise models of seismic instruments, *Seismol. Res. Lett.*, **81**(6), 972–983.

Sánchez-Sesma, F. J. et al., 2011. A theory for microtremor H/V spectral ratio: application for a layered medium, *Geophys. J. Int.*, **186**(1), 221–225.
Saygin, E. et al., 2015. Imaging architecture of the Jakarta Basin, Indonesia with transdimensional inversion of seismic noise, *Geophys. J. Int.*, **204**(2), 918–931.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**(5715), 1615–1618.

Smith, K. & Tape, C., 2019. Seismic noise in Central Alaska and influences from rivers, wind, and sedimentary basins, *J. geophys. Res.*, **124**(11), 11 678–11 704.

Takemura, S., Akatsu, M., Masuda, K., Kajikawa, K. & Yoshimoto, K., 2015. Long-period ground motions in a laterally inhomogeneous large sedimentary basin: observations and model simulations of long-period surface waves in the northern Kanto Basin, Japan, *Earth, Planets Space*, **67**(1), 33, doi:10.1186/s40623-015-0201-7.

Taylor, G., Rost, S., Houseman, G.A. & Hillers, G., 2019. Near-surface structure of the North Anatolian Fault zone from Rayleigh and Love wave tomography using ambient seismic noise, *Solid Earth*, **10**(2), 363–378.

Vaca, S., Vallée, M., Nocquet, J.-M. & Alvarado, A., 2019. Active deformation in Ecuador enlightened by a new waveform-based catalog of earthquake focal mechanisms, *J. South Am. Earth Sci.*, **93**, 449–461.

Vassallo, M., De Metteis, R., Bobbio, A., Di Giulio, G., Adinofi, G., Cantore, L., Cogliano, R., Fodarella, A., Maresca, R., Pucillo, S. & Riccio, G., 2019. Seismic noise cross-correlation in the urban area of Benevento city (Southern Italy), *Geophys. J. Int.*, **217**(3), 1524–1542.

Ventosa, S., Schimmel, M. & Stutzmann, E., 2017. Extracting surface waves, hum and normal modes: time-scale phase-weighted stack and beyond, *Geophys. J. Int.*, **211**(1), 30–44.

Viens, L., Koketsu, K., Miyake, H., Sakai, S. & Nakagawa, S., 2016. Basin-scale Green's functions from the ambient seismic field recorded by MeSO-nett stations, *J. geophys. Res.*, **121**(4), 2507–2520.

Villaseñor, A., Yang, Y., Ritzwoller, M.H. & Gallart, J., 2007. Ambient noise surface wave tomography of the Iberian Peninsula: implications for shallow seismic structure, *Geophys. Res. Lett.*, **34**(11), doi:10.1029/2007GL030164.

Wapenaar, K., 2004. Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation, *Phys. Rev. Lett.*, **93**(25), 254301, doi:10.1103/PhysRevLett.93.254301.

Wapenaar, K., Draganov, D., Snieder, R., Campman, X. & Verdel, A., 2010. Tutorial on seismic interferometry: Part 1—basic principles and applications, *Geophysics*, **75**(5), 75A 195–75A 209.

Wathelet, M., 2008. An improved neighborhood algorithm: parameter conditions and dynamic scaling, *Geophys. Res. Lett.*, **35**(9), doi:10.1029/2008GL033256.

Wathelet, M., Jongmans, D. & Ohrnberger, M., 2004. Surface-wave inversion using a direct search algorithm and its application to ambient vibration measurements, *Near Surf. Geophys.*, **2**(4), 211–221.

Wathelet, M., Chatelain, J., Cornou, C., Giulio, G., Guillier, B., Ohrnberger, M. & Savvaidis, A., 2020. Geopsy: a user-friendly open-source tool set for ambient vibration processing, *Seismol. Res. Lett.*, **91**(3), doi:10.1785/0220190360.

Zheng, S., Sun, X., Song, X., Yang, Y. & Ritzwoller, M.H., 2008. Surface wave tomography of China from ambient seismic noise correlation, *Geochem. Geophys. Geosyst.*, **9**(5), doi:10.1029/2008GC001981.