Articles | Volume 14, issue 5
Research article
23 May 2023
Research article |  | 23 May 2023

Probing environmental and tectonic changes underneath Mexico City with the urban seismic field

Laura A. Ermert, Enrique Cabral-Cano, Estelle Chaussard, Darío Solano-Rojas, Luis Quintanar, Diana Morales Padilla, Enrique A. Fernández-Torres, and Marine A. Denolle

The sediments underneath Mexico City have unique mechanical properties that give rise to strong site effects. We investigated temporal changes in the seismic velocity at strong-motion and broadband seismic stations throughout Mexico City, including sites with different geologic characteristics ranging from city center locations situated on lacustrine clay to hillside locations on volcanic bedrock. We used autocorrelations of urban seismic noise, enhanced by waveform clustering, to extract subtle seismic velocity changes by coda wave interferometry. We observed and modeled seasonal, co- and post-seismic changes, as well as a long-term linear trend in seismic velocity. Seasonal variations can be explained by self-consistent models of thermoelastic and poroelastic changes in the subsurface shear wave velocity. Overall, sites on lacustrine clay-rich sediments appear to be more sensitive to seasonal surface temperature changes, whereas sites on alluvial and volcaniclastic sediments and on bedrock are sensitive to precipitation. The 2017 Mw 7.1 Puebla and 2020 Mw 7.4 Oaxaca earthquakes both caused a clear drop in seismic velocity, followed by a time-logarithmic recovery that may still be ongoing for the 2017 event at several sites or that may remain incomplete. The slope of the linear trend in seismic velocity is correlated with the downward vertical displacement of the ground measured by interferometric synthetic aperture radar, suggesting a causative relationship and supporting earlier studies on changes in the resonance frequency of sites in the Mexico City basin due to groundwater extraction. Our findings show how sensitively shallow seismic velocity and, in consequence, site effects react to environmental, tectonic and anthropogenic processes. They also demonstrate that urban strong-motion stations provide useful data for coda wave monitoring given sufficiently high-amplitude urban seismic noise.

1 Introduction

Near-surface geological structures and soil properties are important determinants of seismic hazard (e.g., Field2000). Shallow, poorly consolidated sediments can strongly amplify long-period seismic waves, and the strong impedance contrast to the underlying bedrock can trap energy and vastly prolong ground motion (e.g., Roten et al.2008; Cruz-Atienza et al.2016). Therefore, considerable effort is invested to determine shallow sediment properties in urban areas (see Foti et al.2019,  for a review). One of the decisive quantities controlling site response is the shallow shear wave velocity, which is also used as a hazard assessment parameter in the form of vs30, or shear wave velocity averaged in the top 30 m depth. Shallow seismic velocities react strongly to environmental variations, as has been documented by time-lapse imaging (Bergamo et al.2016) and ambient noise monitoring (e.g., Steinmann et al.2021). Ambient noise monitoring is based on comparing short-term interferometric measurements like cross-correlation or deconvolution of continuous seismic data to a long-term reference (Obermann and Hillers2019). Under the assumption that the coda of the resulting waveforms is predominantly sensitive to changes in the elastic medium, one can measure subtle relative advances and delays in the current waveform compared to the reference; these are approximately linearly related to relative velocity changes dvv. Using this technique reveals that seismic velocity varies with groundwater level (Sens-Schönfelder and Wegler2006; Lecocq et al.2017; Fokker et al.2021; Rodríguez Tribaldos and Ajo-Franklin2021), precipitation, soil moisture and snow load (Obermann et al.2014; Voisin et al.2016; Wang et al.2017; Donaldson et al.2019; Andajani et al.2020; Mordret et al.2020; Feng et al.2021; Illien et al.2021), ground temperature (Richter et al.2014), thawing of the permafrost (Ajo-Franklin et al.2017; James et al.2019; Lindner et al.2021), and even centimeter-scale layers of soil freezing (Steinmann et al.2021). Droughts can induce longer-term changes of seismic velocity (Clements and Denolle2018; Mao et al.2022), as can soil compaction (Taira et al.2018). Although the reported changes are usually small, on the order of approximately 1 % peak-to-peak amplitude or 0.01 % yr−1–0.1 % yr−1 trend, they clearly show that shallow sediment properties are time dependent.

A second phenomenon relevant to site response is the non-linear behavior of soft near-surface sediments subject to large dynamic strains, including shear modulus reduction and plastic deformation (Wu and Peng2012; Oral et al.2019; Bonilla et al.2019). Numerous recent studies utilizing ambient noise interferometry have reported temporary velocity drops consistent with shear modulus reduction during the shaking from moderate to large earthquakes, generally followed by post-seismic relaxation that is approximately linear with the logarithm of time (e.g., Brenguier et al.2008b; Hadziioannou et al.2011; Wu and Peng2012; Froment et al.2013; Hobiger et al.2016; Wang et al.2017). Most studies found small velocity drops on the order of 0.1 %–1 % following ground shaking. However, much larger reductions on the order of 1 %–10 % and more were reported in studies which explicitly targeted shallow structure on the order of 100 m below the surface, as would be expected for non-linear elasticity (Nakata and Snieder2011; Viens et al.2018). The velocity reduction is reported to be even stronger during the shaking itself (Bonilla et al.2019; Bonilla and Ben-Zion2020); however, most studies lack sufficient time resolution to capture this short-lived effect.

In the present study, we investigate both linear and non-linear changes of the seismic velocity underneath Mexico City. Mexico City has suffered devastating ground shaking, particularly due to the response of lacustrine clay deposits in the center of the Mexico City basin (Anderson et al.1986; Singh et al.1988a; Sahakian et al.2018; Arciniega-Ceballos et al.2018; Cruz-Atienza et al.2016), and continues to face a high seismic hazard. Mexico City is also affected by extremely rapid ground subsidence (Cabral-Cano et al.2008; Chaussard et al.2021). The main motivation for our study is to understand how environmental factors and strong ground motions from earthquakes influence the seismic velocity structure of the shallow to intermediate sediments. Furthermore, we aim to demonstrate that it is possible to monitor such temporal changes using data from continuous recordings of urban seismic noise at a relatively sparse strong-motion sensor network with seismic interferometry. We characterize the observed velocity changes through physics-based modeling and probabilistic inversion of key parameters like velocity drop and sensitivity to surface temperature.

In the following, we describe the data (Sect. 2) and the processing approach we used to overcome the particular challenges of urban seismic noise (Groos and Ritter2009; Schippkus et al.2020, Sect. 3). In Sect. 4, we introduce our model for velocity changes and the probabilistic inversion of model parameters. Finally, we present and discuss the observations and inversions (Sect. 5) and conclude with an outlook and recommendations for further research (Sect. 6).

Figure 1Shaded relief (Farr et al.2007; NASA JPL2013), geotechnical zonation (Gobierno de la Ciudad de México2017) and location of seismic stations in the Valley of Mexico region (Quintanar et al.2018; MASE2007; Roult et al.2010). The lake zone (brown outline) and transition zone (orange outline) include locations where the shallow subsurface is characterized by quaternary lacustrine sediment. Strong-motion and broadband sensors at Geoscope station G.UNM (green triangle) were used for comparing velocity changes assessed with different instruments and for analysis of the velocity changes. Strong-motion sensors of the Red Sísmica del Valle de México (red triangles) and broadband stations of the temporary Tectonic Observatory (TO) deployment (blue triangles) and G.UNM were used for analysis of velocity changes. Strong-motion stations of the Red Acelerográfica del Instituto de Ingeniería (RAII-UNAM, black rhombs) were used for determining additional peak ground acceleration values during the 2017 Mw 7.1 Puebla earthquake.

Figure 2(a) Acceleration spectrogram of urban noise at G.UNM (North component), broadband seismometer. Prior to June 2008, the sensor was operated at lower gain. The spectrogram illustrates high noise levels characteristic for the urban location. Faint changes in the spectrum coincide with the September 2017 earthquake (left cyan line and arrow above the panel), as well as a marked drop in noise level during the Covid-19 pandemic (right cyan line and arrow: first announcement of anti-pandemic measures). (b) Clustering of short-term N–N autocorrelation waveforms in the frequency band 2–4 Hz. Clustering results suggest that autocorrelation waveforms change on a day–night rhythm, a weekly rhythm (not visible here), and an annual switch to and from daylight savings time. Similar urban patterns appear in clustering results at all stations, mostly at frequencies of ≥0.5Hz. The black vertical lines indicate the 2017 earthquake and announcement of anti-pandemic measures. The color scale was chosen for accessibility (Crameri et al.2020).


2 Study area and data

The National Seismological Service of Mexico operates a state-of-the-art seismic network in the Valley of Mexico, the Red Sismica del Valle de México (RSVM; Quintanar et al.2018). Here, we use data from RSVM strong-motion sensors that were mostly installed in 2017. Strong-motion sensors are chosen for dense deployment in some urban areas with high seismic hazard and thus provide a valuable data source for ambient noise monitoring despite their comparably low sensitivity (e.g., Tokyo metropolitan area, Seattle, Berkeley Borehole network in the San Francisco Bay area; Kasahara et al.2009; Viens et al.2018; Bonilla et al.2019). We first compare results obtained from a co-located seismometer and accelerometer to verify that the urban noise at periods of 2 s and shorter is well captured by both types of sensors and that results are consistent. We then focus on continuous strong-motion recordings at 12 locations representative of the different site conditions in Mexico City, including station locations on soft, intermediate and hard sites, as defined by Quintanar et al. (2018). We supplement these with broadband sensor observations from the Geoscope network site G.UNM and four stations of the temporary Tectonic Observatory deployment (MASE2007).

The station G.UNM includes a co-located STS-1 seismometer (broadband sensor) and Kinemetrics EpiSensor accelerometer (strong-motion sensor) and has been recording continuously since 1995 for the broadband sensor and since 2013 for the strong-motion sensor on the main campus of the National Autonomous University of Mexico (green triangle in Fig. 1). The upper panel of Fig. 2 shows a spectrogram obtained from the autocorrelations of the North component of the broadband seismometer between 1995 and 2021. It illustrates that urban-noise levels around 1 Hz surpass Peterson's New High Noise Model of −120 dB (Peterson1993). Furthermore, instrumental effects are clearly visible, such as the change in instrument gain during 2008. The recorded urban noise is not stationary; a group of short-lived spectral peaks appears only in 2015. In 2020, the noise level drops at the onset of the Covid-19 pandemic (see Pérez-Campos et al.2021). We will later on compare observed and modeled results from the two co-located sensors to assess the use of strong-motion stations for coda wave monitoring (see Sect. 5.1).

