Valley of Mexico site response from non-reference site approach Open Access

Site response is fundamentally a ground motion difference between two site conditions. One is a reference condition, generally involving relatively competent materials, while the other is the site condition at a particular site of interest. Reference site and non-reference site approaches can be used to measure these ground motion differences (Field and Jacob, 1995). Site response in the Valley of Mexico, which contains Mexico City, has typically been evaluated in past work using reference site approaches (Borcherdt, 1970, 1994), in which ground motions on various site conditions are normalized by ground motions from reference sites. The site amplification that is derived depends on the attributes of the reference site or sites, which have usually been taken (for applications outside of the Valley of Mexico) as bedrock sites. Non-reference site approaches define the reference condition not for a single site, but for a parameterized site condition as derived from seismological inversions or as used in ground motion models (GMMs); large values of time-averaged 30 m shear wave velocity (V_S30 = 760–1000 m/s) are typically selected.

Current practice is that ground motion estimations for sites in the Valley of Mexico are based on a reference site approach having two basic steps (e.g. Ordaz et al., 2024; Ordaz et al., 1994; Reinoso and Ordaz, 1999). The first step pertains to the estimation of reference site ground motions, which are for a specific location on the UNAM campus in the hill zone adjacent to the Valley of Mexico lakebed (referred to as the CU site). The second step estimates site response for various locations in the Valley of Mexico relative to the reference site. There are several drawbacks to this traditionally applied approach, including the need for a GMM for a single reference site, which necessarily produces data sampling problems (the available data for a single site are limited) and divorces the ground motion estimation for the Valley of Mexico from locations elsewhere in Mexico, which use GMMs derived from much larger databases.

We revisit the well-studied problem of Valley of Mexico site response using a non-reference site approach, which to our knowledge has not been applied previously for this important locale where site response has produced structural damage (e.g. Mayoral et al., 2019a; Singh et al., 2018). The benefit of this alternative approach is that it can produce site factors that are consistent with regional GMMs, thus avoiding the need for distinct ground motion estimation procedures in Mexico City and the remainder of Mexico. Such approaches have been applied elsewhere to derive regional (e.g. Parker and Stewart, 2022) and local (Wang et al., 2022) site response models in subduction zone regions. We first review the lakebed geology and zonation procedures used to characterize spatial variations in site conditions. We then synthesize site data for the purpose of developing site parameters for ground motion stations and for estimating values of site parameters in geotechnical zones. We describe current procedures used for ground motion estimation in the Valley of Mexico and provide example results. We then describe the application of a non-reference site approach with available ground motion data, from which a subregional site response model (i.e. model specific to the Valley of Mexico) is developed that is conditional on V_S30. The model is compared to site response estimated using currently applied reference site procedures—large differences are observed that are related to appreciable site response at the reference site—and is validated against a testing data set from the 1985 M8.0 Michoacán earthquake independent of the training data. We conclude with recommendations on how the present results can be utilized to improve ground motion hazard estimation in Mexico City. A previous version of the material presented in this article appeared in the dissertation of the first author (Contreras, 2022) and a subsequent report based on that dissertation (Contreras et al., 2023a). Portions of the work presented here and in a companion paper (Contreras et al., 2025) have also been the subject of conference presentations (Contreras et al., 2023b; Contreras et al., 2024).

Mexico City is located within the eastern sector of the Trans-Mexican Volcanic Belt, north of the volcanic front (i.e. on the backarc side). Figure 1 depicts the geology of the subregion, as adapted from Flores-Estrella et al. (2007) and Arce et al. (2019). Zones delineating different geotechnical conditions are shown in Figure 2. These zones were originally derived from geologic maps and geotechnical logs (Marsal and Mazari, 1969; Romo et al., 1988). As indicated by Lermo et al. (2020), the zones were refined by Normas Técnicas Complementarias para el diseño por Sismos (NTCS) (2004) in consideration of site period maps by Lermo et al. (1988) and Lermo and Chávez-García (1993, 1994a, 1994b) and have been used for site classification purposes.

Within the Valley of Mexico, the lakebed contains lacustrine sediments that are remnants of ancient Lake Texcoco. The lake occupied the entire Valley of Mexico between the Pleistocene epoch and the last glacial period (about 11,000 years ago). The lake zone is indicated in Figure 1 by the Ql unit and includes the Texcoco Lake (while this was once encompassing the entire Valley of Mexico, as used here it refers to a smaller lake in the post-Aztec era), Xochimilco Lake, and Chalco Lake (Romo et al., 1988). These smaller lakes existed as recently as 1245 (upon arrival of the Aztecs), after which the lakes were impacted by human activity (Alcocer and Williams, 1996). Geotechnical conditions in these lake zones consist of a surface layer of desiccated alluvial deposits overlying thick and soft lacustrine clay interbedded with thin seams of sands, silty sands, volcanic glass, and fossils, which is locally known as the Upper Clay formation. Underlying this formation are hard deposits of silty, weakly cemented sands, a deeper clay stratum known as the Lower Clay formation, and compact lacustrine cemented silty sands and gravels. The Xochimilco–Chalco Lake consists of a clay deposit (somewhat stiffer than the Upper Clay formation in Texcoco) with interbedded seams of silty sands, silts, and sands. This clayey deposit is underlain by a basalt layer (lava flow). The lakebed area has four subzones (IIIa–IIId in Figure 2) that are differentiated in part by the period of the soil column (1.0–1.5 s, 1.5–2.5 s, 2.5–3.5 s, >3.5 s, respectively) (NTCS, 2004).

