Evidence for HV peaks superposition leading to extreme horizontal ground motion amplification revealed by the Xochimilco ambient noise tomography

SUMMARY

Mexico City is one of the largest cities in North America, facing high seismic hazards and water supply problems. This paper presents an ambient seismic noise tomography of the city's south area in Xochimilco, where large amplifications have already been registered during subduction earthquakes. Eighty-four seismic stations have been installed, and their records processed. The tomography method combines the inversion of the horizontal-to-vertical spectral ratios (HVSR) and multimodal dispersion curves. The importance of considering a multimodal approach is justified in light of the complex geological setting. The dispersion curves analysis shows that the surface wave energy is divided over the fundamental and the higher modes, particularly between 50 and 300 m s⁻¹, and in the whole frequency range analysed. We observe a spatially continuous decrease of the dominant peak frequency of the HVSRs toward the lake interior but a heterogeneous amplification. By analysing the velocity profiles associated with the highest amplifications, we discovered that these latter result from the superposition of several resonance peaks. Their coincidence in frequency is due to the overall constant linear gradient velocity in the sedimentary basin crossed by several low-velocity anomalies due to high water content or high-velocity anomalies due to lavas. Although most of the shallow water is trapped in clay sediment, the velocity model also allows for identifying deeper water reserves. All these analyses are of fundamental importance for the correct seismic mitigation in Mexico City but might also be extended to other cities built on top of sedimentary basins.

1. INTRODUCTION

About 22 million people lived in 2020 in the Metropolitan Area of the Valley of Mexico (MAVM), that is, approximately 20 per cent of Mexico's population and 17 per cent of the nation's Gross Domestic Product (INEGI 2020). The Metropolitan area corresponds mainly to the Mexican Basin, which is well known for its earthquake amplification effects (Lermo & Chávez-García 1994; GCDMX 2004; Mayoral et al. 2019) and for its water reserves (Ortega‐Guerrero et al. 1993; González-Morán et al. 1999; Ruiz & Ruiz 2013; Martinez et al. 2015; Neukirch et al. 2021; Palma Nava et al. 2022). Although the MAVM is located 340 km away from the subduction trench and is equipped with a seismic alert (Espinosa Aranda et al. 1995; Iglesias et al. 2007; Suárez 2022), frequent large earthquakes, due to the convergence of the Cocos, Rivera and North American Plates (Pardo & Suarez 1995; Ferrari et al. 2012; González-Flores et al. 2020), have impinged the MAVM with substantial damage (Singh et al. 1985, 1996, 1998, 1999; Lermo et al. 1988; Singh, Mena et al. 1988; Campillo et al. 1989). Additionally, the MAVM belongs to the central southern part of the Trans-Mexican Volcanic Belt (TMVB) (Fig. 1), and several regional and local active faults increase the seismic hazard (Iglesias et al. 2002; García et al. 2004; Ovando-Shelley et al. 2012; Lermo et al. 2016; Singh et al. 2018). The seismic risk management in the MAVM relies on amplification maps that historically have been developed to characterize the site effects and update construction codes (Singh, Lermo et al. 1988; Ordaz & Reyes 1999; Reinoso & Ordaz 2001; GCDMX 2004; Cardona et al. 2010; Auvinet & Juárez 2011; González et al. 2011; Ordaz et al. 2013; Clemente-Chavez et al. 2014). According to the normative, the subsoil of the MVMA is classified into three main geotechnical zones (GCDMX 2004): I, II and III, usually referred respectively to as hills zone, transition zone and lake zone as they are in good agreement with the thickness or absence of the soft sediments from the ancient Lake Texcoco (see Fig. 1). Note that the denomination ‘lake’ refers in this paper to the ancient lake and that at present time, the lakes have all almost entirely disappeared. Zone III amplifies the acceleration not only in intensity but also in duration due to regional effects (Iida 1999; Roullé & Chávez-García 2006) and to the resonance of higher surface wave modes inside the basin (Rivet et al. 2015; Cruz-Atienza et al. 2016a). However, current amplification maps are produced by interpolating scarce measurements and are imprecise in fully representing the subsoil's complexity (Clemente-Chavez et al. 2014). Increasing the measurement points in a city is challenging due to the existing infrastructure, installation permits or security issues, and not always efficient due to the high level of heterogeneities. A detailed seismic model of the Mexican Basin is an alternative to predict the seismic response and amplification levels throughout the city from a deterministic approach. The model can also help better understand the origins of such exceptional seismic amplification effects.

The Mexican Basin has been the subject of long-term research using geological, geophysical and borehole data (Pérez-Cruz 1988; Urrutia-Fucugauchi & Flores-Ruiz 1996; Campos-Enríquez et al. 1997; Arce et al. 2019; González-Flores et al. 2020a). Information from the National Seismological Network, the Mexico City Seismic Network array and some temporary experiments have been used to develop some regional seismic models of the area (Campillo et al. 1996; Iglesias et al. 2001, 2010; Chávez-García & Quintanar 2010; Cruz-Atienza et al. 2010; Spica et al. 2016; Rodríguez-Domínguez et al. 2019; Rosalia et al. 2020). Conceptual models based on geological and geotechnical observations have also been advanced (Vazquez-Sanchez & Jaimes-Palomera 1989; Arce et al. 2013; Auvinet et al. 2017; Mooser 2018; Arce et al. 2019). The main subsoil characteristics consist of a bed of saturated clays (usually referred to as ‘soft sediments’) overlying a lacustrine formation interbedded with pyroclastic deposits. It results in several high shear wave velocities (⁠|${V_s}$|⁠) contrasts and a globally strong gradient in the first tens of metres. Incorporating the geological information together with regional inversions has also allowed the development of more geophysical models to explain the seismic behaviour of the basin (Ramirez-Guzman et al. 2007; Cruz-Atienza et al. 2016a; Aguilar-Velázquez et al. 2023). However, the interstation distances did not allow a vertical resolution thinner than 2 km. Higher resolution in the surface layers is essential because the main seismic peak amplification is related to the first tenth of metres of the shallowest structure (e.g. vs30; Borcherdt 1994; Boore et al. 1997; Boore 2004) and two large water reservoirs have been found within the first 1500 m depth (Chouteau et al. 1994).