For analyzing velocity changes, we focus on urban seismic noise autocorrelations. Autocorrelations have been successfully used for monitoring earthquake damage and climate effects on the near-surface velocities in previous studies (e.g., De Plaen et al.2019; Sánchez-Pastor et al.2019; Feng et al.2021). Moreover, coherency of the ambient noise between stations is poor at frequencies above 1 Hz, which may be due to the near-station sources of high urban noise, the strong attenuation in the basin sediment (e.g., Cruz-Atienza et al.2016) and the comparably low sensor sensitivity.

3 Measuring velocity changes in an urban setting

We remove large global earthquakes of Mw=7 and above according to the global CMT catalogue using magnitude-dependent window lengths following the approach of Ekström (2001). Local earthquakes in a 5 radius from the ISC catalog (International Seismological Centre2022; Bondár and Storchak2011) are removed by cutting out a window starting 10 s before the direct P wave using iasp91 (Kennet1991) and ending 60 s after a 1 km s−1 surface wave. We cut data into 8 h segments; detrend, taper and pre-bandpass-filter them; and remove the instrument responses to correct the data to ground acceleration (including the seismometer for direct comparison of the results). Finally, we correlate 20 min windows of the data overlapping by 10 min. We save all single-window correlations for further processing in hdf5 format (Folk et al.2011). The pre-processing and correlation are performed with the Python module ants, which uses obspy (Beyreuther et al.2010; Krischer et al.2015) for instrument correction, filtering and other seismic data-processing tasks (see the code and data availability section for details on code availability).

3.1 Clustering correlation waveforms and selective stacking

We adopt a novel processing strategy for ambient noise correlations proposed by Viens and Iwata (2020). It is based on the premise that variable incident ambient noise conditions result in different correlation waveforms. By grouping single-window autocorrelations into clusters and stacking them selectively, we aim to increase their temporal coherence and to reduce the effect of temporally varying noise sources. Clustering is performed by applying Gaussian mixture models after reducing the data dimensionality through principal component analysis (PCA), and the Bayesian information criterion (BIC) is used to determine the optimal number of clusters, balancing misfit and model complexity (Viens and Iwata2020).

We modify the original approach in several ways in order to adapt it to long-term urban noise. Details and rationales are provided in the Supplement. The modifications can be summarized as follows: (i) we apply clustering per octave frequency band to account for narrow-band cultural sources; (ii) we pre-determine the principal component axes on a subset of the data due to the large data volume; (iii) we normalize by dividing each waveform by its maximum absolute value; (iv) we fix the number of clusters a priori at an optimum determined in a preliminary clustering run, here four clusters; and (v) we rely on cultural patterns with respect to local time to label the clusters as day, night, noise and other. Here, noise refers to the smallest cluster of sporadically appearing, transient disturbances. The other cluster is required to fulfill the optimum number of four clusters. It mostly coincides with the transition between day- and nighttime. The result for the N–N component of the seismometer at G.UNM is shown in Fig. 2. The urban day–night rhythm emerges for all stations at frequencies of 0.5 Hz and above, as well as for several stations inside the basin for frequencies between 0.25 and 0.5 Hz. We tested the effect of using the amplitude-unbiased phase cross-correlation (Schimmel et al.2011; Sánchez-Pastor et al.2018) on clustering and found that, although the clustering results change, the day–night rhythm is still clearly visible (Fig. S1 in the Supplement). This suggests that processing schemes aimed at equalizing noise amplitudes fail to suppress the effect of temporally varying noise sources in an urban environment. We therefore use the clustering approach for the subsequent analysis. We stack the short-term correlations for the daytime clusters, which are the largest in number and the most consistent with time, over a duration of 10 d. We handle the short-term correlation data, filtering, clustering and stacking with the Python module ruido (see the code and data availability section).

3.2 Stretching measurement

We use the stretching method to measure relative changes in velocities (e.g., Sens-Schönfelder and Wegler2006). A preliminary comparison to the moving window cross-spectral method (e.g., Takano et al.2020) showed good overall agreement. We use a multiple-reference approach due to the lack of long-term waveform coherence in our observations; details are provided in the Supplement. Multiple-reference approaches have been previously used (in somewhat different forms), e.g., by Sens-Schönfelder et al. (2014) and Donaldson et al. (2019). The stretching is performed in four frequency bands (0.5–1, 1–2, 2–4 and 4–8 Hz) and in two windows on each stack, one of which extends from 4 to 10 times the longest period of interest (defined by the frequency band) and the other from 8 to 20 times the longest period of interest (e.g., for a 0.5–1 Hz observation, we measure velocity changes in the 8–16 and 16–40 s windows). Figure 3 shows results at station G.UNM (0.5–1 Hz) using the coda window at 16–40 s. We note that the information contained in the single-station cross-correlations (mixed components) is visually consistent with the pure autocorrelations (Fig. 3). We thus proceed with the pure autocorrelations, mostly because their interpretation is conceptually simpler; this also reduces computational requirements for the following inversion. We also note that data quality as measured by the correlation coefficient between reference trace and current trace after stretching (CC_best; shown by color hue) increases markedly after the sensor at G.UNM was set to high-gain mode in 2008. For further analysis, we discard data points with CC_best <0.6.

Figure 3Velocity changes at station G.UNM, all unique channel pairs, 16–40 s lag. Channels labeled BHE, BHN and BHZ are the East, North and Vertical channels recorded by the seismometer. Color hue shows the correlation coefficient of the current and reference trace after stretching, with light colors indicating low and dark colors indicating high correlation.


4 Modeling velocity changes

Visual inspection of the time series in Fig. 3 suggests that they contain different patterns, specifically seasonal changes, a co-seismic drop coinciding with the 2017 Mw 7.1 Puebla and 2020 Mw 7.4 Oaxaca earthquakes, and post-seismic recovery. We model these processes by considering the influence of surface temperature, precipitation, earthquake shaking and healing on the seismic velocity. In addition, as we note that the seismic velocity generally increases over time, we include a linear trend. For simplicity, we sum the submodels as a linear combination (see also Hobiger et al.2014, 2016; Wang et al.2017; Donaldson et al.2019; Feng et al.2021), although there may be interactions that lead to non-linearity, e.g., between water content and temperature (Sens-Schönfelder and Eulenfeld2019). We use a Markov chain Monte Carlo (MCMC) inversion to infer the unknown model parameters. We assume that the changes observed at the surface are dominated by changes in Rayleigh wave phase velocity, which are mainly sensitive to changes in shear wave velocity at depth and are mostly insensitive to changes in P-wave velocity (see Sect. 5.2).

4.1 Earthquake effects

The 2017 Mw 7.1 Puebla earthquake caused strong ground motion in Mexico City, reaching peak ground accelerations over 15 % g (Alberto et al.2018) and a sharp acceleration of subsidence in various locations (Solano-Rojas et al.2020). We observe velocity drops coinciding with this earthquake, followed by approximately logarithmic recovery at most stations. More subtly, a similar signal is observed at several stations following the 2020 Mw 7.4 Oaxaca (Sta. María Xadani) earthquake – see Fig. 5. Physical processes causing the co-seismic velocity drop may be caused by plastic deformation (fractures opening due to shaking-induced dynamic stresses that exceed the geologic material yield strength; Gassenmeier et al.2016; Bonilla et al.2019) or by a non-linear mesoscopic elasticity that describes the loss and reestablishment of chemical bonds or capillary bridges that change frictional contacts (Sens-Schönfelder et al.2018). Snieder et al. (2017) proposed the following model for the long-term post-seismic recovery described as the superposition of exponential relaxation mechanisms with different relaxation timescales τ, ranging from τmin to τmax, after the time of the earthquake tquake:

(1) c ( t ) - c 0 c 0 = s c 0 ln τ max τ min - 1 τ min τ max 1 τ exp ( - t / τ ) d τ for  t t quake .

Here, s is a negative value indicating the drop of seismic velocity from its previous value c0. For intermediate timescales, (τmintτmax), this model captures the approximately logarithmic time dependence, while it tends towards c0 for large times and is finite for t=0 where a purely logarithmic recovery is not defined. Hypothesizing that the minimum timescale τmin is below the temporal resolution of our measurements, we set τmin=0.1s and fit both τmax and the co-seismic step drop s through the inversion described below. Velocity changes due to non-linear elasticity are depth dependent (e.g., Wang et al.2019). For the earthquake-induced non-linearity, we account for depth dependence only insofar as we invert data from different frequency bands separately (Obermann et al.2014; Sawazaki et al.2016; Wu et al.2016).

4.2 Seasonal effects

The following sections describe how we model the effects of surface temperature (Sect. 4.2.2) and precipitation (Sect. 4.2.3) on Rayleigh wave phase velocity. In both cases, we use analytical solutions to diffusion-like problems for a homogeneous half-space to model their effect on vs from the surface to 2 km depth in terms of (a) thermoelastic stress and (b) pore pressure. In particular, we use (a) the solution of Berger (1975) as formulated by Richter et al. (2014) and (b) the solution of Roeloffs (1988) as formulated by Talwani et al. (2007). In the thermoelastic model, annual and sub-annual periodic surface temperature changes diffuse through the shallow subsurface. In the pore pressure model, rain leads to a sudden increase in pore pressure near the surface, which then diffuses towards depth. We compute surface wave sensitivity kernels with surf, a Python package based on the Takeuchi and Saito (1972) solutions to the surface wave eigenproblem in layered, anisotropic elastic media (Fichtner2020). This is done using station-specific 1-D velocity profiles and assuming an isotropic medium (see Sect. 4.2.1). Finally, we integrate the depth-dependent vs change to obtain the predicted surface wave phase velocity (c) changes.

The solutions at depth in the homogeneous half-space depend on (a) thermal diffusivity and (b) hydraulic diffusivity of the sediments, which are not well known and in the case of (b) can vary by several orders of magnitude (Roeloffs1996). Inverting for these parameters probabilistically would be challenging because it requires re-evaluating the diffusion terms at all depths for each iteration. Instead, we run the inversion repeatedly for six homogeneous half-spaces characterized by (a) three and (b) two trial diffusivities in the ranges of (a) 0.15 to 2mm2 s−1 (Richter et al.2014) and (b) 0.0001 to 4m2 s−1 (Roeloffs1996). We retain the diffusivity for each site that produces the best fit across all components and frequency bands.

4.2.1 Site-specific velocity profiles

We need estimates of shear wave velocity vs, compressional wave velocity vp and density ρ to establish the effect of velocity changes at depth on surface wave velocity changes dcc. To account for the different subsurface geologies at each station, we use a classification of hard, intermediate and soft, as provided by Quintanar et al. (2018), and we consider the aquitard thickness information from Solano-Rojas et al. (2015). We then construct the 1-D profile of vs, vp, and ρ at each site as follows:

  1. assign a depth to bedrock value, namely hard sites 0 m, intermediate sites 100 m and soft sites 300 m, coarsely representing the lower boundary of volcanic and alluvial sediments, which generally deepens towards the center of the basin;

  2. assign a site-dependent lacustrine sediment thickness (clay) to the intermediate and soft sites derived from the thickness of the upper aquitard (Solano-Rojas et al.2015);

  3. assign velocities to clays, volcanic and alluvial sediments, and bedrock according to Table 1, where sediment 1 is the upper half of the sediment column (alluvial), sediment 2 is the lower half (volcanic sediments), bedrock 1 is considered up to 1000 m below the sediment, and bedrock 2 is considered at greater depths.

