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

. 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 broad-band seismic stations throughout Mexico City, including sites with different geologic characteristics ranging from city center locations situated on lacustrine clay to hillsite 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 postseismic 5 changes, as well as a long-term linear trend in seismic velocity. Seasonal variations can be explained by self-consistent models of thermo-elastic and poro-elastic 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 M w 7.1 Puebla and 2020 M w 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 10 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-15 wave monitoring given sufficiently high-amplitude urban seismic noise


Introduction
Near-surface geological structure and soil properties are important determinants of seismic hazard (e.g. Field, 2000). 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). 20 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 v s30 , 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 interfer-25 ometric measurements like cross-correlation or deconvolution of continous seismic data to a long-term reference (Obermann and Hillers, 2019). 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, which are approximately linearly related to relative velocity changes dv v . Using this technique reveals that seismic velocity varies with groundwater level (Sens-Schönfelder and Wegler, 2006;Lecocq et al., 2017;Fokker et al., 2021;Rodríguez Tribaldos and 30 Ajo-Franklin, 2021), 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 Denolle, 2018;Mao et al., 2021), as can soil compaction (Taira et al., 2018). Although the reported changes are 35 usually small, on the order of approximately 1 % peak-to-peak amplitude or 0.01 -0.1 %/year 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 Peng, 2012;Oral et al., 2019;Bonilla et al., 2019). Numerous recent studies utilizing ambient noise interferometry have reported temporary velocity drops consistent with 40 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 Peng, 2012;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 elas-45 ticity (Nakata and Snieder, 2011;Viens et al., 2018). The velocity reduction is reported to be even stronger during the shaking itself Bonilla and Ben-Zion, 2020); 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 50 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 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 55 from continuous recordings of urban seismic noise at a relatively sparse strong motion sensor network with seismic inter-ferometry. 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 Ritter, 2009;Schippkus et al., 2020, sect. 3). In sect. 4, we introduce our model for velocity 60 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).

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 65 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 seconds and shorter is well captured by both types of sensors and that results are consistent. 70 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 broad-band sensor observations from the Geoscope network site G.UNM and four stations of the temporary Tectonic Observatory deployment (MASE, 2007).
The station G.UNM includes a co-located STS-1 seismometer (broadband sensor) and Kinemetrics EpiSensor accelerometer 75 (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 Figure 1). The upper panel of Figure 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 -120dB (Peterson, 1993). Furthermore, instrumental effects are clearly visible, such as the change in instrument gain 80 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 85 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.  (Farr et al., 2007), geotechnical zonation (Gobierno de la Ciudad de México, 2017), and location of seismic stations in the Valley of Mexico region (Quintanar et al., 2018;MASE, 2007;Roult et al., 2010). Lake zone (blue outline) and transition zone (orange outline) include locations where the shallow subsurface is characterized by quaternary lacustrine sediment. Strong motion and broad-band 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 broad-band 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 Mw7.1 Puebla earthquake.
3 Measuring velocity changes in an urban setting 90 We remove large global earthquakes of M w = 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 Centre, 2022;Bondár and Storchak, 2011) are removed by cutting out a window starting ten seconds before the direct P-wave using iasp91 (Kennet, 1991) and ending 60 seconds after a 1 km s surface wave. We cut data into 8-hour segments, detrend, taper, pre-bandpass-filter them, and remove the instrument responses to correct the data to ground accel-95 eration (including the seismometer for direct comparison of the results). Finally, we correlate 20-minute windows of the data overlapping by 10 minutes. We save all single-window correlations for further processing in hdf5 format (Folk et al., 2011).
The pre-processing and correlation are performed with 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 sect. 6 for details on code availability).

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 reduce the effect of temporally varying noise sources. Clustering is performed by applying Gaussian Mixture Models after reducing the 105 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 Iwata, 2020).
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 110 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 4 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 4 clusters. It mostly coincides with the transition between day-and nighttime. The result for the North-North component of 115 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 Hz 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 (Suppl. Fig. S1). This suggests that processing schemes aimed at equalizing noise amplitudes fail to suppress the effect of temporally varying noise sources in an 120 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 numbers and the most consistent with time, over a duration of 10 days. We handle the short-term correlation data, filtering, clustering, and stacking with the python module ruido (see "Code and data availability").

125
We use the stretching method to measure relative changes in velocities (e.g. Sens-Schönfelder and Wegler, 2006). 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 130 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 seconds and 16-40 seconds windows). Figure 3 shows results at station G.UNM, 0.5-1 Hz, using the coda window at 16-40 seconds. 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 135 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.

