## 3 Acquisition and Processing of GPR Data  
For the geophysical measurement acquisition, the GSSI GPR system SIR-3000 with a 400 MHz monostatic shielded antenna was employed (Figure 2). The GPR measurements were conducted in front of Quetzalcoatl Temple covering a surface of 34 \(\times\) 20 m. For 3D acquisition, 34 m-long parallel study lines were surveyed along a south-north orientation with lines spaced every 0.25 m. The trace interval and time window were 0.02 m and 100 ns, respectively. For the initial processing of the GPR data, the software package of Radlab (written in MATLAB) developed at the Institut Terre et Environnement de Strasbourg and dedicated in-house MATLABcodes developed at the Applied Geophysics Lab. of the Technical University of Crete were utilized. The processing sequence involved time-zero shift, low-frequency removal or DC filtering (dewow), a background filter of 200 traces, band-pass filtering of 200-600 MHz, and, finally, amplitude gain (linear and/or exponential). We determined an average electromagnetic (EM) wave velocity of 0.095 m/ns by analyzing different diffraction hyperbolas, which indicated a soil with moderate humidity. For the three-dimensional (3D) presentation of GPR data, we used ReflexW software. Figure 1: (**a**) Stratification planes, joints, and fractures enlarged by the hydrochemical dissolution process; (**b**) cave with anthropic modifications in the ceiling for astronomical observation. The initially processed GPR sections' data were input into a multipath summation sequence approach, as described by [PERSON] et al. [19], and involved the following processing steps for each of the GPR sections: (1) Application of constant-velocity F-K migrations; (2) estimation and assignment of a weighting factor for each migrated section, with which it contributes to the final stack; (3) weighted stacking of the time migrated sections; and (4) post-multipath summation time-varying spectral whitening of the stacked section. Figure 2: In front of Quetzalcoatl Temple, 81 profiles spaced 0.25 m apart were acquired with a 400 MHz antenna. The orientation of the profiles is S–N (**a**). Aerial view of the Xochicalco archaeological site showing the surveyed area in front of Quetzalcoatl Temple, where 81 profiles spaced 0.25 m apart were acquired with a 400 MHz antenna; the orientation of the profiles is S–N (**b**). The velocity range of the constant velocity migrations for the first processing step was set to 0.04-0.24 m/ns. The constant velocity value utilized for each migration was altered every 0.005 m/ns, which concluded with 33 constant velocity sections. Each constant-velocity migrated section was evaluated regarding the focusing of the diffractions. This was implemented by the estimation of the local slopes using the local plane-wave equation [28]. \[\frac{\partial P}{\partial x}+\sigma\frac{\partial P}{\partial t}=0 \tag{1}\]  
where \(P\) is the wave field, \(x\) and \(t\) are the space and time variables, and \(\sigma\) denotes the local slope. Sections with unfocused diffractions were expected to have higher standard deviation (std) values of the local slope distribution in comparison to the efficiently migrated sections, as the former included much more local slope values that are far from zero due to the diffraction branches energy not gathered to the apices. Let the reciprocal value of the std of the local slope distribution, \(s_{i}\), be \(s_{i}^{\prime}\).