As shown in Figure 1, Quaternary alluvium occurs between the lake deposits and hill areas to the west and south, producing a geotechnical transition zone (denoted Zone II by NTCS, 2004; shown in yellow in Figure 2). These alluvial deposits have variable sequences of firm soils, sands, silty sands, and soft clays. The hill areas west of Lake Texcoco consist of Tertiary age epiclastic rock and alluvial fans (Flores-Estrella et al., 2007). Geotechnical materials encountered near the surface generally consist of silty sands with gravels and cemented tuffs. To the west and south of Xochimilco–Chalco Lake, Quaternary basaltic, andesitic, and epiclastic deposits up to 20 m in thickness overlie these formations (Flores-Estrella et al., 2007; Arce et al., 2019). These volcanic deposits were derived from a series of volcanoes in the region, including Iztaccíhuatl, Popocatépetl, La Malinche, Ajusco, and Nevado de Toluca. As shown in Figure 2, the hill zone is denoted as Zone I by NTCS (2004).

We have compiled available information on shear wave velocity (V_S) profiles, site periods (T_S), and depths to a firm horizon (Lower Coarse Grain layer, $zLGC$⁠) in the Valley of Mexico. The V_S profiles and T_S data are from the following technical reports, research papers, and student dissertations: Cárdenas-Soto et al. (2016), Chamorro (2016), De La Rosa et al. (2013), Díaz-Rodríguez (2003), Flores-Guzmán et al. (2014), Luna (2012), Mayoral et al. (2011), Mayoral et al., (2016), Mayoral et al. (2017), Mayoral et al. (2019b), Pérez (2017), Ramos-Martínez et al. (1997), Romo and Garcia (2003), Seed et al. (1987), Stephenson and Lomnitz (2005), Vergara-Huerta and Aguirre (2020), Wood et al., (2019, 2024, 2023). The $zLGC$ information is from a contour map provided by Juárez-Camarena et al. (2016).

Velocity data come from invasive (downhole, crosshole, suspension logging, seismic cone penetration) and noninvasive methods (Multi-channel Analysis of Surface Waves—MASW, Microtremor Array Measurements—MAM, Modified microtremor seismic method using SPatial AutoCorrelation—MSPAC, seismic interferometry—SI, seismic refraction). Site period data are derived from microtremor-based measurements of Horizontal-to-Vertical Spectral Ratios (mHVSR), the lowest frequency peaks of which reveal the fundamental mode site period (Nogoshi and Igarashi, 1970, 1971; Field and Jacob, 1993, 1995; Bonilla et al., 1997, 2002; Cadet et al., 2012; Satoh et al., 2001; Theodulidis et al., 1996). With the exception of 23 sites from Wood et al. (2019, 2024, 2023), site periods are provided without the source mHVSR curves, which means that the procedures used to select site periods may be inconsistent and other parameters (such as the height of peaks) are not available. For sites with V_S profiles that reach a large impedance contrast (i.e. below the lower clay formation in Texcoco Lake or the basalt layer in Xochimilco–Chalco Lake), site period was also derived from V_S profiles using the Simplified Rayleigh Method as given in Urzúa et al. (2017) using soil layers generally down to $zLGC$⁠.

Table A1 in the electronic supplement materials presents the compiled information for sites with a V_S profile, including the measurement locations, zones (from Figure 2), V_S30 values, derived site periods co-located with the V_S profile where available, testing methods, $zLGC$ values, and references. Two columns for site period are provided, one for results of mHVSR testing and one for period derived from the V_S profile. Figure 1 shows the locations of these sites.

Figure 3 shows box and whisker plots of V_S30 and site period by zone. Both sources of site period (measured and derived from the V_S profile) are used; for sites with both types of site period, the mHVSR result was used. Where an mHVSR measurement is not co-located with the ground motion site, we used values of site period from the large compilation of Lermo and Chávez-García (1994a) (409 sites within the Valley of Mexico). Median and mean values of V_S30 and site period are indicated by the blue and red lines in Figure 3, respectively. The lower and upper quartiles are indicated by the boxes and the extremes are indicated by the whiskers (extreme horizontal lines). Results in Figure 3 show that, on average, the three main zones (i.e. Zone I-Hills, Zone II-Transition, and Zone III-Lake) have statistically distinct values of V_S30, as judged by non-overlapping boxes. The subdivision of Zone III shows significant differences between V_S30 values in Zone IIIa and in the other zones (IIIb, c, d); however, the differences between these latter zones are lower and not necessarily significant. In the case of site periods, the zones have statistically distinct site periods. An exception is Zones IIIc–d, which present a similar median. Zone IIIc shows a large variability in comparison with the other zones, which may reflect natural variability or may be a consequence of data quality; further research would be needed to assess and potentially improve site period characterization in this area.

Figure 4a shows the relationship between V_S30 and site period. The negative correlation is strong, with short-period sites having high V_S30, and V_S30 decreasing with increasing $TS$⁠. This result was expected given that about 57% of the data shown in Figure 4a are from measured V_S profiles that were constrained using the measured mHVSR site periods. Figure 4b shows a similarly strong negative correlation between V_S30 and $zLGC$⁠, which is considered below in model development.