The values for geologic structures and seismic properties are based on a synopsis of the works of Pérez Cruz (1988); Chavez-Garcia and Bard (1994); Singh et al. (1995, 1997) and Shapiro et al. (2001), which serve as references for modeling the seismic velocity structure of the Mexico City basin (Cruz-Atienza et al.2016; Asimaki et al.2020). The velocity values suggest that Poisson's ratio approaches 0.5 for the lacustrine sediments. Figure S2 shows that our rule-based velocity model is in reasonable agreement with the shear wave velocity profiles based on well logs from Singh et al. (1995, 1997). Future investigations would benefit from a more detailed 3-D velocity model of the basin and surroundings.

Table 1Compilation of approximate elastic properties based on the synopsis of Pérez Cruz (1988), Singh et al. (1995, 1997) and Shapiro et al. (2001).

Download Print Version | Download XLSX

As the surface wave sensitivity kernels in Fig. S3 illustrate, the model leads to the following behaviors. At hard sites, the sensitivity is spread over the shallowest 1000 m, and the peak sensitivity shifts upwards with frequency from about 450 m in the lowest to about 50 m in the highest frequency band. At soft sites, the sensitivity is concentrated in the shallowest 50 m for the lowest and the shallowest 5 m for the highest frequency band but is, in any case, inside the low-velocity sediment. At intermediate sites, the behavior is similar to hard sites for the lower frequencies up to 2 Hz and similar to soft sites for the higher frequencies. Our single-station analysis presumes a laterally homogeneous structure near each station. While this is a rather crude assumption, we believe that it captures important aspects of the basin, particularly the very shallow sensitivity of our observations at soft sites.

4.2.2 Surface temperature effects

We adopt the approach of Richter et al. (2014) to compute the thermoelastic stress, neglecting variations at greater depth. The shear wave velocity change depends linearly on the temperature T at location x, depth z and time t,

(2) δ v s v s ( x , z , t ) = s t T ( x , z , t ) ,

(see Eq. 14 of Richter et al.2014), where we summarize material parameters into a depth-independent temperature sensitivity as follows:

(3) s t = 2 b α ρ v 2 σ c ,

(see Eq. 12 of Richter et al.2014), where b=1+ν1-ν for S waves; ρv2σc describes the change of shear modulus with respect to stress, i.e., the non-linear elastic rheology; ν is Poisson's ratio; and α is the linear thermal expansion coefficient. The single factors of this sensitivity are not particularly well known, so we fit the product st during the inversion described below. We use surface temperature measured at the meteorologic network of the Universidad Nacional Autónoma de México (UNAM; Instituto de Ciencias de la Atmósfera y Cambio Climatico2022). We expand the temperature curves into Fourier series with five terms using the nearest available meteorologic station to each seismic station. This proved to be important for the model fit because temperature variations with sub-annual periods also affect dcc. Finally, we convert the shear wave velocity change at depth to dcc using the phase velocity sensitivity kernels for the fundamental-mode Rayleigh waves. During the inversions, we use three trial values of thermal diffusivity κt in the range of 1.5×10-7 to 2×10-6m2 s−1 and select the κt at each station that produces the minimum average misfit over all components and frequency bands.

4.2.3 Hydrological effects

Quantitative interpretation of observed velocity changes in terms of hydrology is a matter of current research and can be highly site specific; it requires measuring or estimating the hydrologic changes (usually pore pressure) and estimating the sensitivity of the observed velocity change to these. We compile an exemplary, non-exhaustive selection of recent approaches in Table 2. Here, we estimate the pore pressure changes from precipitation data, use a depth-varying medium response to stress based on granular-media theory and translate changes in vs to dcc with surface wave sensitivity kernels.

We compute the pore pressure change using the impulse response to the hydraulic head change derived by Roeloffs (1988) for a homogeneous half-space, assuming test values for hydraulic diffusivity κh of 1 m2 s−1, representing relatively impermeable unconsolidated sandy sediments, and 10−4m2 s−1, representing permeable unconsolidated clay sediments (Roeloffs1996). Note that the sensitivity of the modeled dcc to hydraulic diffusivity is low because of the depth integration for surface wave phase velocity change and the high values of Poisson's ratio at the sites (Fig. S4). We estimate the change in vs due to pore pressure changes using granular-media theory (Mavko et al.2020). Saul and Lumley (2013) formulated the effective moduli for effective pressure peff, based on which

(4) 1 v s v s p eff = 1 6 p eff ,

(see also Takano et al.2017), where peff=Pov-nPp with overburden pressure Pov and pore pressure Pp. We estimate overburden pressure as lithostatic pressure from the 1-D density models described in Sect. 4.2.1 and unperturbed pore pressure as hydrostatic pressure, using porosities of 0.6 and 0.2 for clay and any other sediment, respectively, following Ortega and Farvolden (1989) and assuming n=1. Granular-media theory was previously used by Rodríguez Tribaldos and Ajo-Franklin (2021) to model velocity changes in an unconfined aquifer. Not all sites in our study have the same hydrogeological properties, but Takano et al. (2017) found that granular-medium theory can approximately explain observations of stress sensitivity over a large range of depths and geologic materials. Equation (4) tends towards infinity for effective pressures approaching zero. We mitigate this by introducing a minimal effective pressure p0 as a parameter in the inversion, i.e., peff=maxp0,Pov-nPp. The resulting dcc is strongly sensitive to p0.


Table 2Comparison of selected approaches to model hydrological effects on seismic velocity (not an exhaustive list). SW6, Sens-Schönfelder and Wegler (2006); L17, Lecocq et al. (2017); W17, Wang et al. (2017); RT21, Rodríguez Tribaldos and Ajo-Franklin (2021); F21, Fokker et al. (2021)

n/a: not applicable

Download Print Version | Download XLSX

4.3 Probabilistic inversion

Considering the superposition of the temperature, precipitation, earthquake and linear trend, we aim to model the phase velocity change dcc as

(5) d c c = f temp ( κ t , s t ) + f rain ( κ h , p 0 ) + f seismic ( Δ c , τ max ) + f lin ( a , b ) ,

where κt,κh are the thermal and hydraulic conductivity, respectively; st is the sensitivity of the shear wave velocity to thermoelastic stress; p0 is the minimal overburden pressure; Δc is the co-seismic drop in observed velocity (which we assume is surface wave phase velocity); τmax is the maximum timescale of exponential recovery in the slow-dynamics model of Snieder et al. (2017); and a and b are the slope and offset of the linear trend. Values for the remaining variables are taken from the literature. In particular, we use fluid volume fractions of 0.6 and 0.2 for lacustrine sediments and all other materials, respectively, from Ortega and Farvolden (1989) as an approximation for porosity. Furthermore, we use approximate bulk moduli of 2.5 and 35 GPa for the pore water and rock matrix, respectively. We use these values to estimate Skempton's B following Roeloffs (1988), and we adopt their estimate of the undrained Poisson ratio. We conduct an MCMC inversion to determine the parameters st, p0, Δc, τmax, a and b using the Euclidean distance between the observed and modeled dcc with the emcee Python package (Foreman-Mackey et al.2013). Details of the initialization and convergence of the emcee runs are given in Sect. S4 in the Supplement. Median models are shown in Fig. 5 for the lag window of 4–8 s of the 2–4 Hz frequency band; results in terms of median models and model ranges for all frequency bands and lag windows are shown in Figs. S5 and S6 in the Supplement.

5 Results and discussion

5.1 Comparison of co-located sensor results

Comparison between observed dcc at the co-located seismometer and accelerometer of the station G.UNM shows overall strong consistency (Fig. 4). Correlation coefficients between observed time series range from 0.83 to 0.99, with the exception of the East component at 0.5 and 2.0 Hz, where a gap in highly coherent observations results in an offset of dcc between the seismometer and accelerometer after the Puebla earthquake. Consequently, we discourage the analysis of single components; a synoptic analysis or weighted average of components, as suggested by, e.g., by Hobiger et al. (2014), is preferable. We conclude that strong-motion stations yield overall good results for dcc studies in urban settings with high-amplitude anthropogenic noise comparable to G.UNM.

Figure 4Comparison of dcc results from the co-located broadband seismometer and accelerometer at station G.UNM. (a) Data for the duration of operation of both sensors (since 2013) for all components (E, N, Z) of autocorrelations and all frequency bands (0.5–1, 1–2, 2–4 and 4–8 Hz). (b) Median posterior models obtained by MCMC inversion, constrained by the data on the left. Both the data and the model are in good agreement between the sensors, with the exception of the East component.


5.2 Observed and modeled velocity changes for RSVM accelerometers and G.UNM seismometer

Results at 2–4 Hz for the strong-motion stations of RSVM and the seismometer at G.UNM (Fig. 5) illustrate the seasonal dcc variations and the drops and recoveries for the 2017 Puebla (19 September 2017) and 2020 Oaxaca (23 June 2020) earthquakes. Despite its simplicity, our model captures the behavior of the observed velocity changes reasonably well. For RSVM and G.UNM at 2–4 Hz, with a lag window of 4–10 s, the mean correlation coefficient (CC) between modeled and observed dcc is 0.77, and the median CC is 0.84. (0.5–1.0 Hz: 0.68 and 0.75; 1–2 Hz: 0.67 and 0.73; 4–8 Hz: 0.69 and 0.76). The inversion failed to converge in approximately 10 % of cases, which we include nevertheless (see Figs. S5 and S6 in the Supplement). We exclude results above a misfit threshold set at the bottom quartile from further analysis. These models are shown by dashed white lines in Fig. 5.

Our model generally fits better for stations inside the basin, including the transition zone and lake zone stations, than for those outside. This is reflected by a consistently higher CC between observed and modeled dcc for stations inside the basin and by a slightly lower average misfit.

We propose three reasons for the better performance of the model inside the basin. First, incident ambient noise may be more stable near the city center due to repetitive traffic patterns and a generally higher amplitude. A detailed analysis of noise source conditions is out of scope of this work but may benefit future urban monitoring studies. Second, the assumption that observations are dominated by scattered surface waves may be more appropriate for basin stations located on stratified sediment with strong impedance contrasts. Third, better a priori information is available about the basin subsurface structure compared to areas outside of the basin thanks to the detailed thickness information of the low-vs aquitard (Solano-Rojas et al.2015) and to shallow profiles of vs from well logs (Singh et al.1995, 1997). This results in what are likely to be more accurate station-specific surface wave sensitivity kernels inside the basin, which illustrates the value of a priori site-specific geologic and hydrologic information for quantitative shallow-structure monitoring with ambient noise.