This paper uses ambient seismic field interferometry (ASFI) to improve vertical resolution. ASFI is a geophysical technique widely used to assess the subsurface elastic properties. It is ideal in urban environments as it avoids the usage of active sources such as explosions or strong vibrations. Two ASFI techniques stand out: The horizontal-to-vertical-spectral-ratio (HVSR) and the ambient noise tomography (ANT). They have already been used to assess the shallow velocity structure in local areas inside Mexico City (Rosas et al. 2011; Cárdenas-Soto 2016; Vergara-Huerta & Aguirre 2020; Cárdenas-Soto et al. 2021).

This study presents an ANT of the Xochimilco area (Fig. 1). Xochimilco is well-known for its tourist canals and is one of the few areas where natural lakes still exist. Due to the lake sediments, the Xochimilco and the nearby Chalco areas present a thick bed of saturated clays responsible for some of the highest amplification in the city (Juárez-Camarena et al. 2016a). Multiple volcanic layers from the Chichinautzin volcanic range interbed the lacustrine deposits (Arce, Layer, Lassiter et al. 2013). We first discuss the difficulty of making an ANT in such a complex setting based on a known local velocity profile. Subsequently, we detail the method: a tomography based on a joint HVSR and multimodal dispersion curves (DC) inversion. Finally, we interpret and evaluate the results in terms of geological knowledge.

2. GEOLOGY OF THE XOCHIMILCO AREA

The Mexican Basin is located in the central southern part of the TMVB. Its formation responds to the regional tectonic and volcanic activities. In the past, it hosted five major lacustrine bodies that provided a continuous accumulation of clays and limes (Lozano-García & Ortega-Guerrero 1998; Siebe et al. 2004; Siebe & Macías 2006; Arce et al. 2013; González Torres et al. 2015). The continuous expansion of the MVMA led to the desiccation of the main water bodies, and most existing rivers were channelled to drain the rainwater out of the city and prevent inundations. As a result, an important part of the shallow subsoil consists of saturated clays, usually called ‘soft sediments,’ with thicknesses ranging from 5 to 100 m (GCDMX 2004; Juárez-Camarena et al. 2016b). These soft sediments significantly contribute to the amplification and subsidence problems in Mexico City. The soft sediments overlay a two kilometres-thick sequence of Oligocene to Miocene (Quaternary) lacustrine deposits interbedded with ‘hard layers’ of volcanic deposits such as sands, tufts, gravels and lava flows (Arce et al. 2019). This particular intercalation of lacustrine deposits and volcanic rocks is the origin of the storage of one of the main aquifers of Mexico City (Ortega‐Guerrero et al. 1993; Campos-Enríquez et al. 1997a; Siebe et al. 2004). Below this sequence stands Cretaceous limestones of marine origin (Singh et al. 1995; Campos-Enríquez et al. 1997a; Arce et al. 2019; González-Flores et al. 2020).

Our study stands in the Xochimilco area, within the Chalco-Xochimilco sub-Basin (Fig. 1c). The sub-Basin is a natural recharge area that contains two large remnant lakes, the Xochimilco Lake and the Chalco Lake (Ortiz-Zamora & Guerrero 2007). Gravity and seismic data with geological and borehole information (Pérez-Cruz 1988; Campos-Enríquez et al. 1997a; González-Flores et al. 2020a) allow the defining of four main units: the first soft sediment layer containing the aquitard, two lacustrine-volcanic units and the deeper limestone formations.

3. DATA DESCRIPTION

A temporary seismic network consisting of 18 triaxial seismometers was deployed in the Xochimilco area between June and September 2017. The seismometers are Guralp middle band (0.03–100 Hz) CMG-6TD models with integrated digitizer recording at 100 samples per second (sps). The stations were oriented with the magnetic pole and buried at 0.6 m when possible. The seismometers recorded data for three to four days and then moved, resulting in six temporary networks with 84 measurement points (Fig. 1c).

4. TOMOGRAPHY PROCEDURE: JOINT INVERSION OF DCs AND HVSRs

We conducted a shear wave velocity tomography using a joint inversion of multimodal surface wave DCs and HVSRs obtained from ambient noise processing. The basis of the processing relies on the fact that the cross-correlation functions (CCFs), using seismic noise records, are proportional to the imaginary part of the Green Function (GF) (Lobkis & Weaver 2001; Weaver 2005; Sánchez-Sesma & Campillo 2006). The CCFs are subsequently used to extract the DCs between two distinct measurement points or the HVSRs in the particular case of the autocorrelation at a single station (Sánchez-Sesma et al. 2011).