We apply protocols for characterizing V_S30 for ground motion sites as described by Ahdi et al. (2022). That study classified different methods of V_S30 prediction for sites without seismic velocity measurements. The data compilation in Table A1 was not available at the time the NGA-Sub database was assembled. Hence, while the general approach is the same as given by Ahdi et al. (2022), the assigned site parameters are different for sites that were already in the database, and many new sites have been added to the database.

The database contains 97 ground motion sites in the Valley of Mexico. Of these, 17 are Code 0 sites per Ahdi et al. (2022; their Table 5), meaning that V_S30 was computed from a V_S profile, which is taken as the exponent of the mean of the log values (⁠$μlnV$⁠) and the natural log standard deviation (⁠$σlnV$⁠) is taken as 0.1. To assign a V_S profile to a ground motion site, we require separation distances ≤300 m and both locations to be in the same zone.

For the 80 sites without a V_S profile, we assigned $μlnV$ based on the zone the site is located in. This is an example of a region-specific prediction model, which is indicated as Code 2MC (Code 2 is as given in Table 5 of Ahdi et al. (2022); the letters “MC” are added for this Mexico City–specific model, which was not available when the NGA-Sub database was developed). The natural log means and standard deviation for each zone are given in Table 1, which are taken from the results in Figure 3.

Current seismic design in Mexico City uses response spectra that are derived using a two-step procedure:

1. Fourier amplitude spectra (FAS) are derived for a reference site in the hill zone (Ciudad Universitaria, CU site).

2. Reference site FAS are then modified using transfer functions for different locations and the resulting site FAS are converted to response spectra using random vibration theory.

