The raw 3D data processing scheme is similar to that proposed by Barone et al.^{17}. Travel times for different frequencies are directly extracted from the phase of the surface wave signal, taking advantage of fine sampling in both spatial directions. However, the phase distribution must be representative of a single mode of propagation (i.e. the fundamental mode of the Rayleigh wave). More generally, any coherent noise (e.g. higher modes, vibrations from stationary noise sources, backscatter, etc.) should be removed from the signal, as it disturbs both phase and amplitude variations with offset, generating periodic patterns, as explained in detail by Baron et al.^{29}.

Preliminary analyzes on the recorded data set are needed to investigate signal characteristics, including the frequency content of Rayleigh waves, the presence of higher-order modes and/or coherent noise, and the degree of heterogeneity in the two study areas, and to identify the LMO velocities used for correction in 3D processing (see Supplementary Text S1). This preliminary work consists of the analysis of f–k spectra calculated on several 2D lines extracted from the 3D dataset. The results showed significant backscattered energy and at least one higher-order propagation mode with non-negligible energy. The frequency band of interest to the fundamental mode extends from 10 to 50 Hz, making this interval suitable for our further analysis. The dispersion curves plotted along the different lines show a high degree of dispersion, revealing very heterogeneous conditions. From this analysis, we can also appreciate the overall lower speeds inside the amphitheater compared to the outside area, which has been clearly illustrated in the “Results” section.

The 3D treatment sequence described below refers to a single plane and a single frequency and consists of three main steps. First, we apply a linear displacement correction (LMO) to the raw traces to (i) reduce the phase jumps related to the 2(pi) periodicity of the phase, and (ii) shift energy from higher modes to the negative wavenumber quadrant of the f–k spectrum. Second, we perform pseudo-2D f–k filtering on azimuth sectors of equal width. Third, we apply a 2D phase unwinding scheme to compensate for residual phase jumps and remove the LMO. The frequencies used in this analysis are 10.00 Hz, 10.98 Hz, 12.16 Hz, 13.64 Hz, 15.52 Hz, 18.00 Hz, 21.43 Hz, 26.47 Hz, 34.62 Hz and 50 Hz. Such an inversely regular frequency sampling results in a quasi-regular sampling in depth. Indeed, the depth of maximum sensitivity of surface waves is proportional to the wavelength, which is inversely proportional to the frequency. The phase maps obtained are finally converted into relative travel time maps (referred to the receiver closest to the source) as:

$$begin{aligned} Delta t = frac{phi – phi _0}{2 pi f}, end{aligned}$$

(1)

where (Delta t) is the relative travel time, (phi) is the phase recorded at any receiver, (phi _0) is the phase recorded at the receiver closest to the source and *F* is the frequency. Two separate processing steps were performed for Grid 1 and Grid 2, using local coordinates. Only close-up shots around each patch, covering all azimuths, were selected: shots 1 to 19 were included in the Grid 1 processing while shots 2 to 7 and 20 to 38 were used for the grid 2. A detailed description of the different processing steps is given in the supplementary information (Supplementary Texts S2 to S4).

Travel time maps for one frequency and all shots are entered into an Eikonal tomography scheme to extract phase velocity maps. Eikonal tomography derives phase slownesses at a frequency from the magnitude of the travel time gradient^{9}. This operation is repeated for each shot point, and the phase velocities are finally calculated as the inverse of the average phase slownesses on the different shots. This method of tomography is very fast because it does not involve true inversion (see discussion in Barone et al.^{17}). Moreover, the regular receiver geometry of our dataset is perfectly suited for this method, as it renders unnecessary the most critical and user-dependent step of Eikonal tomography, which is the 2D interpolation of the travel time on a regular grid. Another relevant aspect of Eikonal tomography is that it does not require absolute travel times (referenced to a source/virtual source position), but it only takes into account travel time differences between adjacent receivers (for this reason, we used relative travel time maps). Finally, this method makes it possible to extract standard deviations of speed, which represent an estimate of the error for the speeds obtained. On the other hand, Eikonal tomography is very sensitive to outliers^{17}. For this reason, we performed a statistical analysis to identify and remove outliers from single-shot phase velocity maps: velocity values outside the range of the spatial mean plus/minus three times the standard deviation were rejected.

The analysis described so far has mainly focused on the phase of the signal, with the aim of extracting velocity information. However, most tomography approaches, including Eikonal tomography, are based on the high-frequency approximation^{30}. For this reason, the achievable lateral resolution is inherently limited by the signal wavelength which, for low frequencies, is much greater than the sensor spacing. Signal amplitude analysis can help detect sudden changes in lateral velocity caused by small objects or discontinuities^{31}. Several methods for measuring amplitude variations are available. We focused on the autospectrum method^{32}which analyzes the autospectral density at each frequency *g*(*F*), defined as:

$$begin{aligned} G(f) = { Im[Y(f)] }^2 + { D[Y(f)] }^2 = { A(f) }^2, end{aligned}$$

(2)

where *Yes* and *A* are respectively the complex spectrum and the amplitude spectrum of a seismic record. The autospectral density is a measurement of the surface wave energy for a certain frequency: significant spatial variations of this parameter could indicate the presence of a scatterer (negative variation) or of an amplifying zone (positive variation)^{31}. For this reason, this study focuses on the analysis of autospectral gradient maps. The procedure followed here is to retrieve autospectral density maps for all planes, which are normalized by their maximum value. Then the autospectral gradient is calculated and an average is calculated from the gradient amplitude maps of different planes. From this analysis we obtain a robust estimate of the spatial variations of energy, even if the positive/negative sign of the anomalies is not preserved. Traces included in the analysis should not be filtered, as they should include backscatter. In addition, geometric spreading effects must be removed before the Fourier transform, and the quasi-shift region must be excluded as for phase velocity analysis.

Depth inversion is required to produce a quasi-3D shear wave velocity model of the near surface. In order to ensure the stability of the inversion, the phase velocity maps were initially smoothed with a Gaussian filter. The length of the filter depends on the frequency: we used (lambda /2) as filter length, with (lambda) being the average wavelength at a specific frequency. In doing so, low-frequency maps, characterized by lower spatial resolution, undergo stronger filtering, while high-frequency maps better preserve their spatial variability. Local dispersion curves were obtained from the superposition of smoothed phase velocity maps at different frequencies. The dispersion curves at each location were inverted independently, with no lateral constraint. For the inversion, we used the dinverdc module in the GEOPSY Dinver framework^{33}which uses the neighborhood algorithm^{34}. Full details on the model space parametrization are given in the Supplementary Information (Supplementary Text S8).