Several previous studies have applied detailed physics-based models to dcc, usually focusing on a small number of selected stations (Tsai2011; Richter et al.2014; Lecocq et al.2017; Fokker et al.2021; Illien et al.2021, 2022) or time series (Rodríguez Tribaldos and Ajo-Franklin2021). Previous array-wide analyses have mostly focused on empirical transfer functions between hydrologic and meteorologic parameters and observed dcc (e.g., Wang et al.2017; Donaldson et al.2019). Here, we take a partially physics-based approach to the scale of a sedimentary basin and metropolitan seismic array, using a comprehensive set of time series in terms of frequency (0.5–8 Hz in octaves) and spatial components (E, N, Z) and integrating multiple processes influencing dcc (poroelastic stress, thermal stress, non-linear elasticity with co-seismic drop and recovery). It is straightforward to extend our approach to other surface wave modes (Love waves, higher modes). Here, we chose to limit our analysis to fundamental-mode Rayleigh waves, since limited information on subsurface structure is available, and we have limited knowledge of the surface wave modal representation in scattered waves.

A current limitation of our model is that we assume a linear superposition of different processes affecting dcc (see Illien et al.2022, for a discussion on this point). Nevertheless, it is the first step to a comprehensive model for dcc based on simplified physics at the basin scale.

Figure 5Time series of relative velocity changes (circles) and median models (white lines) in the frequency band of 2–4 Hz for autocorrelations of the East, North and Vertical components of RSVM stations and G.UNM. Stations are arranged by topographic elevation outside the basin and by aquitard thickness from Solano-Rojas et al. (2015) for stations inside the basin. Stations at a lower elevation (near the basin edges) tend to be located on alluvial sediment, while those at a high elevation are on bedrock. Color coding shows the aquitard depth in shades of yellow to black (minimum = 0 m, maximum = 75 m) and elevation in shades of yellow to green. Vertical dashed red lines indicate the timing of the 2017 and 2020 earthquakes, while the horizontal dashed gray line indicates the boundary between stations at hard sites and stations at intermediate and soft sites. Vertical grid spacing is 1 % velocity change.


5.3 Parameters controlling seasonal velocity variations

The parameters relevant to seasonal variations in our model are (i) thermal diffusivity, (ii) hydraulic diffusivity, (iii) the sensitivity of vs to thermal stress st (which summarizes non-linear and linear elastic properties), and (iv) the minimum effective pressure p0. For practical reasons, only a few values were tested for (i) and (ii). Therefore, we only present the inference of st and p0 by the inversion. We observe variability between the Z, N and E components for both parameters, as well as between the frequency bands. Depending on the station and frequency band, either temperature or precipitation emerges as the dominant seasonal control from our inversion; this is apparent from results shown in the Supplement.

In previous studies, single-station correlations have been interpreted as multiply scattered body waves (e.g., Richter et al.2014; Sánchez-Pastor et al.2018), whereas we interpret them as surface waves. This is supported by elliptical particle motions (Fig. S7). The coda of cross-correlation is usually interpreted in terms of surface waves (Table 2). The main difference between the auto- and cross-correlation coda is receiver separation, which may not suffice to change from observing predominantly scattered surface waves to predominantly scattered body waves. Moreover, Yuan et al. (2021) found in numerical experiments that surface waves tend to dominate single-station correlation sensitivity in scattering media in the presence of depth-varying velocity structure, which is the case in sedimentary basins. For these reasons, we consider an interpretation in terms of surface waves to be sensible. Apart from the study by Yuan et al. (2021), the single-station correlation sensitivity to velocity changes remains scarcely researched in comparison to cross-correlation coda sensitivity. This topic merits further attention.

Assuming that surface waves dominate our observations, we deduce that the variation of parameter values with frequency range is an effect of surface wave dispersion. As a working hypothesis for the differences in terms of parameters between components, we assume that waves recorded on E, N and Z components are dominated by different scattered arrivals or different wave modes and therefore are sensitive to different parts of the surrounding medium.

Figure 6Inferred sensitivity of vs to thermal stress at seismic stations (Fig. 1). Colored triangles show the results for the North, Vertical and East components at each station, where the inverted triangle indicates the station location and shows the value of the Vertical component. White triangles show excluded models that had a misfit above the threshold. Gray triangles show the excluded observed data that did not pass quality control. Symbols for station AOVM are shown further east than the station to accommodate them on the map. The shaded relief is from SRTM (Farr et al.2007; NASA JPL2013), and the geotechnical zonation is from Gobierno de la Ciudad de México (2017).

5.3.1 Sensitivity of vs to thermal stress

Inverted values for the depth-independent sensitivity st of vs to thermal stress range from 10−6 to 0.1K−1 (all sites) and 10−6 to 0.01K−1 (only basin sites; Fig. 6). Inside the basin, st ranges from 10−6 to 0.01K−1. For the frequencies 1–2 and 2–4 Hz, we observe that stations with a high sensitivity to surface temperature variations are mostly located in and near the lake zone, whereas lower st is found mostly in the transition zone. This split disappears at 4–8 Hz.

When considering Eq. (3) and assuming that Poisson's ratio ν ranges approximately from 0.2 to 0.5, thermal expansion α from 10−6 to 10−5K−1 and elastic non-linearity ρv2σc from 10 to 1000 (Richter et al.2014), a physically reasonable range of st is 3×10-5 to 6×10-2K−1. Most inferred values fall into this range. However, at hard sites (i.e., stations outside the transition zone outlined in yellow on the map in Fig. 1), particularly at those at higher elevation, st estimates are beyond what appears to be physically plausible. As discussed in Sect. 5.2, we believe that the lack of knowledge of subsurface structure at these stations leads to inaccurate depth sensitivity kernels. In particular, near-surface sensitivity may be underestimated, as no near-surface low-velocity layer is included in hard site profiles, which in turn could lead to an overestimate of the temperature sensitivity.

A second possible reason is the fact that we neglected the depth term of Berger's thermoelastic solution. If we consider that surface wave sensitivity is greater at depths for sites where the medium does not have a thick low-velocity sediment cover, then it is plausible that this term is more important for bedrock sites than for soft sediment sites.

Figure 7Inferred minimum effective pressure p0 at seismic stations (Fig. 1). Colored triangles show the results for the North, Vertical and East components at each station, where the inverted triangle indicates the station location and shows the value of the Vertical component. White triangles show excluded models that had a misfit above the threshold. Gray triangles show the excluded observed data that did not pass quality control. Symbols for station AOVM are shown further east than the station to accommodate them on the map. The shaded relief is from SRTM (Farr et al.2007; NASA JPL2013), and the geotechnical zonation is from Gobierno de la Ciudad de México (2017). Low values of p0 correspond to a high sensitivity of vs to pore pressure changes – see Eq. (4).

5.3.2 Minimum effective pressure p0

Inverted effective pressures are mostly in the range of 0.01 Pa to 10 kPa. While the spatial pattern of the results varies, the stations with the highest minimum effective pressures are mostly found in and near the lake zone (Fig. 7).

Our observations may not constrain the details of the shallowest sediment, as surface waves may not be sensitive to the topmost layers at all sites and frequency ranges. Confining pressure may reach about 30 to 50 kPa in the first 20 m. Thus, inverted values of p0 of up to 10 kPa are in a reasonable range.

Our current model of poroelastic stress is based on a homogeneous half-space solution (Roeloffs1988), whereas the central Mexico City basin is known to have a complex hydrologic structure with an aquitard and several underlying aquifers between interbedding of lacustrine sediments and tephra deposits (Arce et al.2019). In addition, the built environment strongly influences whether water can enter the sediment; the surface rainwater is retained at the surface or redirected in pipes. A more refined model of hydrological dcc could be constructed if the measured hydraulic heads of groundwater wells at high temporal resolution were available. The results of p0 are nevertheless informative. Low values of p0 (i.e., a high value on the righthand side of Eq. 4) suggest that vs is sensitive to precipitation at a site, while a high p0 indicates the opposite.

Lower effective pressures in the hill and transition zones show that vs is sensitive to precipitation there. This may be because these sites are located where the volcanic and alluvial-pyroclastic aquifers are close to the surface (Vargas and Ortega-Guerrero2004). Two stations in the lake zone (MULU and VRVM – see Fig. 1) also show low p0 at a lower frequency of 0.5–1 Hz. A possible interpretation is that the lower-frequency observations possess some sensitivity to the shallow aquifer, which lies at a depth of approximately 50 m at those sites, while observations at 2–8 Hz are mostly sensitive to the lacustrine clay of the overlying impermeable aquitard.

Assi Hagmaier (2022) recently identified the presence of poroelastic seasonal velocity variations in the Valley of Mexico. Her analysis shows that these are too small to visibly influence the resonance period identified by horizontal-to-vertical spectral ratios (HVSR). However, she noted that an additional resonance peak dominates in the rainy season at lake zone sites but disappears in the dry season. These findings stress that it is important to study seasonal effects on both dvv and HVSR, especially considering that the latter is commonly used for site effect assessment.

Figure 8Inverted relative velocity drop and τmax for the 2017 Puebla earthquake. (a)(c) show the relative velocity drop for the E, N and Z components. (d)(f) show the maximum recovery timescale τmax for the E, N and Z components. The gray-shaded area highlights inverted recovery timescales that surpass the recording duration. For both panels, error bars show the range between the 16th and 84th percentiles.


5.4 Co-seismic damage and recovery

We present the co-seismic drop and maximum relaxation timescale τmax of the 2017 Puebla earthquake with respect to peak ground acceleration (PGA) in Fig. 8. Markers show median models of velocity drop and τmax, and error bars show the range between the 16th and 84th percentiles. Not all strong-motion stations of the RSVM network recorded the Puebla earthquake; where PGA is not available from the RSVM stations, we use the PGA at the nearest available station (taken from horizontal ground motions), including the triggered strong-motion stations shown by black rhombs in (Fig. 1). We add a dither (random shift) of up to 0.025 m s−2 to the PGA values in order to ensure that all markers are visible. Gray-shaded rectangles indicate values of τmax that lie beyond the duration of observation (approximately 2.5 years post-earthquake). Velocity drops strongly by up to ≈10 %, and recovers slowly at most sites, with several τmax on the order of a decade. The base 10 logarithms of both τmax and the velocity drop show significant positive correlation with the logarithm of the PGA (p<0.05), but the observed variances are not well explained (R2=0.28 for τmax and R2=0.12 for the velocity drop). We note that both τmax and, to a lesser extent, the velocity drop are not very well constrained (wide bars).