Modeling velocity changes 140
Visual inspection of the time series in Figure 3 suggests that they contain different patterns: seasonal changes, a co-seismic drop coinciding with the 2017 M w 7.1 Puebla and 2020 M w 7.4 Oaxaca earthquakes, as well as 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., 2014Hobiger et al., , 2016Wang et al., 2017;Donaldson 145 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 Eulenfeld, 2019). 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 mostly insensitive to changes in P-wave velocity (see sect. 5.2).

Earthquake effects
The 2017 M w 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 M w 7.4 Oaxaca (Sta. María Xadani) earthquake, see ranging from τ min to τ max after the time of the earthquake, t quake : Here, s is a negative value indicating the drop of seismic velocity from its previous value c 0 . For intermediate timescales, (τ min ≪ t ≪ τ max ), this model captures the approximately logarithmic time dependence, while it tends to c 0 for large times, and is finite for t = 0 where a purely logarithmic recovery is not defined. Hypothesizing that the minimum timescale τ min is 165 below the temporal resolution of our measurements, we set τ min = 0.1 s 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).

Seasonal effects 170
The following sections describe how we model the effects of surface temperature (4.2.2) and precipitation (4.2.3) on Rayleighwave phase velocity. In both cases, we use analytical solutions to diffusion-like problems for a homogeneous half-space to model their effect on v s from the surface to 2 km depth in terms of a) thermo-elastic 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 thermo-elastic model, annual and sub-annual periodic surface temperature changes 175 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 (Fichtner, 2020). 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 v s -change to obtain the predicted surface-wave phase velocity (c) changes.

180
The solutions at depth in the homogeneous halfspace 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 (Roeloffs, 1996). 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 6 homogeneous half-spaces characterized by a) 3 and b) 2 trial diffusivities in the ranges a) 0.15 to 2 mm 2 /s (Richter et al., 2014) and b) 0.0001 to 4 m 2 /s (Roeloffs, 1996). We retain the 185 diffusivity for each site that produces the best fit across all components and frequency bands.