The classical approach for tomography based on DCs considers only the fundamental mode, assuming it is the most energetic and easily identified by following the line of higher energy in the velocity-frequency domain. However, this approach has three significant disadvantages. (i) It smooths strongly the resulting velocity model (M. Perton et al. 2022). Although smoothing can mitigate minor artefacts and reduce reliance on the initial model, it also limits the ability to retrieve velocity inversions, which are often of primary interest. (ii) The fundamental mode identification may be biased. Dispersion curves may exhibit osculation points (i.e. points where two dispersion curves almost intersect) that significantly reduce the interpreter's capability to separate the fundamental mode from the higher modes. Furthermore, the higher modes may carry more energy than the fundamental mode for specific velocity profiles, as shown in Mexico City (Rivet et al. 2015; Cruz-Atienza et al. 2016b). Both situations can lead to incorrect mode identification and wrong tomography results. (iii) The largest interstation distance determines the maximal depth resolution of the tomographies and intrinsically results in a horizontally smooth inverted velocity model. On the contrary, a multimode procedure may improve lateral resolution and depth of investigation (Xia et al. 2003; Liu & Fan 2012; Spica et al. 2018; Wang et al. 2019; M. Perton et al. 2020; Ma et al. 2022).

Although most multimode tomographies still require a priori mode identification, we use a blind mode tomography procedure here to avoid the risk of mode misidentification. The inversion procedure tests several velocity models, looking for the one that maximizes the coincidences between the dispersion points and the corresponding theoretical dispersion curves. This blind multimode tomography method has successfully been applied in several geological contexts, such as volcanos, onshore natural gas fields and sedimentary basins (Spica et al. 2018; M. Perton et al. 2020, 2022). However, when the surface wave mode density is high, the inversion of the dispersion curve might not be unique (Viens et al. 2023) and should be guided by an additional observation, such as the HVSRs. As shown in these previous papers, the HVSRs are sensitive to high impedance contrasts, that is, relative velocity contrasts, while the DCs are sensitive to absolute phase velocities. Both techniques provide complementary information that helps reduce the non-uniqueness of the solution during the inversion procedure.

To conduct the tomography in an optimal configuration, we preliminarily present and discuss the theoretical DCs and HVSR associated with the Tulyehualco borehole velocity profile (COMESA 2014). These observations justify the importance of using a multimodal DC approach. We then present the accurate calculations of the experimental and theoretical HVSRs, which undergo numerous numerical instabilities for the given geological context and velocity contrasts. Subsequently, we discuss the determination of the surface wave dispersions and their joint inversion with HVSRs to obtain several 1-D velocity profiles. Finally, we explain how we build a 3-D velocity model from these 1-D velocity profiles.

4.1. Preliminary theoretical considerations from the Tulyehualco deep borehole

In Fig. 2(a), we present the velocity profile associated with the Tulyehualco deep borehole. The deeper part (>600 m) corresponds to velocity directly measured by the sonic suspension logging technique. However, this technique does not provide direct measurements in the shallow part because of the very soft and complex media that introduces mixed arrivals. The shallower part results from preliminary inversions realized by other seismologists (Erick Ramos Perez, private communication, 2022) and nearby geotechnical measurements (<100 m). The shear wave velocity varies from 0.05 to 2.5 km s⁻¹, that is, two orders of magnitude ratio between the minimum and maximum velocities. For the first 100 m, the shear wave (S-wave) velocity gradient is 11 s⁻¹. Until the depth of 450 m, the velocities depict a vertical gradient of about 3.2 s⁻¹ for S waves and 4.4 s⁻¹ for the longitudinal waves (P-wave). Since the lacustrine sediments are interbedded with hard layers of volcanic deposits, we also observe velocity inversions and low-velocity anomalies at 100 m depth and between 500 to 1150 m depth.

The effect of velocity gradient over the surface wave mode density has recently been discussed by Viens et al. (2023). However, the velocity profile is more complex here due to the velocity inversions. To evaluate the expected mode density and correctly calibrate our inversion scheme, we calculated in the first step the displacement kernels |$U( {k,f} )$| in the wavenumber-frequency domain for an unbounded layered media with flat and parallel interfaces and a source at the free surface. The calculation is done in 2-D in the |$\{ {x,z} \}$| plane by using the discrete wave number (DWN) method for the five components [P and vertically polarized S (SV) waves on the one hand and horizontally polarized S (SH) wave on the other hand] of the Green's function tensor (GF): |${G_{xx}}$|⁠, |${G_{xz}}$|⁠, |${G_{zx}}$|⁠, |${G_{zz}}$|⁠, |${G_{yy}}$|⁠. The first four components include the P, SV and Rayleigh waves. The last component corresponds to the SH and Love waves. Then, we express these kernels in the phase velocity-frequency domain |$( {U( {c,f} )} )$| by using the relation |$k = 2\pi f$|⁠. Finally, we take the absolute value and compress the amplitude of the images to help visualize all the modes (Figs 2c, d, e, f and g). As the surface waves correspond to the poles of the kernels, we introduced a small constant imaginary part in the frequency to avoid infinite amplitude. The surface waves appear in bright shades of grey, carrying the most energy. The shear bulk velocities correspond to horizontal asymptotes and are of weak energy. Only the fundamental surface waves show horizontal asymptotes with high energy. To help the mode identification, we superimposed the theoretical dispersion curves in |${G_{xx}}$|⁠, |${G_{yy}}$|⁠. The |${G_{xz}}$|⁠, |${G_{zx}}$|⁠, |${G_{zz}}$| have the same DCs as |${G_{xx}}$|⁠, but the distribution of the energy strongly differs between the different components. We can observe at least ten dispersion curves and several osculation points. The fundamental mode is the strongest only at low frequencies, below 0.7 Hz. The energy is mainly at low velocities in |${G_{zz}}$|⁠. On the other hand, the other Rayleigh components have a higher sensitivity to higher velocities and then to the deepest part of the velocity profile. The Love DCs show a uniform energy distribution and are consequently essential to avoid local minima during the inversion. Any method that only considers the fundamental DC mode or the vertical component of the seismic noise would have limited information about the subsoil structure.