Due to limited data coverage of the 2020 Oaxaca earthquake, we will not interpret results of its co-seismic velocity drop and recovery apart from stating that it caused a sudden drop in phase velocity at several sites in Mexico City which was smaller than the drop during the 2017 Puebla earthquake.

A velocity decrease is observed across all the studied locations, including those on lacustrine clay, supporting the conclusions by Singh et al. (1988b) and Beresnev and Wen (1996), who stated that the lacustrine sediment behaves non-linearly.

Slow dynamics, the approximately logarithmic recovery of material mechanical properties after sudden changes induced by transient strains, is a subject of active research (TenCate2011; Sens-Schönfelder and Eulenfeld2019; Ostrovsky et al.2019, e.g.,). However, it is difficult to test current hypotheses in the metropolitan and geologically complex area of Mexico City given our sparse seismic network. To the extent that interpretation is possible, our results suggest that larger perturbations (in terms of higher PGA) may lead to a longer recovery time τmax. This is consistent with other in situ observations, notably of Viens et al. (2018), who observed a slower recovery of velocities in the Kanto basin after the 2011 Tohoku-Oki earthquake at sites that experienced a higher strain rate during the mainshock. Illien et al. (2022) applied the Snieder et al. (2017) model to the 2015 Gorkha earthquake and aftershocks and found a τmax of 846 d for the main shock, which is much longer than for the aftershocks (155 d). These observations run counter to laboratory experiments of slow dynamics, which suggest that the recovery of a particular material is independent of the perturbation amplitude (Shokouhi et al.2017). The τmax on the order of 10 years inferred from our study are uncharacteristically long compared to other studies, which reported 100 d for the Kanto basin (Viens et al.2018), 250 d in Nepal (Illien et al.2022), and 3 years in Parkfield (Brenguier et al.2008a; Wu et al.2016).

A possible explanation for the inferred long τmax is that the relaxation functions used to model dcc do not account for permanent damage. Permanent changes in the interferometric waveforms may lead to poor estimates of the relaxation timescale. Permanent damage can be assessed using the decorrelation (Larose et al.2010)

(6) Dcorr = 1 - CC i

in Fig. 9. Here, CCi is the correlation coefficient of the stretched traces with respect to the average trace for 2017. Stations that were not operational during the 2017 Puebla earthquake are omitted, and data with poor overall CCs are excluded (amounting to 40 % of data at 0.5 Hz, 12 % at 1 Hz, 3 % at 2 Hz and 6 % at 4 Hz). In the frequency bands 1–2, 2–4 and 4–8 Hz, decorrelation increases after the time window containing the earthquake (gray rectangle). This increase is more pronounced for stations not located on the clay aquitard, as defined in Solano-Rojas et al. (2015), although those stations experienced, on average, a lower PGA. At 0.5–1 Hz, decorrelation is less marked and is similar between stations on and off the aquitard. Using decorrelation as a proxy for permanent damage, this suggests that, compared to other sediments, the shallow clay of the lake zone would suffer less permanent damage.

Figure 9Observed decorrelation for stations located at sites with (basin) and without (hill) lacustrine sediments. The gray rectangle indicates the timing of the 2017 Puebla earthquake (time resolution of our stacks 20 d).


Inferred seismic velocity drops in our model are high and weakly, but significantly, correlated in double-logarithmic space with PGA, as well as with the strain rate proxy PGA/vs10. We use the first 10 m instead of the first 30 m here, as this distinguishes more accurately between hard sites and intermediate and soft sites. A dependence of the magnitude of the velocity change on peak ground motion amplitudes is expected based on laboratory studies (Shokouhi et al.2017) and is consistent with earlier field studies (Viens et al.2018).

Figure 10Slope of the linear term in the dcc model versus the downward vertical displacement measured by InSAR. The moderate correlation suggests that the seismic velocity increase in Mexico City may be caused by compaction associated with groundwater extraction.


5.5 Residual linear trend

Following the visual appearance of the dcc time series, we introduced a linear term in the dcc model (see Sect. 4). The inversion results in strong slopes for this term, with phase velocity increases of up to 0.75 % per year (see Fig. 10). A linear term in a dvv model was previously used by Gassenmeier et al. (2016), who similarly introduced it on a heuristic basis as the data appeared to require it. Slope values were comparable to those we observe (0.27 % per year) and were hypothesized to relate to recovery from an earthquake that had occurred 5 years prior to the start of their observations, with modified Mercalli intensity (MMI) IV to V at the investigated site. Prior to the 2017 Puebla earthquake, the most recent earthquake that inflicted severe damage on Mexico City was the 1985 Michoacán earthquake (Singh et al.1988a), which caused an MMI of IX in Mexico City (Arciniega-Ceballos et al.2018). Despite its large intensity, we find it unlikely that the subsurface would still be recovering 30 years later (in a logarithmic regime, a recovery rate of 0.5 % per year would require an initial velocity drop of at least 50 %, which is unrealistic). However, we cannot entirely rule out the possibility that the subsurface is recovering from the cumulative effect of multiple events.

We propose an alternative hypothesis. Mexico City has been undergoing rapid subsidence since more than a century ago due to groundwater extraction and sediment compaction (Solano-Rojas et al.2015; Chaussard et al.2021; Cigna and Tapete2021; Cabral-Cano et al.2008, and references therein). It has been suggested that another consequence of the subsidence process is the reduction of the fundamental resonance periods of the sediment strata throughout the basin (Ovando-Shelley et al.2007; Arroyo et al.2013), which may be caused by an increase in seismic velocity, as well as the compaction and resulting reduction of sediment strata thickness. Figure 10 shows the correlation between the slope of the linear trend of our models and the vertical displacement from InSAR, which we interpret as the effects of ground subsidence. Although the datasets cover different time ranges, subsidence rates have been approximately constant for decades (Cabral-Cano et al.2008; Chaussard et al.2021; Osmanoglu et al.2011) so that the rates can be directly compared. Both rates are moderately and significantly correlated (R2=0.29, p<0.0001), indicating that sediment compaction may indeed be causing a velocity increase in Mexico City's underlying stratigraphy. Similar results from Taira et al. (2018) at the Salton Sea (California) showed a steady velocity increase at a much smaller rate of less than 0.1 % yr−1, which they interpreted as an effect of poroelastic contraction. Because of the importance of resonance frequency changes as estimated by Arroyo et al. (2013) for Mexico City's seismic hazard, the ongoing rapid subsidence and its associated hazard (Fernández-Torres et al.2022), and the magnitude of the changes, this topic merits further and more detailed investigation.

6 Conclusion

We presented a comprehensive study of seismic velocity changes underneath Mexico City. Our study has several innovative aspects. (i) We used clustered autocorrelations of urban noise recorded by a strong-motion network. (ii) We modeled array-wide velocity change time series using a linear superposition of mostly physics-based terms, namely exponential relaxation for slow dynamics, poroelastic changes, thermoelastic changes and a heuristic (not physics-based) linear trend. (iii) We conducted a probabilistic inversion for the unknown model parameters.

We find that autocorrelations at strong-motion stations can be used for coda wave monitoring, at least in urban high-noise settings where results are comparable between a strong-motion and a co-located broadband sensor.

We estimated that observed velocity changes for frequencies above 2 Hz at soft sites atop very low-vs lacustrine sediments and at intermediate sites are mostly related to shear wave velocity changes in the top 100 m, relevant to site effects. Our model performs best in this region of the array, likely due to the larger amount of prior knowledge on the shallow subsurface structure there. Observed seasonal velocity changes in this region and at these depths reach 1 % peak-to-peak amplitude. At most sites, observed seasonal velocity changes show clear differences between the East, North and Vertical components.

We showed that poroelastic and thermoelastic effects can be modeled in a self-consistent manner with physically reasonable results. Stations on thick lacustrine deposits show greater sensitivity to surface temperature than stations on shallow lacustrine deposits overlying the alluvial-pyroclastic aquifer. Sensitivity to precipitation appears to be greater for sites outside of the former lake perimeter on volcanic or alluvial-pyroclastic aquifers, while it is low atop the thick lacustrine deposit. Future research should investigate the spatial sensitivity of different component autocorrelations.

Velocity drops during the Puebla 2017 and the Oaxaca 2020 earthquakes, followed by logarithmic recovery, indicate that sediments throughout the array show non-linear elastic behavior, with transient strong velocity changes on the order of 1 %–10 %. We conclude from stronger de-correlation at hard sites that permanent damage is more common there than on the soft lacustrine sediment sites. Future studies modeling slow dynamics using field measurements should account for permanent damage (e.g., through the added parameter of a static offset after the earthquake) and should investigate the cumulative effect of multiple earthquakes with longer-duration observations.

Finally, we observe an increasing trend in surface wave phase velocity that is positively and significantly correlated with vertical displacements from InSAR, while InSAR measurements show the signal of compaction of the aquitard and aquifer due to groundwater extraction. As this trend is strong (up to approximately 1 % yr−1) and may provide additional depth-sensitive information about the compaction processes, it certainly merits further investigation, particularly in light of previously reported resonance frequency changes.

Code and data availability

Data from the Geoscope station G.UNM are public and can be obtained through FDSN web services (; Institut de physique du globe de Paris and Ecole et Observatoire des Sciences de la Terre de Strasbourg1982). Data from the Red Sísmica del Valle de México are collected and curated by the Mexican National Seismological Service (SSN). Seismic peak ground acceleration data at stations SCT2, CUP5, CCCS and TACY were provided by the Accelerographic Network of the Institute of Engineering (RAII-UNAM) and are a product of the instrumentation, processing and distribution of the Seismic Instrumentation Unit. The data are distributed through the Accelerographic Database System: (Instituto de Ingeneria, UNAM, Grupo de Procesamiento y Analisis Sismico2021; Quaas et al.1995). Topographic data used for shaded reliefs were obtained via PyGMT (, last access: 12 May 2023) and (Uieda et al.2021; Tozer et al.2019) and are based on NASA SRTM data (Farr et al.2007; NASA JPL2013). Copernicus Sentinel-1 SAR data were retrieved from ASF DAAC, processed by ESA.

The processing and correlation code ants is available at (Ermert and Kaestle2023). The NoisePy module used for measuring stretching is available at (Jiang and Denolle2020). The ruido Python module used for data handling, urban noise clustering and stacking is available at (last access: 12 May 2023) and (Ermert2023).


The supplement related to this article is available online at:

Author contributions

Conceptualization: MAD, LAE. Data curation: ECC, LQ, LAE, EAFT, DMP. Formal analysis: LAE. Funding acquisition: MAD (supported LAE), EAFT, ECC. Methodology and investigation: LAE, MAD. Project administration: MAD, ECC. Resources: MAD, ECC, LQ. Software: LAE, MAD. Supervision: MAD. Validation: LAE, EC, ECC, DSR. Visualization: LAE, ECC. Writing – original draft: LAE. Writing – review and editing: LAE, MAD, EC, ECC, LQ, DSR.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank an anonymous reviewer and Alicia Hotovec-Ellis for their positive comments and detailed feedback. We are also very grateful to editor Michal Malinowski for handling and commenting on the paper.