These procedures are implemented for applications in an online software package, SASID (System of Seismic Actions for Design; https://sasid.unam.mx, last accessed October 2022, registration is required).

Four GMMs have been developed for the CU reference site (Arroyo et al., 2024; Jaimes et al., 2015; Jaimes et al., 2006; Ordaz et al., 1994). The number of recordings that are considered are equal to the number events and range from 21 to 23 for the first three models and 67 for the fourth (40 for “coastal” events, mainly interface, and 27 for intermediate depth events). These GMMs are referred to here as “single-site GMMs” to differentiate them from conventional GMMs that are derived considering data from many sites. It is currently unclear whether the GMM that is applied in the SASID software is a modified version of Ordaz et al. (1994) or the newer Arroyo et al. (2024) model. Contreras et al. (2023a) provided further details on the single-site GMMs and how they are applied to derive reference site spectra.

The transfer functions applied to modify the reference site FAS use data from recording sites in the Valley of Mexico and the CU reference site. The precise manner in which these calculations are performed in the current version of SASID is not well documented. The available documentation is provided by Reinoso and Ordaz (1999) and Ordaz et al. (2024), who describe procedures, whereby transfer functions are developed for a series of sites within each of the zones defined in Figure 2 (hill, transition, lakebed). As described by Reinoso and Ordaz (1999) and Ordaz et al. (2024), transfer functions are computed for each instrument location, and for locations between instruments, an interpolation procedure is applied.

Issues with these procedures that are addressed in the remainder of this article include:

1. The amount of data used to constrain the single-site GMM is small, which does not allow for the evaluation of important effects such as magnitude–scaling relationships that include region-specific break points, near-source saturation effects, region-specific anelastic attenuation effects that distinguish forearc and backarc attenuation, and differences between event types (interface vs intraslab).

2. The single-site GMM produces different ground motion estimates for the CU site than would be provided by regional models (e.g. Contreras et al., 2025; García-Soto and Jaimes, 2017; Jaimes and García-Soto, 2020; Kuehn et al., 2023), producing an artificial discontinuity in seismic hazard around the margins of Mexico City.

3. The site transfer function does not use site information (other than location) from the site of interest. As a result, the site information is assumed to be adequately represented by the averaging of transfer functions from neighboring sites. This does not allow for a local deviation of site condition, for example, from a filled stream channel or similar.

4. Nonlinearity in site response is not considered.

5. Due to limited documentation, the procedure is not reproducible. Analysts can apply the procedure through the SASID software but cannot independently apply or check it.

We illustrate site responses derived using the reference site approach as provided by the SASID software for 65 sites in the Valley of Mexico at the locations shown in Figure 2 that recorded three large events in 2017–2018. The SASID program requires only location as input. Effects of site period are implicit to the interpolation procedure. Nonlinearity is not considered, so parameters related to the strength of shaking are not required. As shown in Figure 5, response spectra amplification for sites in each of the five zones relative to Zone I follow a consistent pattern, whereby amplification decreases slightly from 0.1 s to about 0.2–0.3 s, increases to a peak that is variable in amplitude and peak period by zone, and then decreases at longer periods. The peak amplification values range, on average, from 1.0 in ln units (∼2.7) for Zone II to approximately 1.9–2.5 in ln units for lakebed zones (∼7–12). The period of the peaks increases with zone number, being about 0.8 s for Zone II and ranging from 1.5 to 3.0 s for Zone III. Within-zone dispersion across sites is modest.

In this section, we first apply a non-reference site approach to estimate site response at ground motion recording sites in the Valley of Mexico. The estimated site responses are then used to develop a sub-regional site response model.

The ground motion data used in this study are free-field acceleration recordings that include both body and surface waves, each of which contributes to the spectral accelerations used to infer site response. The analysis of ground motion data to evaluate site response effects follows a non-reference site approach (adapted from Field and Jacob (1995)). In this approach, residuals are computed between recorded data and a GMM, which are then partitioned to evaluate site effects. By using a GMM to develop the comparison spectra, our application of the non-reference site approach differs from Field and Jacob (1995), who used an inverted Fourier amplitude spectrum, although at a conceptual level, the approaches align. The GMM-based application of the non-reference site approach used here has extensive precedent (e.g. Parker and Stewart, 2022; Seyhan and Stewart, 2014). The data considered in this investigation are an expanded version of the NGA-Sub database—the expansion is to include additional events that had not originally been considered, including three large events in 2017–2018. Information on the expanded database is provided in Contreras et al. (2023a) and in a companion paper (Contreras et al., 2025) and is not repeated here. The GMM that is applied is the Parker et al. (2022) model developed in the NGA-Sub project with adjustments to account for additional anelastic attenuation in the backarc region (Contreras et al., 2023a, 2025).

The total (unpartitioned) residuals are computed as the difference between data and model,

where (⁠$Yij$⁠) is a ground motion intensity measure for event $i$ recorded at site $j$⁠, $μln,ijr$ is the mean ground motion prediction for reference rock site conditions in natural log units from the adjusted GMM (with the arguments of moment magnitude, $M$⁠, event-type parameter, $FS$⁠, and rupture distance $Rrup$⁠), and $FV,j$ is an ergodic site amplification model (Parker and Stewart, 2022) with linear and nonlinear components,

The superscript $v$ in Equation 1 indicates that the residuals were computed using a site amplification model. Total residuals are partitioned using mixed effects analyses to identify systematic event terms (⁠$ηE,i$⁠), site terms (⁠$ηS,jv$⁠), and remaining, essentially random, residual components (⁠$εij$⁠),

The mixed-effects analyses were performed in MATLAB using the fitlme command (documentation available at https://www.mathworks.com/help/stats/linear-mixed-effects-models.html).

For the investigation of site response nonlinearity in the next section, it is useful to consider within-event residuals relative to a reference rock condition, $δWijr$⁠. These can be computed as,

where $Flin$ and $Fnl$ are the two components of the ergodic site response model (⁠$FV$⁠); note that while both components depend on the site, the nonlinear component also depends on the event because different events produced different shaking intensities and hence have different levels of nonlinearity.

The addition of the ergodic site amplification model in the $δWijr$ calculation (equation 4) offsets the subtraction of this same model in the original residuals calculation (equation 1). This is desirable because (1) the global model may be anticipated to be ineffective in the Valley of Mexico because of its unusual geologic conditions and because applications require extrapolation of the model to V_S30 values below the operable range (<150 m/s) and (2) it provides residuals referenced to a 760 m/s site condition, which is needed to quantify site response effects as used in model development.

Figure 6 shows example plots of $δWr$ (representing site amplification) against the strength of shaking on reference rock (denoted $PGAr$⁠) for two intensity measures (PGA and PSA(2.0 s)) and sites within Zones I and IIIa. The results in Figure 6 are from three well-recorded earthquakes (EQ1: 2017 M8.3 Tehuantepec, EQ2: 2017 M7.2 Puebla-Morelos, EQ3: 2018 M7.2 Pinotepa Nacional). Two of these events (Tehuantepec and Pinotepa Nacional) produced weak shaking conditions (⁠$PGAr≈0.0015−0.0025g$⁠) and one (Puebla-Morelos) produced relatively strong shaking (⁠$PGAr≈0.06−0.09g$⁠). The results suggest appreciably lower levels of PGA site response for the Puebla-Morelos event than for the weaker shaking events, with the level of reduction being greater for the softer site condition (Zone IIIa) than for the firmer condition (Zone I). For PSA(2.0 s), the site responses for the three events are largely similar. These features are consistent with expectation when site response is nonlinear, which is known to reduce short-period site response, has little effect for long-period site response, and is of increasing importance as site conditions soften. These site response features, while consistent with many prior observations worldwide (e.g. Parker and Stewart, 2022; Seyhan and Stewart, 2014), are not incorporated into current design procedures (SASID) for the Valley of Mexico, which assume site response to be linear. This study is not the first to identify nonlinear site response in Mexico City; Singh et al. (1988b) and Sánchez-Sesma et al. (1988) identified nonlinear site response during the Michoacán earthquake at specific sites (Central de Abastos Oficina, CDAO and Central de Abastos Frigorífico, CDAF), although other investigators have disagreed (Chávez-García and Bard, 1994) and the topic has remained controversial (Flores-Estrella et al., 2007).

The evidence for nonlinearity shown in Figure 6 was derived from response spectral ordinates, which have some dependence on spectral shape. As shown by Stafford et al. (2017), spectral shape effects are important for small magnitudes (generally < 4) for which linear site response in Fourier amplitudes can appear to be nonlinear in PSA due to sensitivities of earthquake corner frequencies to magnitude within the frequency range of oscillator response saturation. However, given the M range of the data considered (7.2–8.3), such effects are very unlikely to explain the nonlinearity identified in this study.

For Zone I (hills), most of the recording stations (15) are located to the southwest of Mexico City (Figures 1 and 2). However, there are five stations in other areas of the city mapped as hills, three of which presented significant differences in site response. One of these stations is located on the slopes of an isolated hill within the lake zone (Cerro del Peñón) and two are located on the slopes of hills located in the north of the city (Cerro Tepeyac and Cerro de Guerrero). Site responses for these three stations are shown in Figure 6 with hollow triangles. The Cerro del Peñón station shows higher within-event residuals than the average for Zone I (possibly influenced by the lakebed), whereas Cerro Tepeyac and Cerro de Guerrero present lower within-event residuals than average for Zone I. Due to this anomalous behavior, the data from these three stations were not used to evaluate the overall nonlinear site response of Zone I.

To model the nonlinear components of site response, we fit the data in Figure 6 using Equation 5,

where f_1,z to f_3,z are coefficients for zone $z$⁠, and $PGAr$ is as defined previously. Parameter $f1z$ represents the linear component of the site response for each zone z, whereas the second term in the sum is $Fnl$⁠. These fits were performed setting $f3z=0.01g$ and then regressing for $f1z$ and $f2z$ using least-squared procedures. Almost any reasonable value of $f3z$ could have been selected for the present application. The reason we selected 0.01 g was because it produced values of $f2z$ that are generally compatible with prior experience from subduction zone site responses (for similar V_S30 values and IMs). This is shown in Figure 7, where we plot the regressed values of $f2z$ as a function of period for the six zones. Predictions from Parker and Stewart (2022) and smoothed values of the Valley of Mexico coefficients are also shown. The Parker and Stewart (2022) model uses a value of $f3z=0.05g$⁠, which is constrained from global data generally from firmer site conditions than in the Valley of Mexico. Parameter $f3z$ represents the level of reference site shaking level where nonlinear effects begin to become appreciable. Given the very soft soil conditions in the Valley of Mexico, larger strains could be expected for modest shaking conditions than for relatively firm sites; as a result, a reduction of $f3z$ to account for this effect has a physical basis.

The results in Figure 7 show strongly negative values of $f2z$ at short periods for all zones. For softer sites (e.g. Zones III), the values are more negative than for relatively firm sites (Zone I) and the negative values extend to longer periods (e.g. at 1.0 s, $f2z≈0$ for Zone I but remains negative for all Zone III groups). A smooth fit to regressed coefficients was provided by eye as shown in the figure. The resulting period-dependent values of $f2z$ are given in Table 2.

To model the linear component of site response, it is necessary to correct within-event residuals for nonlinearity. This adjustment is made using Equation 6,

where $(Fnl)ijz$ is dependent on event $i$ and site $j$⁠, because of their effect on $PGAr$⁠, and on zone $z$⁠, which determines the coefficients that are applied. As shown by Contreras et al. (2023a), when Equation 6 is applied, the adjusted within-event residuals are effectively independent of $PGAr$⁠. These within-event residuals are then used in mixed-effects analyses to compute site terms referenced to a site condition of V_S30 = 760 m/s, which are denoted $ηSr$⁠,

Equation 7 does not include a constant term, which forces all of the site response bias into the $ηSr$ term and $εij$ is the remaining residual. The physical meaning of $ηSr$ is the mean difference between linearized ground motions for a given site and what would be expected for the site if its site condition was V_S30 = 760 m/s. Subsequently, $ηSr$ is referred to as the “linear site response.”

Typical characteristics of Valley of Mexico site responses are illustrated with example results for seven sites from different zones (site locations in Figure 2). The linear site responses for these seven sites are shown in Figure 8. The Zone I site (CUP5; this is the modern station code associated with the traditional CU site) has appreciable amplification, ranging from 0.5 to 1.7 (natural log units), with the maximum amplification being essentially constant for T = 1.5–5s. The Zone II site (AO24) has higher amplification for PGA and PSA, with a peaked response at 0.8 s that is much stronger than for CUP5 (∼2.5), which then decays sharply for T > 1 s. The Zone III sites (IIIa, IB22; IIIb, SCT2; IIIc, XP06; IIId, AE02) also have larger amplification than CUP5 for PGA and PSA. At the longer periods, the site responses are peaked at amplitudes ranging from 3.2 to 3.9 (factors of 25–50 approximately) at periods of 1.5 s for IIIa, 1.5–2 s for IIIb, 2.5 s for IIIc, and 4–5 s for IIId.

Results in Figure 8 for the CUP5 site are particularly notable, because that site has been widely used as a reference site in prior studie of Valley of Mexico site response (e.g. Reinoso and Ordaz, 1999). The V_S30 for this site is 295 m/s (Table A1, see ID#7) and its site response is appreciable. Strong site response at the CUP5 site was previously recognized by Singh et al. (1995) and attributed to complex shallow structure. Reference sites used in other regions have generally had much stiffer geologic conditions (e.g. 750–1000 m/s; Borcherdt, 1994, 2002).

Site responses relative to a 760 m/s condition are compared to those derived using SASID in Figure 9. The patterns of site responses between zones and across periods are broadly consistent, although the present site responses are larger by amounts ranging from 0.7 to 1.4 in natural log units (factors of 2.0–4.0) due to the different reference site conditions. These substantial differences are largely a consequence of the significant site response at the CU site. This is illustrated by the “adjusted median” lines in Figure 9, which are computed as the sum of the SASID medians and the CU site response. These adjusted medians approximately convert the SASID site response from the CU reference condition to a 760 m/s reference condition. The results are much closer to the median non-reference site response, but are generally somewhat higher (∼0.1–0.5, or factors of ∼1.1–1.6) for T > 1.0 s and lower for T < 1.0 s.

Due to the strong correlation of V_S30 with depth and site period T_S (Figure 4), our model development approach begins by modeling amplification trends based solely on V_S30 and subsequently checking dependencies on $zLCG$ and T_S using residuals analysis. Figure 10 shows site response for all Valley of Mexico sites as a function of V_S30 for the intensity measures of PGA, PGV, PSA (0.3 s), and PSA (3.0 s). At short periods (e.g. PGA), the site response increases as V_S30 decreases from about 400 to 150 m/s. For softer sites(< 150 m/s), the site response does not scale with V_S30. This general pattern is retained up to a period of 2.0 s, but for longer periods, the saturation at low V_S30 gradually disappears, instead continuing to increase as V_S30 decreases to its minimum value of 50 m/s. The data trends in Figure 10 are distinct from the Parker and Stewart (2022) ergodic model when that model is extrapolated below 150 m/s. This is notable because that model is generally compatible with site response trends elsewhere in Central America and Mexico (CAM; Contreras et al., 2025).

To capture the subregional response, the following model is applied:

where the VM superscript indicates it is a Valley of Mexico model, the $vm$ subscripts indicate that the coefficients are specific to that subregion (otherwise the same variable names are used as in the Parker and Stewart (2022) global model), V_ref = 760 m/s, $V1vm$ is a break velocity, $s1vm$ is the V_S30-scaling for $VS30<V1vm$⁠, and $s2vm$ is the V_S30-scaling for $VS30>V1vm$⁠. This model is similar to that used by Parker and Stewart (2022), with the following changes: (1) site response is only modeled for V_S30 < V_ref (higher velocities are not applicable in the Valley of Mexico) and (2) an amplification shift parameter f_VM is introduced that allows for higher site response across all V_S30, which is needed across the full period range. Figure 10 shows the fit of the model in Equation 8 to the data, with model coefficients indicated. The selected function fits the data well.

Figure 11 shows the period-variations of model coefficients. Negative values of slope parameters indicate that amplification increases as V_S30 decreases. For short periods (0.01–0.2 s), the slope parameters have relatively little change with period, with $s1vm~−0.2$ and $s2vm$ ranging from −1.25 to −0.9. The strongest negative slopes are at long periods (>2.0 s for $s1vm$⁠, >0.7 s for $s2vm$⁠), which are needed to model the large long-period amplification. Positive values of $s1vm$ occur for periods of 0.5–2.0 s, indicating a decrease of amplification as V_S30 decreases in this range. Break velocity $V1vm$ ranges from 120 m/s at short periods to 70 m/s for T > 2 s. Amplification shift parameter f_VM ranges from −0.5 to 1, presenting negative values at short periods (0.01–0.2 s) and large, positive values at long periods (>2.0 s). The adopted period-dependent values of these parameters are given in Table 3.

In this subsection, we evaluate the performance of the V_S30-scaling model proposed for the linear site response of the Valley of Mexico, extend the model to consider effects of sediment depth based on trends observed in residuals analyses, check the sensitivity of residuals to site period T_S, and develop models for within-event components of aleatory variability.

To evaluate model performance, we take the linearized within-event rock residuals from Equation 6 as “data” and subtract the proposed model prediction (Equation 8) to evaluate within-event residual $δWv,lin$⁠, where the “$v$” superscript indicates that V_S30-dependent site response is considered:

Figure 12 shows that the model removes the principle trends of within-event residuals with V_S30 for PGA, PGV and PSA at 0.3, 1.0, 2.0, and 3.0 s. This is evident from the binned means of these residuals generally being nearly zero and not exhibiting a trend with V_S30. Some non-zero binned means (⁠$<±0.5$⁠) occur at periods between 1.5 and 3.0 s (e.g. the sites in Zone IIIc with V_S30∼57 m/s have a positive binned mean for PSA at 3.0 s). This is caused by these sites having a natural site period of T_S∼2.5–3.5 s that approximately matches the oscillator period. Figure 12 also shows that Zone I (V_S30∼300 m/s) sites have larger dispersion than other zones (although the binned means are centered at zero). This indicates a relatively large variability of site response in the Hills Zone, likely reflecting heterogeneity of geotechnical and topographic conditions.

Sediment depth has been found for multiple regions to have predictive capacity for site response beyond that provided by V_S30 (e.g. Day et al., 2008; Nweke et al., 2022). Figure 4b shows that depth to the Lower Coarse Grain layer is negatively correlated with V_S30, which makes it poorly suited as an additional predictive parameter. To avoid this problem, we define a model for mean depth conditioned on V_S30 ( $z¯LCG$⁠) as shown in Figure 4b and use depth differential (⁠$δz=zLCG−z¯LCG$⁠) as a parameter against which to evaluate model performance. Figure 13a and b shows trends of residuals (⁠$δWijv,lin$⁠) with $δz$⁠, which are significant for periods T > 1s for $δz<−5$m and $δz>5$ m (means of residuals are zero for $δz=−5to5m$⁠). We fit the data using Equation 10,

Figure 13c shows the trends with period of $d1$ and $d2$⁠, which we smooth as shown in the figure.

Although V_S30 is used as the primary site parameter in the proposed model (equation 8), the site parameter traditionally associated with the different Valley of Mexico zones is site period, $TS$ (Lermo et al., 2020). Accordingly, we investigate whether the proposed model, including the depth correction in Equation 10, is able to capture the known period dependency of Valley of Mexico site response. This is evaluated by computing linearized within-event residuals from Equation 11

Residuals from Equation 11 are plotted against $TS$ in Figure 14 (plots without adjustment for depth effects are provided by Contreras et al. (2023a)). The results for PGA, PGV, and PSA for T≤ 1.5 s show no trend with $TS$⁠. This likely occurs because the intensity measures shown either do not correspond with a particular oscillator period (PGV) or their period is lower than those for Valley of Mexico sites, even in the case of Zone I.

At longer periods, the residuals patterns in Figure 14 illustrate site resonance effects that are not captured by the V_S30-scaling model. Before discussing the results, it is useful to recall the median values of site period in the different zones (e.g. ∼1 s for Zone II and IIIa, ∼2 s for IIIb, and ∼3 s for IIIc). Residuals plots for a given oscillator period often have peaks when that period aligns with the zone median site period; for example, a peak occurs at $TS$ = 2 s in the plots for PSA (2 s) for the Zone IIIb sites, and a peak occurs at $TS$ = 3 s in the plots for PSA (3 s) for the Zone IIIc sites. These peaks indicate that there are resonance features in the site response that are not captured by the V_S30-scaling model. We considered supplementing the V_S30-scaling model in Equation 8 with a pulse model that would introduce peaks at a site period (e.g. as done by Buckreis et al. (2024) for a different region). While such a feature would improve data fit, the level of improvement would be marginal. Consider, for example, Zone IIIb—the overall level of amplification at the 2 s peak is approximately 3.5 in natural log units (Figure 9), whereas the height of the binned means in Figure 14 is approximately 0.25–0.3. Accordingly, the overwhelming majority of the site response is captured by the V_S30-scaling model and a resonance model would introduce a relatively minor modification (as expected given the parameter correlation, Figure 4). Moreover, the development of a model conditioned on site period, which is derived from mHVSR, is impractical at the present time due to the limited availability of mHVSR curves (23 of 89 sites), as indicated previously. For these reasons, we decided against adding this component.

Residuals analyses using the Valley of Mexico data along with the recommended model (i.e. reference site regional GMM from Contreras et al. (2025), which is modified from Parker et al. (2022); linear site response from Equation 8; nonlinear site response from Equation 5) provides within-event residuals that can be partitioned using mixed-effects analysis into site terms (⁠$ηS,jv$⁠) and remaining residuals (⁠$εij$⁠). This final set of residuals includes updated estimates of event terms that account for the effects of the Valley of Mexico site response model. Mixed effects analyses also provide corresponding standard deviation terms, which are the site-to-site standard deviation (⁠$ϕS2S$⁠) and single-station within-event standard deviation (⁠$ϕSS$⁠). These are important components of the overall aleatory variability model. Their values could be expected to differ from values from global models as a result of the specific location and geology of the Valley of Mexico.

Figure 15 shows subregional $ϕS2S$ and $ϕSS$ terms derived for Valley of Mexico sites for the sampled oscillator periods along with a smoothed model. Also shown for comparison purposes are the corresponding standard deviation terms from the Parker et al. (2022) global model exercised for conditions generally appropriate for the Valley of Mexico data (V_S30 = 150 m/s, R_rup = 100 km). The subregional dispersion terms are lower than those from the global model, with the differences being particularly large for $ϕS2S$⁠. These lower dispersions are expected because the variations in path and site response in the Valley of Mexico subregion are clearly less than those present globally.

In this subsection, we validate the site response model for the Valley of Mexico presented in the previous section. The validation uses a data set from the 1985 M8.0 Michoacán earthquake. Ground motion data from 26 stations that recorded this event are included in the NGA-Sub database, 10 of which are located in Mexico City, which is known to have experienced substantial site response (Ordaz et al., 1988; Seed et al., 1987; Singh et al., 1988a, 1988b). However, the 10 Mexico City recordings were not considered in the development of the site response model because of data screening criteria that were applied (matching those from Parker et al. (2022)). The criteria that led to the Mexico City sites being screened mainly relate to maximum usable distance (i.e. the maximum distance for which recordings can be considered without introducing sampling bias, which is denoted $Rmax$⁠) and lack of necessary site metadata (V_S30). The $Rmax$ value for this event was 137 km and the Mexico City sites were at distances of 233–246 km. The V_S30 screen no longer applies because in this study V_S30 values have been provided for all stations using the procedures given in the Site Parameters section.

While the elimination of the Mexico City data was necessary in order to maintain consistent data selection criteria, the strong subregional site response likely causes the amplitudes to exceed noise levels to a sufficient degree that the recordings are representative of the ground shaking that was experienced. Instead of making a “one-time” exception to the screening criteria, we chose to use the data as an independent testing data set for model validation. Because data from the Michoacán event at closer distances (not in the Valley of Mexico) were used by Contreras et al. (2023a, 2025) in the development of the Parker et al. (2022) GMM adjustments, event terms (⁠$ηE,i$⁠) for the modified model are available (Figure 16). The events terms are nearly 0 up to 1.5 s and then increase, indicating that the adjusted GMM is unbiased for T < 1.5 s.

The validation is performed by summing (in natural log units) the adjusted reference rock GMM, the event terms from Figure 16, and the site amplification models (V_S30-scaling only and combined model that includes differential depth effects). $PGAr$ values for Valley of Mexico sites range from 0.013 to 0.015 g (using the event-adjusted Parker et al. (2022) model), which is large enough for modest effects of nonlinearity to appear. The results of these calculations are compared to data in Figure 17 for six sites (of the 10 available) from all of the zones except IIIa, for which no data are available. To examine the impact of the nonlinear model, predicted spectra are shown in Figure 17 both for the linear models and the linear models adjusted by the nonlinear model (Equation 5). The results shown for the V_S30-scaling model applied with the nonlinear model are bracketed by bands of $±ϕSS$⁠. The depth adjustments only affect the SCT and Tlahuac Bombas sites because the others have $|δz|<5m$⁠. Overall, the proposed site response model performs well. The misfit of the model to the data is reduced in 60% of cases by the use of the nonlinear model. This exercise provides confidence that the proposed model is consistent with observations from the landmark 1985 event. SASID results are not provided in Figure 17 because it is not configured to provide results for scenario events (it provides hazard results across multiple events).

This research examined site response within the Valley of Mexico, which contains Mexico City, using a non-reference site approach. The work presented here on site response is part of a broader study of ground motions from subduction earthquakes in CAM. The extension of the NGA-Sub project database (Mazzoni et al., 2022) to include additional events and stations in Mexico is presented in a companion paper (Contreras et al., 2025). The supplemented database provides multiple ground motion recordings at sites across the Valley of Mexico that enabled the present application of non-reference site analyses to measure site response and its uncertainty.

We have summarized the unique geotechnical conditions in the Valley of Mexico and compiled available site parameters (V_S30, T_S, $zLCG$⁠) from site-specific investigations and other literature. These data support the zonation of the lakebed and surrounding regions (Figures 1 and 2) as currently applied in forward applications of site response using a reference site approach (Ordaz et al., 2024; Reinoso and Ordaz, 1999). This approach relies on normalization of ground motion Fourier amplitudes within the lakebed relative to those at a hill zone site (CU).

We have shown that the CU site, commonly used as a reference, has an appreciable level of site response across a broad period range (based on data from the CUP5 station). The site response results presented in this article are broadly similar to previous findings with regard to relative amplification levels between lakebed zones and the period ranges where these amplifications are maximized. However, there are two notable aspects of our findings: (1) amplification levels are higher due to the referencing of amplification to V_ref = 760 m/s in this work, and (2) by examining data from multiple events, we find short period site amplifications exhibit nonlinearity. The use of the non-reference site approach for the derivation of amplifications is particularly significant from a ground motion modeling perspective, because it removes the need to derive single-station reference GMMs for the CU site as has been done in prior work (Arroyo et al., 2024; Jaimes et al., 2015; Jaimes et al., 2006; Ordaz et al., 1994). Doing so substantially increases the level of rigor with which the reference site GMM can be defined.

A nonlinear model conditioned on V_S30 is presented for site response in the Valley of Mexico. The linear portion of the model is able to capture the main trends of site response variations across the lakebed. The linear model includes the option of applying a differential depth parameter derived from $zLCG$ (this option can be turned off by setting δz = 0). While some resonance effects are seen in the data that are not captured, these effects are small compared to what is captured from the V_S30-based model and are not included herein. We acknowledge that models conditioned solely on mHVSR parameters, which were not developed in this research, could be similarly effective as the proposed V_S30-based model. However, a challenge related to developing such models is that mHVSR curves are only available for 23 of the 89 sites considered in this research, and 19 of those 23 sites are in only two zones (II and IIIa). The nonlinear portion of the proposed model is admittedly based upon limited data, but those data demonstrate nonlinearity unambiguously. The development of a more rigorous nonlinear model would require calibrated simulations in which the nonlinear soil behavior is considered, which was beyond the scope of this investigation.

The proposed modeling approach for sites in the Valley of Mexico involves the use of a regionally customized GMM and a subregional site amplification model. The GMM is an adaptation of the Parker et al. (2022) model for CAM, with the modifications being solely related to backarc attenuation effects. The site amplification model has linear and nonlinear components (Equation 2). The V_S30 portion of the linear component is given by Equation 8 with the coefficients in Table 3. The optional modification based on δz is given by Equation 10 with coefficients in Table 3. The nonlinear component is given by Equation 5 with the coefficients in Table 2. Corresponding values for the two components of within-event variability (⁠$ϕS2S$ and $ϕSS$⁠) are given in Figure 15. Coded versions of the subregional site response model presented in this article and the regionally customized GMM presented by Contreras et al. (2025) are provided on the NGA-Subduction project web portal (https://www.risksciences.ucla.edu/nhr3/nga-subduction).

The findings and recommendations presented here provide an opportunity to improve ground motion hazard analyses in Mexico due to the following factors:

1. Ground motion estimation in the Valley of Mexico (Mexico City) and elsewhere in Mexico is developed using consistent approaches.

2. The models are transparent, thus allowing for future research where individual components can more readily be interrogated and improved.

3. Site response is conditioned upon the site-specific V_S30 parameter—this encourages site characterization and avoids ambiguities that can arise for sites near zone boundaries.

4. Nonlinearity in site response is considered.

The next steps are to evaluate the impact of these models on seismic hazards across Mexico and see whether consensus can be developed regarding their introduction into national and regional seismic design standards. Such analyses have not been undertaken in the present work.

The authors are grateful to the NGA-Sub project for providing the database and models from which this project was developed. SSN data was obtained by the Servicio Sismológico Nacional (México); station maintenance, data acquisition, and distribution are possible thanks to its personnel. The views expressed herein are those of the authors and do not necessarily reflect the views of the CTBTO Preparatory Commission. Daniel De La Rosa assisted with data collection for sites in Mexico, particularly geologic maps. We appreciate the helpful feedback from the Associate Editor and two anonymous reviewers, which improved the paper.