In this study, a 3D physics-based numerical approach, based on the spectral element numerical code SPEED, is used to simulate seismic wave propagation due to a local earthquake in the Mexico City area. The availability of detailed geological, geophysical, geotechnical, and seismological data allowed for the creation of a large-scale (60 km × 60 km in plan, 10 km in depth) heterogeneous 3D numerical model of the Mexico City area, dimensioned to accurately propagate frequencies up to about 1.3 Hz. The results of numerical simulations are validated against the ground motion recordings of the July 17, 2019, Mw3.2 earthquake, with peak ground acceleration exceeding 0.3 g about 1 km away from the epicenter. A good agreement with records is found, quantitatively evaluated through goodness-of-fit checks. Furthermore, for the lake zone, the simulated decay trend of the peak ground velocity with epicentral distance is reasonably close to the observations, for both horizontal and vertical components. In spite of some limitations, the simulations are successful to provide a realistic picture of seismic wave propagation both in the hill and in the lake zones of Mexico City, including the onset of long-duration quasi-monochromatic ground motion in the basin, with strong amplification at low frequencies (between 0.4 and 0.7 Hz). The numerical results also suggest that surface waves, with predominant prograde particle motion at the ground surface and large ellipticities, dominate the wavefield in the lake zone. Based on these positive outcomes, we conclude that this numerical model may be useful for both a better insight into the seismic response of the Valley of Mexico and the simulation of ground motions during larger-magnitude earthquakes, to generate improved seismic damage scenarios in Mexico City.