Site-specific velocity profiles
We need estimates of shear wave velocity v s , compressional wave velocity v p , and density ρ to establish the effect of velocity changes at depth on surface-wave velocity changes dc c . To account for the different subsurface geology at each station, we use a classification as "hard", "intermediate", and "soft" provided by Quintanar et al. (2018) and we consider the aquitard thickness  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 at greater depths.
The values for geologic structure and seismic properties are based on a synopsis of the work of Pérez Cruz (1988); Chavez-Garcia and Bard (1994); Singh et al. (1995Singh et al. ( , 1997 and Shapiro et al. (2001), which serve as references for modeling the seismic 200 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. Supplementary Fig. 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. (1995Singh et al. ( , 1997. Future investigations would benefit from a more detailed 3-D velocity model of the basin and surroundings. As the surface-wave sensitivity kernels in Supplementary Fig. S3 illustrate, the model leads to the following behaviors. At   (1988); Singh et al. (1995Singh et al. ( , 1997; Shapiro et al. (2001). for the higher frequencies. Our single-station analysis presumes a laterally homogeneous structure near each station. While this 210 is a rather crude assumption, we believe that it captures important aspects of the basin, in particular the very shallow sensitivity of our observations at soft sites.

Surface temperature effects
We adopt the approach of Richter et al. (2014) to compute the thermo-elastic 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, (see eq. 14 of Richter et al., 2014), where we summarize material parameters into a depth-independent temperature sensitivity: (eq. 12 of Richter et al., 2014), where b = 1+ν 1−ν for S-waves, ∂ρv 2 ∂σc describes the change of shear modulus with respect to stress, 220 i.e. the non-linear elastic rheology, ν is Poisson's ratio, and α the linear thermal expansion coefficient. The single factors of this sensitivity are not particularly well known, so we fit the product s t 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 Climatico, 2022). We expand the temperature curves into Fourier series with five terms using the nearest available meteorologic station to each seismic station. This proved important for the model fit, as the 225 temperature variations with sub-annual period affect dc c . Finally, we convert the shear-wave velocity change at depth to dc c using the phase velocity sensitivity kernels for the fundamental-mode Rayleigh waves. During the inversions we use 3 trial values of thermal diffusivity κ t in the range of 1.5 · 10 −7 to 2 · 10 −6 m 2 /s and select the κ t at each station that produces the minimum average misfit over all components and frequency bands.

Hydrological effects 230
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 v s to dc c with surface wave sensitivity kernels. 235 We compute the pore pressure change using the impulse response to hydraulic head change derived by Roeloffs (1988) for a homogeneous halfspace, assuming test values for hydraulic diffusivity κ h of 1 m 2 /s representing relatively impermeable unconsolidated sandy sediments and 10 −4 m 2 /s representing permeable unconsolidated clay sediments (Roeloffs, 1996). Note that the sensitivity of the modelled dc c 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 ( Supplementary Fig. S4). We estimate the change in v s due 240 to pore pressure changes using granular media theory (Mavko et al., 2020). Saul and Lumley (2013) formulated the effective moduli for effective pressure p eff , based on which: (see also Takano et al., 2017), where p eff = P ov − nP p with overburden pressure P ov and pore pressure P p . We estimate overburden pressure as lithostatic pressure from the 1-D density models described in sect. 4.2.1, and unperturbed pore pressure 245 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 to infinity for effective pressures approaching zero. We 250 mitigate this by introducing a minimal effective pressure p 0 as a parameter in the inversion, i.e., p eff = max (p 0 , P ov − nP p ).
The resulting dc c is strongly sensitive to p 0 .

Probabilistic inversion
Considering the superposition of the temperature, precipitation, earthquake, and linear trend, we aim to model the phase velocity change dc c as where κ t , κ h are the thermal and hydraulic conductivity, respectively, s t is the sensitivity of the shear wave velocity to thermoelastic stress, p 0 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  This study Surface pp from rain (Roeloffs, 1988) Stress-dependent vs using granular medium theory Rayleigh wave Kvs 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 GPa and 35 GPa for the pore water and rock matrix, respectively. We use these values to estimate Skempton's B following Roeloffs (1988), and adopt their estimate of the undrained Poisson ratio. We conduct an MCMC inversion to determine the parameters s t , p 0 , ∆c, τ max , a, and b using the Euclidean distance between the observed and modeled dc c with the emcee python package Comparison between observed dc c 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 Hz and 2.0 Hz, where a gap in highly coherent observations results in an offset of dc c between 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, e.g., by Hobiger et al. (2014) is preferable. We conclude 275 that strong-motion stations yield overall good results for dc c studies in urban settings with high-amplitude anthropogenic noise comparable to G.UNM.

Observed and modeled velocity changes for RSVM accelerometers and G.UNM seismometer
Results at 2-4 Hz for strong-motion stations of RSVM and the seismometer at G.UNM (Fig. 5)  Our model generally fits better for stations inside the basin, including transition zone and lake zone station, than those outside. This is reflected by a consistently higher CC between observed and modeled dc c for stations inside the basin, and 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 290 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-v s aquitard (Solano-Rojas et al., 2015) and to shallow profiles of v s from well logs (Singh et al., 1995(Singh et al., , 1997. This results in likely more accurate station-specific 295 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 dc c , usually focusing on a small number of selected stations (Tsai, 2011;Richter et al., 2014;Lecocq et al., 2017;Fokker et al., 2021;Illien et al., 2021Illien et al., , 2022 or time series (Rodríguez Tribaldos and Ajo-Franklin, 2021). Previous array-wide analyses have mostly focused on empirical transfer func-300 tions between hydrologic and meteorologic parameters and observed dc c (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 dc c (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 305 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 dc c (see Illien et al., 2022, for a discussion on this point). Nevertheless, it is the first step to a comprehensive model for dc c based on simplified physics at the basin scale.

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 v s to thermal stress, s t (which summarizes non-linear and linear elastic properties), and (iv) the minimum effective pressure p 0 . For practical reasons, only few values were tested for (i) and (ii). Therefore, we only present the inference of s t and p 0 by the inversion. We observe variability between the Z, N, and E components for both parameters, as well as between 315 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 supplementary material.
In previous studies, single-station correlations have been interpreted as multiply scattered body waves (Richter et al., 2014;Sánchez-Pastor et al., 2018, e.g.), whereas we interpret them as surface waves. This is supported by elliptical particle motions ( Supplementary Fig S7). The coda of cross-correlation is usually interpreted in terms of surface waves ( Table 2).The main 320 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 sensible. Apart from the study by Yuan et al. (2021), the single-station correlation sensitivity to velocity 325 changes remains scarcely researched in comparison to cross-correlation coda sensitivity. This topic merits further attention.
Assuming that surface waves dominate our observations, we consider that the variation of parameter values with frequency range is an effect of surface-wave dispersion. As a working hypothesis for differences of parameters between components, we assume that waves recorded on E, N, and Z-component are dominated by different scattered arrivals or different wave modes and therefore are sensitive to different parts of the surrounding medium.

Sensitivity of v s to thermal stress
Inverted values for the depth-independent sensitivity s t of v s to thermal stress range from 10 −6 to 0.1 K -1 (all sites) and 10 −6 to 0.01 K -1 (only basin sites) (Fig. 6). Inside the basin, s t ranges from 10 −6 to 0.01 K -1 . For the frequencies 1-2 and 2-4 Hz, we observe that stations with high sensitivity to surface temperature variations are mostly located in and near the lake zone, 335 whereas lower s t is found mostly in the transition zone. This split disappears at 4-8 Hz.
When considering equation 3 and assuming that Poisson's ratio ν ranges approximately from 0.2 to 0.5, thermal expansion α from 10 −6 to 10 −5 K -1 , and elastic non-linearity ∂ρv 2 ∂σc from 10 to 1000 (Richter et al., 2014), then a physically reasonable range of s t is 3 · 10 −5 to 6 · 10 −2 K -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, s t estimates are beyond 340 what appears 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 sensivitity.
A second possible reason is the fact that we neglected the depth term of Berger's thermo-elastic solution. If we consider that 345 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.

Minimum effective pressure p 0
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).

350
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 p 0 of up to 10 kPa are in a reasonable range.
Our current model of poro-elastic stress is based on a homogeneous halfspace solution (Roeloffs, 1988), whereas the central Mexico City basin is known to have a complex hydrologic structure with an aquitard and several underlying aquifers between 355 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 dc c could be constructed if measured hydraulic head of groundwater wells at high temporal resolution were available. The results of p 0 are nevertheless informative. Low values of p 0 (i.e. a high value on the right hand side of equation 4) suggest that v s is sensitive to precipitation at a site, while a high p 0 indicates the opposite.

360
Lower effective pressures in the hill and transition zones show that v s 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-Guerrero, Figure 6. Inferred 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. Grey 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. Shaded relief from SRTM (Farr et al., 2007), geotechnical zonation from Gobierno de la Ciudad de México (2017). White triangles show excluded models that had a misfit above the threshold. Grey 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. Shaded relief from SRTM (Farr et al., 2007), geotechnical zonation 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 equation 4. 2004). Two stations in the lake zone (MULU, VRVM, see Fig. 1) also show low p 0 at 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 365 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 dv v and HVSR, especially 370 considering that the latter is commonly used for site effect assessment. We present the co-seismic drop and maximum relaxation timescale τ max of the 2017 Puebla earthquake with respect to peak ground acceleration (PGA) in Figure 8. Markers show median models of velocity drop and τ max , error bars show the range between the 16th and 84th percentile. Not all strong-motion stations of the RSVM network recorded the Puebla earthquake;

375
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 (Figure 1). We add a dither ( are not well explained (R 2 = 0.28 for τ max and R 2 = 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 385 than the drop during the 2017 Puebla earthquake.
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. earthquake and aftershocks, and found a τ max of 846 days for the main shock, much longer than for the aftershocks (155 days).
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 days for the Kanto basin (Viens et al., 2018), 250 days in Nepal (Illien et al., 2022), and 3 years in Parkfield (Brenguier et al., 2008a;Wu et al., 2016).

400
A possible explanation for the inferred long τ max is that the relaxation functions used to model dc c do not account for permanent damage. Permanent changes in the interferometric waveforms may lead to poor estimates of the relaxation time scale.
Permanent damage can be assessed using the decorrelation (Larose et al., 2010): in Figure 9. Here, CC i is the correlation coefficient of the stretched traces with respect to the average trace for 2017. Stations that 405 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, 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 410 a proxy for permanent damage, this suggests that, compared to other sediments, the shallow clay of the lake zone would suffer less permanent damage.
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 / v s10 . We use the first 10 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 415 field studies (Viens et al., 2018). Figure 10. Slope of the linear term in the dc c model versus 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.

Residual linear trend
Following the visual appearance of the dc c time series, we introduced a linear term in the dc c 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 dv 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 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 425 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 due to groundwater extraction and sediment compaction (Solano-Rojas et al., 2015;Chaussard et al., 2021;Cigna and Tapete, 2021; 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 435 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 (R 2 = 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 %/year, which they interpreted as an effect of poroelastic contraction. Because 440 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.

Conclusion
We presented a comprehensive study of seismic velocity changes underneath Mexico City. Our study has several innovative 445 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 set-450 tings 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-v s lacustrine sediments and at intermediate sites are mostly related to shear wave velocity changes in the top 100 meters, 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 455 most sites, observed seasonal velocity changes show clear differences between 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 460 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. 465 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 %/year) and may provide additional depth-sensitive in-470 formation about the compaction processes, it certainly merits further investigation also 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 (doi:10.18715/GEOSCOPE.G). Data from the Red Sísmica del Valle de México are 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 Ac-475 celerographic 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 is distributed through the Accelerographic Database System: https://aplicaciones.iingen.unam. mx/AcelerogramasRSM/. Topographic data used for shaded reliefs was obtained via PyGMT (Uieda et al., 2021) and is based on NASA SRTM data (Farr et al., 2007). Copernicus Sentinel-1 SAR data was retrieved from ASF DAAC, processed by ESA.