Following Sánchez-Sesma et al. (2011), we also calculated the HVSR for the same velocity profile (Fig. 2b). We also present independently the H and V components. Given the frequency of the prominent peak, it is related to the fundamental Rayleigh and Love modes and reaches an amplification factor of 10. This high amplification is predominantly due to a low vertical energy density rather than a high horizontal energy. This observation is fundamental for correctly processing the HVs from the seismic signals, as discussed in the next section. The HVSR also presents several oscillations at high frequency and a low amplitude just after the prominent peak. As discussed in (García-Jerez et al. 2016; Piña-Flores et al. 2020, 2023), these characteristics are essential for correctly determining the near-surface subsoil.

4.2. HVSR technique

The HVSR is a classic technique to determine the dominant frequency at a specific location. It consists of the ratio of the horizontal over the vertical directional energies of the seismic signals (Nakamura 1989). Under the diffuse field assumption, Perton et al. (2009) demonstrated that the directional energy densities correspond to the imaginary part of the GFs, allowing the technique for subsoil properties assessment.

4.2.1. Observed HVSR

The energy density of the ambient noise seismic signal along the direction |$i = \{ {x,y,z} \}$| is obtained through the mean of its autocorrelation. In the frequency domain, the autocorrelation is the square absolute value of the signal, and the observed HVSR is the middle term of the following equation:

The spatial and frequential dependence |$( {x,\omega } )$| is only stated in the left term but should implicitly be considered for the |${u_i}$| and |${G_{ii}}$| terms. As the time derivation is equivalent to a multiplication by |$i\omega $|⁠, with i the complex number, accelerometer or velocimeter signals can be equally considered since these factors would cancel out in the division. The brackets ˂˃ represent the average over multiple time windows. The signals are detrended and filtered between 0.01 and 10 Hz. Since the spectral range almost covers two orders of magnitude, we divide it into four frequency bands to optimize the processing. The overlap between the bands is about 10 per cent and, in each band, we used a time window duration |$( {T = 6/{f_c}} )$| inversely proportional to the central frequency |$( {{f_c}} )$|⁠. We also applied a spectral whitening. The averaging process converges when the HV curves stop varying significantly for new windows. As seen in Fig. 2(b), the V component is expected to be very low at a specific frequency. Then, it is mandatory not to apply a smoothing filter, and the spectral whitening process is made adaptatively with a different width in each frequency band. A further advantage of adopting the diffuse field assumption lies in the initial averaging of each component, followed by the computation of the H/V ratio. Thus, smoothing procedures like the one proposed by Konno & Ohmachi (1998) to prevent V from approaching zero (which can occur when computing H/V solely from surface waves) are unnecessary. Indeed, the H and V components are energy densities of the entire field wave, which are always supposed to be strictly positive.

We only considered the signals at night in an urban environment to improve the calculation of the HVSRs. We finally checked their spatial coherence by rearranging them in a 4-D volume (latitude, longitude, frequency and amplitude, see Fig. 3). We plot the isosurfaces of the same amplitudes with unique colours. The low-amplitude isosurfaces are slightly transparent to allow the highest amplitudes to be observed. As expected, the HVSR frequencies evolve slowly between nearby stations because they stand above similar subsoil structures. The frequency of the dominant peak is about 0.2 Hz (i.e. 5 s) below the lake zone at the northeast and increases by almost a magnitude order towards the hill zone at the west (see Fig. 1). We can also observe a small secondary peak to the east at 3 Hz, which may indicate surface changes in the shallow subsurface. However, the amplitudes vary rapidly from one position to another. We discussed this unexpected behaviour later.

To better visualize the different shapes, maximum amplitudes and frequency range of the HVSR dominant peaks and the frequency shift according to the station positions, we present Fig. 4. All the green curves are from seismic stations belonging to the same N10 array and illustrate clearly the gradual frequency shift from 0.6 to 1.3 Hz and the amplification variation of the dominant peak. The black curves, associated with stations 132 and 138, only separated by 1.5 km, represent the highest (∼30) and lowest (∼10) amplification in this same N10 array. In the other arrays, the main peak frequencies might even be lower, but not necessarily with higher amplification.

4.2.2. Theoretical HVSR

The right term of eq. (1) is evaluated for an unbounded layered medium using the DWN (Sánchez-Sesma et al. 2011). This method requires the integration of the U (k, f) kernel between |$k = [ { - {k_{\rm max}},{k_{\rm max}}} ]$| to obtain |$U( {x,f} )$|⁠. The limit |${k_{\rm max}} = {\omega _{\rm max}}/{c_{\rm min}}$| depends on the smallest shear bulk wave |$( {{c_{\rm min}} = 20\, \rm m\,s^{-1}} )$| encountered in the velocity profile, which is a hundred times smaller than the highest velocity in the present case. We then require fine sampling along the k direction to consider all the wave contributions and for an unusually large |${k_{\rm max}}$|⁠. This results in evaluating the |$U( {k,f} )$| kernel along a large number of points, dramatically increasing the computation time.