This research was enabled by continuous seismic data provided by the RSVM (Red Sísmica del Valle de México) network. RSVM is maintained by the personnel of the Servicio Sismológico Nacional (SSN, Universidad Nacional Autónoma de México, Instituto de Geofísica). Financial support was provided by Consejo Nacional de Ciencia y Tecnología (CONACyT, National Council for Science and Technology) and the Government of Mexico City through agreement nos. SECTEI/194/2017, CM-SECTEI/263/2021 and CM-SECTEI/156/2022. Seismic peak ground acceleration data at stations SCT2, CUP5, CCCS and TACY were provided by the Accelerographic Network of the Institute of Engineering (RAII-UNAM) and are a product of the instrumentation, processing and distribution of the Seismic Instrumentation Unit. We gratefully acknowledge the data collected by Geoscope and the Tectonic Observatory, which are archived on IRIS. Laura A. Ermert and Marine A. Denolle have been supported by MD's fellowship from the David and Lucile Packard Foundation. Enrique Cabral Cano has been supported by DGAPA-PAPIIT project no. IN107321. Enrique A. Fernández Torres is supported by a fellowship from CONACyT, México, under grant agreement no. CVU 863837. Darío Solano Rojas and Diana Morales Padilla were supported by grant nos. UNAM-PAPIIT IA105921 and LANCAD-UNAM-DGTIC-380. InSAR processing was performed at UNAM's Miztli supercomputer under grant no. LANCAD-UNAM-DGTIC-362 to Enrique Cabral Cano. Laura A. Ermert thanks the Pacific Northwest Seismic Network for supporting this research with caffeine and Mouse Reusch and Paul Bodin for sharing their insights on telemetry gaps and the site conditions in the lake zone. We also acknowledge the valuable discussions with Tim Clements about porolastic effects on seismic velocity, as well as with Naiara Korta Martiartu and Kurama Okubo about probabilistic inversions.

Financial support

