Receiver functions (RFs) from two permanent stations (CUIG and ZAIG) of the Mexican broadband network were calculated to infer the crustal structure at both sites. The data were stacked to allow uncertainty estimations and inverted using two global optimization algorithms: simulated annealing and genetic algorithms. Strong and systematic azimuthal dependence was observed in the stacked functions at CUIG site, which is located in the Mexican volcanic belt (MVB). Assuming a one-dimensional five-layer model per azimuthal stack, these variations were explained with the same upper crustal structure including a superficial ∼2 km-thick low-velocity layer associated with the MVB and a well-defined Conrad interface around 15 km depth, different depths for the Moho discontinuity, and variations in both the S-wave velocity and Poisson’s ratio for the upper mantle. Our inversion of observed RFs at CUIG and theoretical considerations suggest an eastward crustal thickening of at least 5 km below the Valley of Mexico. We tested the sensitivity of the modeling procedure by applying it to the recordings at ZAIG station. In this case, no azimuthal dependence was noticed in the stacked functions. After a simultaneous inversion of three different azimuthal stacks, a single four-layer model satisfactorily explained all data, thus revealing a flat one-dimensional crustal model below the Mexican Mesa Central, where ZAIG station is located. These results support the resolution and robustness of the applied method at both sites.