4.3. Surface wave dispersion curves

The extraction of multimode dispersion points of the surface wave has already been presented extensively in Perton et al. (2020, 2022) and follows the next main steps. First, we compute the CCFs (Bensen et al. 2007) between all the concurrently operating stations. Then, we apply a frequency-time analysis (FTAN) (Dziewonski et al. 1969; Landisman et al. 1969) to identify the group velocity dispersion points between each station pair. Finally, we cluster the coherent information between station pairs with a straight path crossing the same area and a length close to one wavelength. These steps are detailed hereafter.

4.3.1. Cross-correlation

Since each of the six arrays worked at different times, we computed the cross-correlation functions for a total of 978 station pairs with interstation distances ranging from 0.2 to 5 km. For computational performance, we group the possible paths in five bins according to the interstation distances: [0–1], [1–2], [3–4], [4–5] km. The window duration for the cross-correlation was set equal to |$6{d_{\rm max}}/{c_{\rm min}}$|⁠, where |${d_{\rm max}}$| is the maximum distance of the bin. The computation of the cross-correlation functions is performed in the spectral domain. Each signal is detrended, re-sampled at 20 sps, bandpass filtered from 0.05 to 10 Hz, one-bit normalized and spectral whitened (Bensen et al. 2007). To improve the signal-to-noise ratio (SNR), we compute the cross-correlation with 75 per cent-time windows overlap and apply a temporal filter to remove spurious arrivals, that is, those that do not enter in the expected velocity arrivals [0.02 1.5] km s⁻¹. The effect of this time filter is shown in supplementary material in Fig. S1. Because several surface wave modes are expected, we do not observe a single pulse propagating through the network but rather several wave packets. We prefer not to discard CCFs at this point, that is, in the time domain, even if the SNR strongly decays with the interstation distance and does not allow us to observe arrivals distinctively above 2.5 km. We can indeed rescue some information by filtering the noisy signals appearing in a limited frequency band. To verify the convergence of the CCFs, we also compare them by considering different signal durations: two hours, one night and 4 or 5 nights (Fig. S2).

For each station pair, we first compute the full cross-correlation tensor [NN, NE, NZ, EN, EE, EZ, ZN, ZE, ZZ] components and rotate the tensor into the components [ZZ, ZR, ZT, RZ, RR, RT, TZ, TR, TT], where R is for radial and T for transverse to separate the contribution of the Rayleigh and Love waves. The rotationater Commission) with black dots and their pr angles are optimized to allow errors in the orientation of the seismometers to be corrected (M. Perton et al. 2022).

4.3.2. FTAN

The group arrival times of the Love (resp. Rayleigh) waves between each station pair are directly obtained in the FTAN of the TT (resp. ZZ, RR, RZ and ZR) components by selecting all local maxima at each frequency. These maxima are the dispersion points. To gather the information given by the four components of the Rayleigh waves, we do not consider the logarithmic stacking technique. In this common technique, the FTANs are stacked to improve the SNR (Campillo et al. 1996). As discussed in M. Perton et al. (2022) and observed in Fig. 2, the energy distribution of the dispersion curves of the four Rayleigh components is different, and stacking them may destroy the information. Instead, we merge the dispersion points of the four FTANs. At the end of this process, we save the dispersion points and do not gather them in dispersion curves because we cannot assign them to a specific mode.

4.3.3. Clustering of the dispersion points

As we retain all dispersion points without a selection criterion, the previous step may select artefacts that do not correspond to surface waves. To remove them, we identify the spatial coherence of the dispersion points. At each station position, we gather the dispersion points of all the individual station pairs with a path close to the station and a length of about one wavelength. We only conserve them if there are several occurrences (M. Perton et al. 2020, 2022), that is, several paths show similar dispersion points. Most artefacts were at low frequency and low velocity and associated with station pairs with large interstation distances and paths that span roads. The traffic along these roads introduces incoherent arrivals that are difficult to filter.

4.3.4. Forward calculation of DCs

As discussed in M. Perton et al. (2020, 2022), the dispersion points cannot be related, at this step, to their respective modes. As shown in Fig. 2, it is natural that the mode energy is not continuous along a dispersion curve. Worst, the osculation points in the plot of |${G_{yy}}$| make the energy continuous over different modes and not along the same mode. Then, we use a blind mode inversion to compare the dispersion points with theoretical dispersion curves. These latter are obtained assuming an unbounded horizontally layered structure below each of the seismometers and using the scheme presented in Perton & Sánchez-Sesma (2016).

4.4. Joint inversion

As shown for the Tulyehualco case, the surface wave mode density is high, and the inversion of the dispersion points has to be guided by the joint inversion of the HVSRs. As both of these methods are poorly sensitive to the compressional wave velocity |${V_P}$|⁠, we use a polynomial relation between |${V_P}$| and |${V_S}$|⁠. However, the commonly used relations proposed by Brocher (2005) have been obtained with rocks of various origins but not soils, as in the MAVM. Therefore, we used an ad-hoc polynomial relation that is valid for the MAVM. This latter has been obtained by gathering data from boreholes, cross-holes, penetration points, geotechnical and ultrasound measures (Zuñiga-Torres & Perton 2024):