This research has been supported by the David and Lucile Packard Foundation (Marine A. Denolle's Fellowship); the Consejo Nacional de Ciencia y Tecnología, Sistema Nacional de Investigadores (grant nos. CVU863837, SECTEI/194/2017, CM-4SECTEI/263/2021 and CM-SECTEI/156/2022), and the Universidad Nacional Autónoma de México (Dirección General de Asuntos del Personal Académico: grant nos. IN107321 and UNAM-PAPIIT IA 105921; Dirección General de Cómputo y de Tecnologías de Información y Comunicación: grant nos. LANCAD-UNAM-DGTIC-380 and LANCAD-UNAM- DGTIC-362).

Review statement

This paper was edited by Michal Malinowski and reviewed by Alicia Hotovec-Ellis and one anonymous referee.


Ajo-Franklin, J., Dou, S., Daley, T., Freifeld, B., Robertson, M., Ulrich, C., Wood, T., Eckblaw, I., Lindsey, N., Martin, E., and Wagner, A.: Time-lapse surface wave monitoring of permafrost thaw using distributed acoustic sensing and a permanent automated seismic source, in: 2017 SEG International Exposition and Annual Meeting, September 2017, Houston, Texas, OnePetro, paper number SEG-2017-17774027, 2017. a

Alberto, Y., Otsubo, M., Kyokawa, H., Kiyota, T., and Towhata, I.: Reconnaissance of the 2017 Puebla, Mexico earthquake, Soils and Foundations, 58, 1073–1092,, 2018. a

Andajani, R. D., Tsuji, T., Snieder, R., and Ikeda, T.: Spatial and temporal influence of rainfall on crustal pore pressure based on seismic velocity monitoring, Earth Planet. Space, 72, 1–17, 2020. a

Anderson, J. G., Bodin, P., Brune, J. N., Prince, J., Singh, S. K., Quaas, R., and Onate, M.: Strong Ground Motion from the Michoacan, Mexico, Earthquake, Science, 233, 1043–1049,, 1986. a

Arce, J. L., Layer, P. W., Macías, J. L., Morales-Casique, E., García-Palomo, A., Jiménez-Domínguez, F. J., Benowitz, J., and Vásquez-Serrano, A.: Geology and stratigraphy of the Mexico Basin (Mexico City), central Trans-Mexican Volcanic Belt, J. Maps, 15, 320–332,, 2019. a

Arciniega-Ceballos, A., Baena-Rivera, M., and Sánchez-Sesma, F. J.: The 1985 (M8. 1) Michoacán Earthquake and its effects in Mexico City, in: Living under the threat of earthquakes, Springer, 65–75, 2018. a, b

Arroyo, D., Ordaz, M., Ovando-Shelley, E., Guasch, J. C., Lermo, J., Perez, C., Alcantara, L., and Ramírez-Centeno, M. S.: Evaluation of the change in dominant periods in the lake-bed zone of Mexico City produced by ground subsidence through the use of site amplification factors, Soil Dyn. Earthq. Eng., 44, 54–66, 2013. a, b

Asimaki, D., Mohammadi, K., Ayoubi, P., Mayoral, J. M., and Montalva, G.: Investigating the spatial variability of ground motions during the 2017 Mw 7.1 Puebla-Mexico City earthquake via idealized simulations of basin effects, Soil Dyn. Earthq. Eng., 132, 106073,, 2020. a

Assi Hagmaier, H.: Seasonal variations of seismic velocities and vibration periods in the Valley of Mexico observed with seismic noise processing, Master's thesis, Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE), 2022. a

Beresnev, I. A. and Wen, K.-L.: Nonlinear soil response—A reality?, B. Seismol. Soc. Am., 86, 1964–1978, 1996. a

Bergamo, P., Dashwood, B., Uhlemann, S., Swift, R., Chambers, J. E., Gunn, D. A., and Donohue, S.: Time-lapse monitoring of climate effects on earthworks using surface wavesTime-lapse seismic monitoring with SW, Geophysics, 81, EN1–EN15, 2016. a

Berger, J.: A note on thermoelastic strains and tilts, J. Geophys. Res., 80, 274–277, 1975. a

Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., and Wassermann, J.: ObsPy: A Python Toolbox for Seismology, Seismol. Res. Lett., 81, 530–533,, 2010. a

Bondár, I. and Storchak, D.: Improved location procedures at the International Seismological Centre, Geophys. J. Int., 186, 1220–1244,, 2011. a

Bonilla, L. F. and Ben-Zion, Y.: Detailed space–time variations of the seismic response of the shallow crust to small earthquakes from analysis of dense array data, Geophys. J. Int., 225, 298–310,, 2020. a

Bonilla, L. F., Guéguen, P., and Ben‐Zion, Y.: Monitoring Coseismic Temporal Changes of Shallow Material during Strong Ground Motion with Interferometry and Autocorrelation, Bull. Seismol. Soc. Am., 109, 187–198,, 2019. a, b, c, d

Brenguier, F., Campillo, M., Hadziioannou, C., Shapiro, N. M., Nadeau, R. M., and Larose, É.: Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations, Science, 321, 1478–1481, 2008a. a

Brenguier, F., Shapiro, N. M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O., and Nercessian, A.: Towards forecasting volcanic eruptions using seismic noise, Nat. Geosci., 1, 126–130, 2008b. a

Cabral-Cano, E., Dixon, T. H., Miralles-Wilhelm, F., Díaz-Molina, O., Sánchez-Zamora, O., and Carande, R. E.: Space geodetic imaging of rapid ground subsidence in Mexico City, Geol. Soc. Am. B., 120, 1556–1566, 2008. a, b, c

Chaussard, E., Havazli, E., Fattahi, H., Cabral-Cano, E., and Solano-Rojas, D.: Over a Century of Sinking in Mexico City: No Hope for Significant Elevation and Storage Capacity Recovery, J. Geophys. Res., 126, e2020JB020648,, 2021. a, b, c

Chavez-Garcia, F. J. and Bard, P.-Y.: Site effects in Mexico City eight years after the September 1985 Michoacan earthquakes, Soil Dyn. Earthq. Eng., 13, 229–247, 1994. a

Cigna, F. and Tapete, D.: Present-day land subsidence rates, surface faulting hazard and risk in Mexico City with 2014–2020 Sentinel-1 IW InSAR, Remote Sens. Environ., 253, 112161,, 2021. a

Clements, T. and Denolle, M. A.: Tracking groundwater levels using the ambient seismic field, Geophys. Res. Lett., 45, 6459–6465, 2018. a

Crameri, F., Shephard, G. E., and Heron, P. J.: The misuse of colour in science communication, Nat. Commun., 11, 1–10, 2020. a

Cruz-Atienza, V. M., Tago, J., Sanabria-Gómez, J., Chaljub, E., Etienne, V., Virieux, J., and Quintanar, L.: Long duration of ground motion in the paradigmatic valley of Mexico, Sci. Rep.-UK, 6, 1–9, 2016. a, b, c, d

De Plaen, R. S., Cannata, A., Cannavo, F., Caudron, C., Lecocq, T., and Francis, O.: Temporal changes of seismic velocity caused by volcanic activity at Mt. Etna revealed by the autocorrelation of ambient seismic noise, Front. Earth Sci., 6, 251,, 2019. a

Donaldson, C., Winder, T., Caudron, C., and White, R. S.: Crustal seismic velocity responds to a magmatic intrusion and seasonal loading in Iceland’s Northern Volcanic Zone, Sci. Adv., 5, eaax6642,, 2019. a, b, c, d

Ekström, G.: Time domain analysis of Earth's long-period background seismic radiation, J. Geophys. Res., 106, 26483–26493,, 2001. a

Ermert, L.: lermert/ruido: release 0: testing (v0.0.0-alpha), Zenodo [code],, 2023. a

Ermert, L. and Kaestle, E.: ants_2 – a lightweight ambient noise processing tool (Version prelim), Zenodo [software],, 2023. a

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D.: The shuttle radar topography mission, Rev. Geophys., 45, 2005RG000183,, 2007. a, b, c, d

Feng, K.-F., Huang, H.-H., Hsu, Y.-J., and Wu, Y.-M.: Controls on Seasonal Variations of Crustal Seismic Velocity in Taiwan Using Single-Station Cross-Component Analysis of Ambient Noise Interferometry, J. Geophys. Res., 126, e2021JB022650,, 2021. a, b, c

Fernández-Torres, E. A., Cabral-Cano, E., Novelo-Casanova, D. A., Solano-Rojas, D., Havazli, E., and Salazar-Tlaczani, L.: Risk assessment of land subsidence and associated faulting in Mexico City using InSAR, Nat. Hazards, 112, 37–55, 2022. a

Fichtner, A.: surf, GitHub [code],, 2020. a

Field, E. H., T. S. P. I. W. G.: Accounting for Site Effects in Probabilistic Seismic Hazard Analyses of Southern California: Overview of the SCEC Phase III Report, Bull. Seismol. Soc. Am., 90, S1–S31,, 2000. a

Fokker, E., Ruigrok, E., Hawkins, R., and Trampert, J.: Physics-Based Relationship for Pore Pressure and Vertical Stress Monitoring Using Seismic Velocity Variations, Remote Sens., 13, 2684,, 2021. a, b, c

Folk, M., Heber, G., Koziol, Q., Pourmal, E., and Robinson, D.: An overview of the HDF5 technology suite and its applications, in: Proceedings of the EDBT/ICDT 2011 Workshop on Array Databases, 36–47, ACM, 2011. a

Foreman-Mackey, D., Hogg, D. W., Lang, D., and Goodman, J.: emcee: the MCMC hammer, Publ. Astron. Soc. Pac., 125, 306,, 2013. a

Foti, S., Aimar, M., Ciancimino, A., and Passeri, F.: Recent developments in seismic site response evaluation and microzonation, in: Proceedings of the XVII ECSMGE, 223–248, European Conference on soil mechanics and geotechnical engineering, September 2019, Reykjavik, Iceland,, 2019. a

Froment, B., Campillo, M., Chen, J., and Liu, Q.: Deformation at depth associated with the 12 May 2008 Mw 7.9 Wenchuan earthquake from seismic ambient noise monitoring, Geophys. Res. Lett., 40, 78–82, 2013. a

Gassenmeier, M., Sens-Schönfelder, C., Eulenfeld, T., Bartsch, M., Victor, P., Tilmann, F., and Korn, M.: Field observations of seismic velocity changes caused by shaking-induced damage and healing due to mesoscopic nonlinearity, Geophys. J. Int., 204, 1490–1502, 2016. a, b

Gobierno de la Ciudad de México: Normas técnicas complementarias para diseño y construcción de cimentaciones, Gaceta uficial de la Ciudad de México, 10–43, 15 December 2017. a, b, c

Groos, J. and Ritter, J.: Time domain classification and quantification of seismic noise in an urban environment, Geophys. J. Int., 179, 1213–1231, 2009. a

Hadziioannou, C., Larose, E., Baig, A., Roux, P., and Campillo, M.: Improving temporal resolution in ambient noise monitoring of seismic wave speed, J. Geophys. Res., 116, B07304,, 2011. a

Hobiger, M., Wegler, U., Shiomi, K., and Nakahara, H.: Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake, Geophys. J. Int., 198, 90–109,, 2014. a, b

Hobiger, M., Wegler, U., Shiomi, K., and Nakahara, H.: Coseismic and post-seismic velocity changes detected by passive image interferometry: comparison of one great and five strong earthquakes in Japan, Geophys. J. Int., 205, 1053–1073, 2016. a, b

Illien, L., Andermann, C., Sens-Schönfelder, C., Cook, K., Baidya, K., Adhikari, L., and Hovius, N.: Subsurface moisture regulates Himalayan groundwater storage and discharge, AGU Advances, 2, e2021AV000398,, 2021. a, b

Illien, L., Sens-Schönfelder, C., Andermann, C., Marc, O., Cook, K. L., Adhikari, L. B., and Hovius, N.: Seismic velocity recovery in the subsurface: transient damage and groundwater drainage following the 2015 Gorkha earthquake, Nepal, J. Geophys. Res., 127, e2021JB023402,, 2022. a, b, c, d

Institut de physique du globe de Paris (IPGP) and Ecole et Observatoire des Sciences de la Terre de Strasbourg (EOST): GEOSCOPE, French global network of broad band seismic stations, Institut de physique du globe de Paris (IPGP), Université de Paris [data set],, 1982. a

Instituto de Ciencias de la Atmósfera y Cambio Climatico: Programa de Estaciones Meteorológicas del Bachillerato Universitario, (last access: 13 May 2023), 2022. a

Instituto de Ingeneria, UNAM, Grupo de Procesamiento y Analisis Sismico: Archivo Estandar de Aceleracion [data set]., last access: 21 July 2021. a

International Seismological Centre: On-line Bulletin, (last access: 15 April 2022), 2022. a

James, S., Knox, H., Abbott, R., Panning, M., and Screaton, E.: Insights into permafrost and seasonal active-layer dynamics from ambient seismic noise monitoring, J. Geophys. Res., 124, 1798–1816, 2019. a

Jiang, C. and Denolle, M.: NoisePy: a new high-performance python tool for seismic ambient noise seismology, Seismol. Res. Lett., 91, 1853–1866,, 2020 (code available at:, last access: 12 May 2023). a

Kasahara, K., Sakai, S., Morita, Y., Hirata, N., Tsuruoka, H., Nakagawa, S., Nanjo, K. Z., and Obara, K.: Development of the Metropolitan Seismic Observation Network (MeSO-net) for detection of mega-thrust beneath Tokyo metropolitan area, Bull. Earthq. Res. Inst. Univ. Tokyo, 84, 71–88, 2009. a

Kennet, B.: IASPEI 1991 seismological tables, Terra Nova, 3, 122–122, 1991. a

Krischer, L., Megies, T., Barsch, R., Beyreuther, M., Lecocq, T., Caudron, C., and Wassermann, J.: ObsPy: a bridge for seismology into the scientific Python ecosystem, Computational Science & Discovery, 8, 014003,, 2015. a

Larose, E., Planes, T., Rossetto, V., and Margerin, L.: Locating a small change in a multiple scattering environment, Appl. Phys. Lett., 96, 204101,, 2010. a

Lecocq, T., Longuevergne, L., Pedersen, H. A., Brenguier, F., and Stammler, K.: Monitoring ground water storage at mesoscale using seismic noise: 30 years of continuous observation and thermo-elastic and hydrological modeling, Sci. Rep.-UK, 7, 1–16, 2017. a, b, c

Lindner, F., Wassermann, J., and Igel, H.: Seasonal Freeze-Thaw Cycles and Permafrost Degradation on Mt. Zugspitze (German/Austrian Alps) Revealed by Single-Station Seismic Monitoring, Geophys. Res. Lett., 48, e2021GL094659,, 2021. a

Mao, S., Lecointre, A., van der Hilst, R. D., and Campillo, M.: Space-time monitoring of groundwater fluctuations with passive seismic interferometry, Nat. Commun., 13, 4643,, 2022. a

MASE: Meso America Subduction Experiment,, 2007. a, b

Mavko, G., Mukerji, T., and Dvorkin, J. (Eds.): The rock physics handbook, Cambridge University Press,, 2020. a

Mordret, A., Courbis, R., Brenguier, F., Chmiel, M., Garambois, S., Mao, S., Boué, P., Campman, X., Lecocq, T., Van der Veen, W., and Hollis, D.: Noise-based ballistic wave passive seismic monitoring–Part 2: surface waves, Geophys. J. Int., 221, 692–705, 2020. a

Nakata, N. and Snieder, R.: Near-surface weakening in Japan after the 2011 Tohoku-Oki earthquake, Geophys. Res. Lett., 38, L17302,, 2011. a

NASA JPL: NASA Shuttle Radar Topography Mission Global 1 arc second, NASA EOSDIS Land Processes DAAC [data set],, 2013. a, b, c, d

Obermann, A. and Hillers, G.: Seismic time-lapse interferometry across scales, Adv. Geophys., 60, 65–143, 2019. a

Obermann, A., Froment, B., Campillo, M., Larose, E., Planès, T., Valette, B., Chen, J., and Liu, Q.: Seismic noise correlations to image structural and mechanical changes associated with the Mw 7.9 2008 Wenchuan earthquake, J. Geophys. Res., 119, 3155–3168, 2014. a, b

Oral, E., Gélis, C., and Bonilla, L. F.: 2-D P-SV and SH spectral element modelling of seismic wave propagation in non-linear media with pore-pressure effects, Geophys. J. Int., 217, 1353–1365, 2019. a

Ortega, A. and Farvolden, R. N.: Computer analysis of regional groundwater flow and boundary conditions in the Basin of Mexico, J. Hydrol., 110, 271–294, 1989. a, b

Osmanoglu, B., Dixon, T. H., Wdowinski, S., Cabral-Cano, E., and Jiang, Y.: Mexico City subsidence observed with persistent scatterer InSAR, Int. J. Appl. Earth Obs., 13, 1–12,, 2011. a

Ostrovsky, L., Lebedev, A., Riviere, J., Shokouhi, P., Wu, C., Stuber Geesey, M. A., and Johnson, P. A.: Long-time relaxation induced by dynamic forcing in geomaterials, J. Geophys. Res., 124, 5003–5013, 2019. a

Ovando-Shelley, E., Ossa, A., and Romo, M. P.: The sinking of Mexico City: Its effects on soil properties and seismic response, Soil Dyn. Earthq. Eng., 27, 333–343, 2007. a

Pérez-Campos, X., Espíndola, V. H., González-Ávila, D., Zanolli Fabila, B., Márquez-Ramírez, V. H., De Plaen, R. S. M., Montalvo-Arrieta, J. C., and Quintanar, L.: The effect of confinement due to COVID-19 on seismic noise in Mexico, Solid Earth, 12, 1411–1419,, 2021. a

Pérez Cruz, G.: Estudio sismológico de reflexíon del subsuelo de la ciudad de México, Master's thesis, U.N.A.M, Faculty of Engineering, 1988. a, b

Peterson, J.: Observations and modeling of seismic background noise, USGS Open File Report, 93-322, 94 pp., 1993. a

Quaas, R., Alcántara, L., Espinosa, J. M., Mena, E., Otero, J. A., Javier, C., López, B., Perez, C., Vázquez, R., Flores, J. A., and González, F.: Base Mexicana de datos de sismos fuertes. Un Sistema que integra la información acelerográfica registrada en México en los últimos 35 años, Revista de Ingenería Sísmica, No. 51, 1–12, 1995. a

Quintanar, L., Cárdenas‐Ramírez, A., Bello‐Segura, D. I., Espíndola, V. H., Pérez‐Santana, J. A., Cárdenas‐Monroy, C., Carmona‐Gallegos, A. L., and Rodríguez‐Rasilla, I.: A Seismic Network for the Valley of Mexico: Present Status and Perspectives, Seismol. Res. Lett., 89, 356–362,, 2018. a, b, c, d

Richter, T., Sens-Schönfelder, C., Kind, R., and Asch, G.: Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry, J. Geophys. Res., 119, 4747–4765, 2014. a, b, c, d, e, f, g, h, i

Rodríguez Tribaldos, V. and Ajo-Franklin, J. B.: Aquifer monitoring using ambient seismic noise recorded with distributed acoustic sensing (DAS) deployed on dark fiber, J. Geophys. Res., 126, e2020JB021004,, 2021. a, b, c, d

Roeloffs, E.: Poroelastic techniques in the study of earthquake-related hydrologic phenomena, Adv. Geophys., 37, 135–195, 1996. a, b, c

Roeloffs, E. A.: Fault stability changes induced beneath a reservoir with cyclic variations in water level, J. Geophys. Res.-Sol. Ea., 93, 2107–2124, 1988. a, b, c, d, e, f

Roten, D., Fäh, D., Olsen, K. B., and Giardini, D.: A comparison of observed and simulated site response in the Rhône valley, Geophys. J. Int., 173, 958–978,, 2008. a

Roult, G., Montagner, J.-P., Romanowicz, B., Cara, M., Rouland, D., Pillet, R., Karczewski, J.-F., Rivera, L., Stutzmann, E., and Maggi, A.: The GEOSCOPE Program: Progress and challenges during the past 30 years, Seismol. Res. Lett., 81, 427–452, 2010. a

Sahakian, V. J., Melgar, D., Quintanar, L., Ramírez-Guzmán, L., Pérez-Campos, X., and Baltay, A.: Ground motions from the 7 and 19 September 2017 Tehuantepec and Puebla-Morelos, Mexico, earthquakes, Bull. Seismol. Soc. Am., 108, 3300–3312, 2018. a

Sánchez-Pastor, P., Obermann, A., and Schimmel, M.: Detecting and locating precursory signals during the 2011 El Hierro, Canary Islands, submarine eruption, Geophys. Res. Lett., 45, 10–288, 2018. a, b

Sánchez-Pastor, P., Obermann, A., Schimmel, M., Weemstra, C., Verdel, A., and Jousset, P.: Short-and long-term variations in the Reykjanes geothermal reservoir from seismic noise interferometry, Geophys. Res. Lett., 46, 5788–5798, 2019. a

Saul, M. and Lumley, D.: A new velocity–pressure–compaction model for uncemented sediments, Geophys. J. Int., 193, 905–913, 2013. a

Sawazaki, K., Saito, T., Ueno, T., and Shiomi, K.: Estimation of seismic velocity changes at different depths associated with the 2014 Northern Nagano Prefecture earthquake, Japan (MW 6.2) by joint interferometric analysis of NIED Hi-net and KiK-net records, Prog. Earth Planet. Sci., 3, 1–15, 2016. a

Schimmel, M., Stutzmann, E., and Gallart, J.: Using instantaneous phase coherence for signal extraction from ambient noise data at a local to a global scale, Geophys. J. Int., 184, 494–506, 2011. a

Schippkus, S., Garden, M., and Bokelmann, G.: Characteristics of the ambient seismic field on a large-N seismic array in the Vienna basin, Seismol. Soc. Am., 91, 2803–2816, 2020. a

Sens-Schönfelder, C. and Eulenfeld, T.: Probing the in situ elastic nonlinearity of rocks with earth tides and seismic noise, Phys. Rev. Lett., 122, 138501,, 2019. a, b

Sens-Schönfelder, C. and Wegler, U.: Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia, Geophys. Res. Lett., 33, l21302,, 2006. a, b, c

Sens-Schönfelder, C., Pomponi, E., and Peltier, A.: Dynamics of Piton de la Fournaise volcano observed by passive image interferometry with multiple references, J. Geophys. Res., 276, 32–45, 2014. a

Sens-Schönfelder, C., Snieder, R., and Li, X.: A model for nonlinear elasticity in rocks based on friction of internal interfaces and contact aging, Geophys. J. Int., 216, 319–331,, 2018. a

Shapiro, N. M., Singh, S. K., Almora, D., and Ayala, M.: Evidence of the dominance of higher-mode surface waves in the lake-bed zone of the Valley of Mexico, Geophys. J. Int., 147, 517–527,, 2001. a, b

Shokouhi, P., Rivière, J., Guyer, R. A., and Johnson, P. A.: Slow dynamics of consolidated granular systems: Multi-scale relaxation, Appl. Phys. Lett., 111, 251604,, 2017. a, b

Singh, S., Quaas, R., Ordaz, M., Mooser, F., Almora, D., Torres, M., and Vasquez, R.: Is there truly a “hard” rock site in the Valley of Mexico?, Geophys. Res. Lett., 22, 481–484, 1995. a, b, c, d

Singh, S., Santoyo, M., Bodin, P., and Gomberg, J.: Dynamic deformations of shallow sediments in the Valley of Mexico, Part II: Single-station estimates, Bull. Seismol. Soc. Am., 87, 540–550, 1997. a, b, c, d

Singh, S. K., Lermo, J., Domínguez, T., Ordaz, M., Espinosa, J. M., Mena, E., and Quaas, R.: The Mexico Earthquake of September 19, 1985—A Study of Amplification of Seismic Waves in the Valley of Mexico with Respect to a Hill Zone Site, Earthquake Spectra, 4, 653–673,, 1988a. a, b

Singh, S. K., Mena, E. A., and Castro, R.: Some aspects of source characteristics of the 19 September 1985 Michoacan earthquake and ground motion amplification in and near Mexico City from strong motion data, Bull. Seismol. Soc. Am., 78, 451–477, 1988b. a

Snieder, R., Sens-Schönfelder, C., and Wu, R.: The time dependence of rock healing as a universal relaxation process, a tutorial, Geophys. J. Int., 208, 1–9, 2017. a, b, c

Solano-Rojas, D., Cabral-Cano, E., Hernández-Espriú, A., Wdowinski, S., DeMets, C., Salazar-Tlaczani, L., Falorni, G., and Bohane, A.: La relación de subsidencia del terreno InSAR-GPS y el abatimiento del nivel estático en pozos de la zona Metropolitana de la Ciudad de México, Boletín de la Sociedad Geológica Mexicana, 67, 273–283, 2015. a, b, c, d, e, f

Solano-Rojas, D., Cabral-Cano, E., Fernández-Torres, E., Havazli, E., Wdowinski, S., and Salazar-Tlaczani, L.: Remotely triggered subsidence acceleration in Mexico City induced by the September 2017 Mw 7.1 Puebla and the Mw 8.2 Tehuantepec September 2017 earthquakes, P. Int. Ass. Hydrol. Sci., 382, 683–687, 2020. a

Steinmann, R., Hadziioannou, C., and Larose, E.: Effect of centimetric freezing of the near subsurface on Rayleigh and Love wave velocity in ambient seismic noise correlations, Geophys. J. Int., 224, 626–636, 2021. a, b

Taira, T., Nayak, A., Brenguier, F., and Manga, M.: Monitoring reservoir response to earthquakes and fluid extraction, Salton Sea geothermal field, California, Sci. Adv., 4, e1701536,, 2018. a, b

Takano, T., Nishimura, T., and Nakahara, H.: Seismic velocity changes concentrated at the shallow structure as inferred from correlation analyses of ambient noise during volcano deformation at Izu-Oshima, Japan, J. Geophys. Res., 122, 6721–6736, 2017. a, b

Takano, T., Brenguier, F., Campillo, M., Peltier, A., and Nishimura, T.: Noise-based passive ballistic wave seismic monitoring on an active volcano, Geophys. J. Int., 220, 501–507, 2020. a

Takeuchi, H. and Saito, M.: Seismic surface waves, Methods in Computational Physics, 11, 217–295, 1972. a

Talwani, P., Chen, L., and Gahalaut, K.: Seismogenic permeability, ks, J. Geophys. Res., 112, B07309,, 2007. a

TenCate, J. A.: Slow dynamics of earth materials: An experimental overview, Pure Appl. Geophys., 168, 2211–2219, 2011. a

Tozer, B., Sandwell, D. T., Smith, W. H. F., Olson, C., Beale, J. R., and Wessel, P.: Global bathymetry and topography at 15 arc sec: SRTM15+, Earth Space Sci., 6, 1847–1864,, 2019 (data available at:, last access: 12 May 2023). a

Tsai, V. C.: A model for seasonal changes in GPS positions and seismic wave speeds due to thermoelastic and hydrologic variations, J. Geophys. Res., 116, 1–9,, 2011.  a

Uieda, L., Tian, D., Leong, W., Toney, L., Schlitzer, W., Grund, M., Newton, D., Ziebarth, M., Jones, M., and Wessel, P.: PyGMT: A Python interface for the Generic Mapping Tools (v0.9.0), Zenodo [code],, 2021. a

Vargas, C. and Ortega-Guerrero, A.: Fracture hydraulic conductivity in the Mexico City clayey aquitard: Field piezometer rising-head tests, Hydrogeol. J., 12, 336–344, 2004. a

Viens, L. and Iwata, T.: Improving the retrieval of offshore-onshore correlation functions with machine learning, J. Geophys. Res., 125, e2020JB019730,, 2020 (code available at:, last access: 15 May 2023). a, b

Viens, L., Denolle, M. A., Hirata, N., and Nakagawa, S.: Complex near-surface rheology inferred from the response of greater Tokyo to strong ground motions, J. Geophys. Res., 123, 5710–5729, 2018. a, b, c, d, e

Voisin, C., Garambois, S., Massey, C., and Brossier, R.: Seismic noise monitoring of the water table in a deep-seated, slow-moving landslide, Interpretation, 4, SJ67–SJ76, 2016. a

Wang, Q.-Y., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T., and Aoki, Y.: Seasonal crustal seismic velocity changes throughout Japan, J. Geophys. Res., 122, 7987–8002, 2017. a, b, c, d, e

Wang, Q.-Y., Campillo, M., Brenguier, F., Lecointre, A., Takeda, T., and Hashima, A.: Evidence of changes of seismic properties in the entire crust beneath Japan after the M w 9.0, 2011 Tohoku-oki earthquake, J. Geophys. Res., 124, 8924–8941, 2019. a

Wu, C. and Peng, Z.: Long-term change of site response after the Mw 9.0 Tohoku earthquake in Japan, Earth Planets Space, 64, 1259–1266, 2012. a, b

Wu, C., Delorey, A., Brenguier, F., Hadziioannou, C., Daub, E. G., and Johnson, P.: Constraining depth range of S wave velocity decrease after large earthquakes near Parkfield, California, Geophys. Res. Lett., 43, 6129–6136,, 2016. a, b

Yuan, C., Bryan, J., and Denolle, M.: Numerical comparison of time-, frequency-and wavelet-domain methods for coda wave interferometry, Geophys. J. Int., 226, 828–846, 2021. a, b

Short summary
Mexico City is built on a unique ground containing the clay-rich sediments of the ancient lake Texcoco. Continuous imperceptible shaking of these deposits by city traffic and other sources allows us to monitor changes in the subsurface seismic wave speed. Wave speed varies seasonally, likely due to temperature and rain effects; it temporarily drops after large earthquakes then starts to recover. Throughout the studied period, it increased on average, which may be related to soil compaction.