## 2 Data and Methodology of Inversion
### Data Testing  
To evaluate the random noise in the data, we performed an odd/even test in which we also evaluated the stability of the 1D model. Then, we compared the inversion results obtained with the same parameters. If the images are very different, it may indicate problems in the data or in the reconstruction algorithm. In our case, the results of the odd/even test for the P- and S-wave velocity anomalies at 15 km depth are presented in Figure 5 and a good correlation of the structures can be observed. The spatial resolution of the model was evaluated by synthetic tests in which we reproduced the same processing conditions of the experimental data. To define the modeling rays, we used the same source-receiver pairs as in the experimental dataset, with the event locations taken from the final solution of the main model. The synthetic model was built by overlaying a set of 3D anomalies on a 1D reference model. Synthetic travel times were computed in the 3D model using a ray-tracing bending algorithm. To simulate realistic conditions, we introduced random noise with predefined standard deviations of 0.1 and 0.2 for the P-wave and S-wave data, respectively. These noise levels were selected to achieve a variance reduction comparable to that of inverting experimental data. The number of  
Figure 4: Distributions of nodes of two parameterization P-wave grids (red dots) with different basic orientations (0 degrees, left and 45 degrees, right) used for the inversions. The background is the normalized ray density in the depth range 12.5-17.5 km. iterations and all control parameters during the synthetic model recovery were identical to those used to compute the main model. To optimize the recovery quality, we adjusted the control parameters during the synthetic tests and then applied these optimal values to the inversion of the experimental data in the main model. This process involved running the synthetic modeling and the experimental data inversion concurrently. Figure 6 shows the results of a checkerboard test in which the synthetic anomalies were represented as 40 \(\times\) 40 \(\times\) 20 km. Within these cubes, alternating anomalies with amplitudes of \(\pm\)5% were introduced. For dVp and dVs, anomalies were assigned opposite signs to enhance the contrast of the Vp/Vs ratio. The inversion of the P- and S-wave data produced a variance reduction of about 55% and 50% for the P- and S-wave data, compared to those of the inversion of the experimental data, respectively, 40% for P and 54% for S, which shows a more accurate model than the original one. The recovered results in Figure 6 are displayed in three depth sections, corresponding to the central regions of the anomalies in different layers. The results show that in the first two layers, the recovered anomaly amplitudes are close to the original values for both P- and S-waves. However, in the deepest layer, the P-wave amplitudes appear slightly smaller. This suggests that a single set of damping parameters may not achieve optimal amplitude recovery across the entire survey volume. Increasing the damping to enhance recovery in deeper sections can result in amplitude loss in shallower regions, potentially obscuring important anomalies. To balance this, the damping value was selected to maintain the stability of anomaly shapes while  
Figure 5: Odd/even events test (black line is the border line and the blue lines are rivers). minimizing noise artifacts. Features smaller than \(\sim\)40 km laterally and \(\sim\)20 km vertically were not considered reliably resolvable in our dataset. To estimate the vertical resolution, we performed a synthetic test with anomalies defined along the vertical profiles. Typically, the vertical resolution is weaker than the horizontal resolution due to the trade-off between source depth and velocity distribution in tomography studies. Figure 7 shows the vertical checkerboard P-wave vertical tests associated with the vertical profiles used to present the main results.