where both velocities are expressed in km s⁻¹. When |${V_S} > 0.8\,{\rm{km\,s^{-1}}}$|⁠, the relation is similar to the one proposed by Brocher (2005), but below that limit, |${V_P}$| tends toward 1.5 km s⁻¹, as in water, and is always greater than that value because all the materials are fluid-saturated. It is important to underline that no direct measure of |${V_P} < 1.5\,{\rm{km\,s^{-1}}}$| have been reported. The mass density, which is not constrained by the inversion, is also deduced from Brocher's polynomial relations. We considered a medium purely isotropic for the DCs and HVSRs calculation. The inversion procedure has been described by Perton et al. (2020). We estimate the elastic properties below each station by minimizing the same objective function defined in this reference, that is, by comparing the observed and the theoretical HVSRs and the observed dispersion points with the theoretical dispersion curves for Love and Rayleigh surface waves. The only difference is that we considered only the group velocities until the 12th mode. The inversion uses a nonlinear scheme (Byrd et al. 1999; Matlab 2023) and is repeated iteratively by inserting new layers (Spica et al. 2016). Because of the considerable velocity contrast expected in the shallow subsurface observed in the Tulyehualco velocity profile (Fig. 2a), we do not constrain the velocity variation between the contiguous layers here. We started the inversion at the station closest to the Tulyehualco borehole, as it corresponds to the best-known velocity profile in the area. Generally, the procedure leads to 20 layers until the misfit becomes acceptable. Once we obtain a satisfactory fit, we follow the inversion of the elastic properties below the nearest station by considering the last velocity profile as the new starting profile.

5. RESULTS

We obtained a |$Vs$| velocity profile under 80 stations and discarded the other four because of their poor-quality data. Fig. 5 shows inversion results at one station of each array.

In the first column, we plot the HVSR targets in black, the convergence maximum and minimum curves for half the windows in red, and the best-fitting HVSRs in cyan. The plots are in log–log representation. We used the same frequency and amplitude ranges for the six arrays to help their comparison. We observe a fine agreement between the observed and final theoretical HVRSs.

We show the dispersion points for Rayleigh and Love waves and the best-fitting dispersion curves in the second and third columns. The velocity axes are on a logarithmic scale to better appreciate the fit at low velocities; the frequency axes are linear due to the harmonic behaviour of the frequency cuts. With this representation, we observe two separated areas, one at low velocity, that is, below 150 m s⁻¹, where the dispersion points are principally horizontally aligned, and another, above 300 m s⁻¹ with denser packets of scattered points. For several stations, we do retrieve dispersion points at low velocities and high frequency (>3 Hz). This area corresponds to short wavelengths, retrieval of which requires short interstation distances. In particular, we observe no dispersion point in this area for station 103, which belongs to the N06 network where the interstation distances are the largest. The two separated areas agree well with the theoretical dispersion curves. For the low-velocity area, the dispersion curves follow horizontal asymptotes. On the contrary, for the high-velocity area, the dispersion curves exhibit vertical asymptotes at regular cutoff frequencies followed by several rebounds. Between the two areas, the dispersion curves only transit once, explaining the dispersion points' scarcity. In addition to these qualitative agreements, which underline the importance of considering the multimode approach, we present the error maps made on HVSR and dispersion point fits in the supplementary material. Although we conducted a joint inversion, we present the error maps separately to interpret better the misfit errors, that is, average velocity misfit for the dispersion points and average amplitude misfit for the HVSR. These error maps are good indicators of the lateral resolution of the tomography. As discussed in Perton (2022), traditional checkerboard tests can only be conducted when the modes are identified and, when waves travel several wavelengths. These two hypotheses are not verified here. We also think the resulting maps should not be considered, as is usually the case, as velocity resolution maps but rather maps indicating where some rays might be missing to be confident of the regionalized velocities. More importantly, if a mode is wrongly identified, several modes are mixed inside the maps of surface wave velocity to be regionalized, and the tomography will be completely erroneous even if the checkerboard maps are good.

The inverted velocity profiles are shown in the fourth column with a log–log representation appropriate for the thin layers in the shallow part and thick layers in the deep part. They reach a depth of about 700 m. The first layers are a few metres thick and have higher velocities than those underneath. From the surface to the 50-m depth, we observe several huge velocity contrasts, inversions and generally a negative velocity gradient. The gradient is positive below that depth.

Finally, the station positions as the paths between station pairs contributing to the dispersion points are indicated with triangles and dashed lines in the last column.

We finally built a 3-D velocity model by interpolating all the 1-D velocity profiles with a natural neighbour Delaunay triangulation of class C₁ (Matlab 2023), which helps to increase the spatial coherence of the result. Fig. 6(a) shows this 3-D velocity model until 200 m depth, with a linear depth representation. Figs. 6(a), (b) and (c) show three vertical sections below the lines AA’, BB’ and CC’ shown in Fig. 6(a). The sections are shown up to 700 m in depth with a logarithmic depth representation, allowing the identification of numerous intercalations of slow and fast layers. We also show the 2017 static levels of five water wells provided by COANGUA (National Water Commission) with black dots and their projected positions.

6. DISCUSSION

We first discuss the HVSRs results and draw important considerations for sedimentary basin characterization and their effect of amplification. Then, we discuss the velocity model with the geological and hydrological models.

6.1. HVSRs

Fig. 4 shows a 4-D representation of the HVSRs in latitude, longitude, frequency and amplitude, allowing us to follow the dominant peak frequency's spatial coherence and validate the HVSRs' processing. However, the amplitudes do not vary as smoothly in space as the frequencies, and very localized anomalies have tremendous amplifications above ten. Fig. 7 shows a more classical spatial variability of the dominant peak amplification, which is of paramount importance in civil engineering. The whole area follows consequent amplifications above four. The map's southeast (SE) part has a stable amplification above 14, but the northwest (NW) part presents more heterogeneities, with amplifications between 6 and 20. These two areas also seem separated by a lower amplification line that follows the 0.5 Hz isocontour of the frequency of the dominant peaks. This line corresponds perfectly with the lake border shown in Fig. 1, drawn initially according to the water's presence at the surface.

