## 2 Methods
### A Simple Example  
We took a circular cave model (Figure 2a) as an example to illustrate the workflow of MVSS beamforming. The seismic data of the 10 th shot obtained by finite-different time-domain (FDTD) numerical simulations ([PERSON] et al., 2018 [29]) are shown in Figure 2b. The receivers were arranged from 0 m to 30 m with an interval of 0.2 m. There were 151 receivers in total.  
For example, we took point A (16 m, 7 m) as the target point (Figure 2). The delay time was calculated as Equation (A1), and the results are shown in Figure 3. The 151 black dots represent the delay time of each receiver, and the 20 red circles represent the delay time of each shot. The first source shared the same spatial position and delay time with the 36 th receiver. We moved the target point and repeated the delay time calculation using Equation (A1).  
As per the definition of a subarray in Equation (A1), we divided the array of 151 receivers into 76 subarrays, and each subarray was composed of 75 receivers. Taking the first subarray as an example, we calculated the weight of each receiver in the subarray (Figure 4b) according to the DL and MV method using Equation (A13). We multi  
Figure 1: Workflow of MVSS beamforming imaging.  
plied the 75 delayed signals (Figure 4a) with the weight of each receiver (Figure 4b) and then added them together to obtain the amplitudes on the axis line \(x=x_{A}\) (Figure 4c). Two reflection interfaces were observed from the curve. The delay time was reduced by \(t=t_{m}+t_{n}=0.0186\) s to exclude the delay time of point A itself. Before the 0.0186 s was excluded, the delay time is not absolutely correct. If we see point A as the target point, the delay time of the receiver which is above point A should be 0. In delay time calculations according to Equation (A1), the reference point is not point A (16,7) but it is the point on x axis (16,0). This results in the 0.0186 s of delay time of point A itself. So, the 0.0186 s should be excluded from the delay time and the correct time-axis position of point A can be obtained. This can also be found in Equation (A3) as \((t-\tau_{m})\). Then, we repeated the process above with each subarray to obtain the amplitude curve on the axis line \(x=x_{A}\) (Figure 4d). In this certain case, the amplitude of a reflected wave was 0.005 times that of a surface wave in seismic signal. Eventually, the amplitude was 0.05 times that of a surface wave on the superposed curve for all subarrays. Then, we picked the peak amplitude and the corresponding estimated time (0.0093 s) as the amplitude result of point A.  
We repeated the delaying, weighting, and superposition processes with the subarrays at every point in the imaging area to obtain the preliminary image for this shot (Figure 5a). Then, the CF (Figure 5b) was obtained according to Equation (A15). The CF matrix was used as another weight to enhance the SNR by multiplying it with the preliminary image. In this step, surface wave artifacts were suppressed because the adopted methods were developed for reflected waves to only amplify the reflected signals with focus. Other interferences were suppressed for the same reason. Then, we obtained the final image of the 10 th shot.  
Figure 3: Time delay of each receiver and source.  
Figure 2: (a) Circular cave model. (b) Signal of the 10 th shot with 20 stacked traces.  
The time-domain image was converted into a depth-domain image according to the background velocity. The same processing scheme, including beamforming imaging and time-depth conversion, was applied to the remaining shots. The final image was a stacked result of all shots (Figure 6).