## 4 Numerical Experiments  
We performed numerical simulations using the finite-difference method ([PERSON] and [PERSON], 2011) [42] and utilized the resulting data to assess the effectiveness of the proposed method. Figure 3a shows the P-wave velocity model used in our numerical simulation, with a size of 1000 m by 1000 m (among them, the color scale ranges from blue to red, indicating that the P-wave velocity increases from low to high). The upper boundary of the model is set as a free surface, while the left, right, and bottom boundaries are configured as tapered absorbing boundaries to prevent reflections. Considering that seismic interferometry requires cross-correlation calculations between every two receivers, we opted for a larger receiver spacing to reduce computational load. However, to ensure accurate reconstruction of the surface wave virtual shot gather, the receiver spacing cannot be too large. Therefore, the receivers are uniformly distributed with a 5 m interval along the free surface. To compensate for the short simulation time, which makes it difficult to reconstruct effective signals, we increased the number of noise sources. A total of 8000 noise sources are randomly distributed between depths of 10 m and 190 m, with a maximum frequency of 30 Hz. The amplitude of the sources follows a Gaussian distribution, and the duration of each source is 120 s multiplied by a random percentage. The noise sources are triggered at random times during the 120 s simulation period. For the numerical simulation in the elastic medium, we set the S-wave velocity to 60% of the P-wave velocity. First, we select a geophone located at a horizontal distance of 0 m as the virtual source position. As the data are from numerical simulations, we did not conduct normalization and spectral whitening operations. Instead, we directly performed cross-correlation on the vertical component (Z-component) of the passive seismic data. The results are shown in Figure 4a, where clear surface wave dispersion is observed, along with reflected waves at 0.5 s, 0.7 s, and 1.0 s. To verify the reliability of the cross-correlation results, we conducted a similar active numerical simulation using the same model. The same source position is chosen at a horizontal distance of 0 m, with the geophone arrangement identical to that of the passive source. The difference lies in the selection of the source dominant frequency, which is set to 20 Hz. The vertical component shot gather is shown in Figure 4b. Comparedto the virtual shot gather in Figure 4a, consistent surface wave dispersion is observed, along with reflected waves indicated by the arrows. Figure 4: Comparison of passive seismic data cross-correlation results with active source results: (**a**) virtual shot gather reconstructed from passive seismic data in an elastic medium; (**b**) shot gather recorded from direct forward modeling of an active source in an elastic medium. Figure 3: Passive numerical simulation: (**a**) P-wave velocity model; (**b**) raw passive data. We performed the same experiment using the above model in an acoustic medium, and the results are shown in Figure 5. Figure 5a shows the virtual shot seismic record, obtained by cross-correlating the passive simulated data in the acoustic medium, while Figure 5b shows the direct seismic shot gather from the active source simulation in the acoustic medium. In both shot gathers, clear reflections from interfaces at 200 m and 600 m are indicated by the yellow arrows, corresponding to those observed in the elastic medium. Furthermore, a comparison of the horizontal component (X-component) seismic records also confirms that the red arrows in the shot gathers of Figure 4a,b represent converted waves. From the results shown in Figure 4a, it can be observed that, in the virtual shot gather directly obtained by cross-correlating the passive seismic data in an elastic medium, surface waves and body waves are mixed together, making it difficult to effectively process surface and body wave data further.