To better understand the origin of the amplitude heterogeneity, despite the slow variation of the frequencies of the dominant peak, we present in Fig. 8 the theoretical HVSRs of several stations in the south area, all belonging to the network N05. All the curves having the same colour present similar HVSRs and velocity profiles. Consequently, we only present one representative inverted velocity profile for each coloured HVSR class. For all of the stations, the theoretical HVSRs agree well with the observed ones (small error in Figs S3 and S4), so the velocity profiles adequately explain the high amplification observed. The maximum distance between the stations is only 2.2 km, and the frequencies of the dominant peaks vary from 0.2 to 0.4 Hz, but the amplitudes vary from 8 to 18. Amplifications increase while the station positions become closer to the lake. Interestingly, the bases of the peaks (i.e. the amplitude of the curves below 0.2 Hz and above 0.4 Hz) are similar, and the highest amplifications seem to be built on top of these bases as an additional superimposed thin peak. This observation suggests that the extraordinary amplitudes observed in some stations result from the constructive addition of several HVSR peaks coinciding in frequency.

To validate this hypothesis, we analysed the velocity profiles, which all present several shear wave velocity inversions (Fig. 8). Generally, velocities are higher at the surface than at 20 m depth due to soil compaction. The velocities at 20 m depth are between 40 and 80 m s⁻¹, and in the underneath layer, above 100 m s⁻¹ (velocity ratio between two and four). After this minimum velocity, the general behaviour observed is that velocities increase with depth at a significative linear gradient of |${V_{s,z}} = 2.8\,{{\rm{s}}^{ - 1}}$|⁠. Several deeper velocity inversions are also observed, and high-velocity contrasts are then introduced because the velocities tend to retrieve the overall velocity gradient after the velocity inversion. These high contrasts are fundamental to trapping the waves and producing high amplification in the HVSRs. The velocity profiles associated with the highest HVSR amplifications (green curves) present more velocity inversions, particularly a strong one between 35 and 60 m depth. We then computed a theoretical HVSR considering the velocity profile of one of the highest HVSR (green velocity profile) but without this velocity inversion (cyan velocity profile) and obtained an HVSR similar to the one of the orange class, demonstrating consequently that this velocity inversion is responsible for an additive contribution to the large peak. Then, the highest HVSR values result from the superposition of several peaks produced by several velocity contrasts occurring at different depths.

The fact that the peaks all coincide with almost the same frequency is due to the nearly constant velocity gradient. Indeed, for the simplest case of a layer on top of a half-space, the classical formula is |$f = {V_s}/4h$| with f the peak frequency and |${V_s}$| and h the shear velocity and thickness of the top layer, respectively. In this expression, any multiplication of |${V_s}$| and h by the same factor leads to the same frequency. This first shows that we cannot independently retrieve the depth and velocity with the HVSR technique (non-uniqueness of the solution), and secondly, that the peak frequency does not depend on the properties of the underneath medium. The significant parameter is then the ratio |${V_s}/h$|⁠, which is constant for a medium with a constant velocity gradient. However, this formula is only valid for a layer on top of a half-space, and to correctly represent a medium with a velocity gradient, we can consider, in a first approximation, the mean velocity of the gradient profile in the top layer place of |${V_s}$|⁠: |$V_S^m = \frac{1}{h}\int_0^h {{V_{s,z}}{\rm{zdz}}} = {V_{s,z}}h/2$|⁠. Consequently, for a linear velocity gradient with |${V_{s,z}} = 2.8{{\rm{s}}^{ - 1}}$|⁠, the peak frequency is |$f = {V_{s,z}}/8 = 0.35\,{\rm{Hz}}$|⁠, approximately the frequency observed for all the HVSR curves in Fig. 8. This good agreement validates the hypothesis of the effect of a constant velocity gradient origin in the frequencies observed.

The physical mechanism leading to an addition of peak amplitude remains an open question requiring more theoretical development. Nonetheless, the HV corresponds to the division of the horizontal to the vertical energies (Perton et al. 2009; Sánchez-Sesma et al 2011), that is, positive values. It seems then natural that the peaks interfere constructively. Also, the highest HV amplitudes observed here are obtained at frequencies where the fundamental Rayleigh or Love group dispersion curves have almost vertical asymptotes and high energy (Fig. 2). Then, in a narrow frequency band, the modes carry short to large wavelengths that are sensitive to different depths. Contrary to what these peaks suggest, the HV peaks do not correspond to resonances, as they are not a direct measure of the energy of a displacement component. In fact, as shown in Fig. 2, the high values are due to low vertical energies. In this case, the HV peaks correspond to an energy projection, that is, most energy is taken from the vertical component and projected along the horizontal component, at surface. It is then a good hazard indicator because horizontal motion has more significant consequences on building than vertical motion.

If the frequency of the dominant peak varies spatially slowly in Fig. 3, it is first due to the slow variation of the overall gradient velocity. The depths of the several velocity contrasts are then of little importance. It is then of fundamental importance to consider the almost continuous variation of the velocity through depth in sedimentary basins for their correct characterization and then avoid using the HVSR technique to assess an interface depth between a homogeneous sediment layer and a hard basement.

In summary, the extraordinary HVSR amplitudes result from the constructive addition of several HVSR peaks coinciding in frequency because of a constant velocity gradient, typical of sedimentary basins. As other cities, for example Los Angeles, Vancouver, Tokyo, Seattle and parts of the San Francisco Bay area, sit on top of large sedimentary basins, these results are significant in understanding the origin of the high amplification.

6.2. Geological interpretation of the 3-D velocity model

In what follows, we interpret the 3-D |${V_S}$| model and comment on its agreement with other studies in the Chlaco nearby sub-Basin and hydrologic data.

As described in the previous subsection, the top layer presents high velocity, between 400 and 600 m s⁻¹. It is generally associated with basalts and andesites, as for the Tulyehualco-1 section (COMESA 2014), but it also results from the soil compaction and stiffening associated with tree roots and constructions. It is locally known as a ‘superficial crust’ (GCDMX 2004). Below, low velocities (<300 m s⁻¹) volumes appear as several horizontal layers. The first one is between 8 and 30 m depth and presents the lowest velocities (V_S < 100 m s⁻¹), stands east of the studied area, where lakes and canals remain. Then, we identify three deeper layers until 60 m. This limit corresponds to the piezometric level (black points in Fig. 6) of several wells reported by CONAGUA (National Water Commission). Because no water can be extracted above this depth, the permeability of the low-velocity anomaly is low. We associate the media between 8 and 60 m with the clay formation (as already reported in Campos-Enríquez et al. (1997), COMESA (2014) and González-Flores et al. (2020), which stands between 8 and 200 m of depth in the south of the MAVM. Clay-rich sediments are indeed material presenting extremely high void ratios (6–11) and water contents [> 400 per cent (Warren & Rudolph 1997; E. Ovando-Shelley et al. 2003)] and are generally considered as the aquitard (Lesser 2003). A deeper layer from 70 to 150 m, with velocities between 150 and 400 m s⁻¹, can also be identified below the hydraulic head and is associated with the aquifer. The aquifer may be semiconfined or free when water extraction wells are present (Lesser 2003; Pérez Centeno 2009). In the latter case, the hydraulic head is directly at the aquifer level and then these measurements restrain the interface aquitard/aquifer.

Beneath the depth of 150 m, the velocities are above 700 m s⁻¹ with a high-velocity anomaly of 1500 m s⁻¹ toward the N-W. The highest velocities are consistent with the presence of hard volcanic layers due to the Chichinautzin and Santa Catarina ranges and Cerro de la Estrella (Ortega‐Guerrero et al. 1993; Campos-Enríquez et al. 1997a; Juárez-Camarena et al. 2016a).

Bordering the North of our model, Montiel Rosado (1989) and González-Flores et al. (2020b) interpreted the result of a seismic reflection line. The same geological units (a consolidated top layer overlaying an aquitard of about 60 m thick, followed by an aquifer corresponding to volcanic permeable sand and a hard basement) have been identified, and the soft sediment layer thickness reduction toward the west is also observed.

The seismic tomography model only brings limited information to build a hydrological model and does not allow the separation of freshwater from saline waters. Resistivity, gravimetry and geological studies are of fundamental importance to better interpreting our results, as done in the nearby Chalco sub-Basin (Campos-Enríquez et al. 1997b; Vergara-Huerta & Aguirre 2020). Nonetheless, this 3-D shear wave velocity model gives an unprecedented vision of the complexity of the subsurface properties and may help locate the position of future wells by identifying the aquifer's deepest part.

7. CONCLUSIONS

We present an ambient seismic noise tomography of the Xochimilco area based on a joint HVSR and multimodal dispersion curves (DC) inversion. We first discuss the difficulties this method faces due to the complex geological setting. Subsequently, we underline the proper signal processing to obtain the high amplification peak of the HVSRs adequately. Similarly, we discuss the proper calculation of the theoretical counterparts, which requires integrating densely until large wavenumbers due to the unusual presence of extremely low shear wave velocities (< 50 m s⁻¹). We observe a spatially continuous decrease of the dominant peak frequency toward the lake interior but a heterogeneous amplification. By analysing the velocity profiles associated with the highest amplifications, we discovered that they result from the superposition of several resonance peaks coinciding in frequency. This particular superposition is caused by several low-velocity anomalies in a nearly constant linear gradient velocity. We also discovered that the overall gradient velocity in the subsurface varies slowly in the area. We also discussed that these low-velocity anomalies are due to a high water content both in the aquitard and aquifer. These are made of interbedded layers of soft and hard materials. This tomography also demonstrates the imperative necessity of conducting a multimode dispersion curves approach. According to the theoretical kernels, the energy is spread over the fundamental and higher modes of the surface waves, particularly between 50 and 300 m s⁻¹ in the whole frequency range. Below 0.7 Hz, the fundamental mode is preponderant, but at higher frequencies, the first higher mode is more energetic. All these observations are of fundamental importance for the correct seismic mitigation.

ACKNOWLEDGEMENTS

Project IISGCONV-181–2016 ‘Comprehensive Seismic Monitoring System and Generation of Real-Time Intensity and Damage Maps for Mexico City’ (in Spanish). SECTEI, CDMX. Technical manager Dr Leonardo Ramírez Guzmán.

DATA AVAILABILITY

The data underlying this article were provided by the Engineer Institute of the National Autonomous University of Mexico by permission. Data will be shared on request to the corresponding author with the permission of Dr Leonardo Ramírez Guzmán, the manager in charge.

REFERENCES