Investigation of the Time-Lapse Changes with the DAS Borehole Data at the Brady Geothermal Field Using Deconvolution Interferometry

: The distributed acoustic sensing (DAS) has great potential for monitoring natural-resource reservoirs and borehole conditions. However, the large volume of data and complicated wavefield add challenges to processing and interpretation. In this study, we demonstrate that seismic interferometry based on deconvolution is a convenient tool for analyzing this complicated wavefield. We also show the limitation of this technique that it still requires good coupling to extract the signal of interest. We extract coherent wave from the observation of a borehole DAS system at the Brady geothermal field in Nevada. The extracted waves are cable or casing ringing that reverberate within a depth interval. These ringing phenomena are frequently observed in the vertical borehole DAS data. The deconvolution method allows us to examine the wavefield at different boundary conditions and separate the direct waves and the multiples. With these benefits, we can interpret the wavefields using a simple 1D string model and monitor its temporal changes. The velocity of this wave varies with depth, observation time, temperature, and pressure. We find the velocity is sensitive to disturbances in the borehole related to increasing operation intensity. The velocity decreases with rising temperature. The reverberation can be decomposed into distinct vibration modes in the spectrum. We find that the wave is dispersive, and the fundamental mode propagate with a large velocity. This interferometry method can be useful for monitoring borehole conditions or reservoir property changes using densely-sampled DAS data.


Introduction
The fiber-based sensors have been applied in the oil and gas industry for borehole monitoring since early 90's [1]. Since then, the distributed temperature sensor (DTS) has been routinely deployed for monitoring well temperatures. The distributed acoustic sensor (DAS) has gained popularity in seismology more recently. The DAS measures strain rate, and thus, records the seismic wavefield like a string of one-component geophones. The advances in fiber material and computer technologies allow us to obtain data with higher quality and analyze them with array processing techniques.
The DAS has been used in boreholes environments for a variety of applications. These include flow monitoring [2][3][4], wellbore diagnostics [2,4,5], vertical seismic profiling [VSP; [6][7][8][9], hydraulic fracture characterization [10,11], and microseismicity detection [12]. The DAS is suitable for borehole monitoring for several reasons [13,14]: First, the DAS fiber has higher endurance in high temperature, high pressure, and corrosive environments compared to geophones. Second, it provides a dense 1D receiver arrays along the wellbore. Finally, the cost of DAS borehole deployment is relatively low, although the interrogator and the data storage can be expensive. Once installed, the fiber can be left in the well for long-term monitoring without changing locations. This resolves one of the main difficulties for conventional 4D (3D and time) surveys.
A big challenge of analyzing the DAS wavefields is they are often complicated, especially in the borehole environment. Transient borehole processes such as fluid flows and operation activities cause disturbances in the borehole. Optical noises from the DAS interrogator create noise that occurs simultaneously and same amplitude on all sensor channels [6]. For the DAS data we analyze in this study, cable and casing ringing populate a large portion of the data [15]. The ringing is a common phenomenon for DAS in a vertical borehole due to poor coupling between the DAS cable and the casing, or between the casing and the formation [16][17][18][19][20]. It appears as bouncing waves that reverberate within a depth interval. For VSP applications, the ringing is a noise that analysts want to get rid of [21]. Here, we treat these ringing waves as signals and analyze their time-lapse changes. This allows us to interpret the dominant energy sources in the system and understand if the cable and the casing are sensitive to certain processes.
The DAS data we analyze are from a vertical borehole at the Brady geothermal field in Nevada. They were obtained during the PoroTomo project [22,23]. The PoroTomo project was a four-week experiment conducted during March 2016, in which the team performed vibroseis experiments under varying pumping operations and collected a variety of geophysical data including surface DAS (DASH), borehole DAS (DASV), nodal geophones, InSAR, GPS, pressure, and temperature (DTS) data. The DASV data were available from  Previous studies have analyzed the DASV, DTS, and pressure data. Patterson et al. [24] and Patterson [25] analyzed the borehole DTS and pressure data at different stages of operations. Trainor-Guitton et al. [26] imaged features on two nearby steeply dipping faults using a portion of the DASV data. Miller et al. [15] investigated the DASV data to find the signatures of earthquakes, vibroseis sweeps, and responses to different borehole processes. In addition, they suggest the reverberations on the upper half of the DASV is due to ringing of the casing and the DAS cable. We follow their results and further investigate the time-lapse changes of these reverberations.
We use deconvolution seismic interferometry to extract coherent signals along the 1D receivers of the borehole DASV array. The coherent signals are governed by the same wave physics (i.e., wave equation) [27]. Thus, we can understand the property of the structure by examining this wave. This deconvolution method is useful because it modifies the boundary conditions [27][28][29]. Thus, we can convert the wavefield to a favored boundary condition for interpretation. For example, Snieder and Safak [30], Nakata et al. [29], and Nakata and Snieder [31] used this to isolate the ground coupling effect and analyzed the vibration modes of a building. Sawazaki et al. [32], Yamada et al. [33], Nakata and Snieder [34], and Bonilla et al. [35] applied similar methods to obtain near-surface velocity changes in different time scales. In this study, we use this deconvolution method to help us examine the wavefield. It allows us to interpret the wavefield using a simple model. Furthermore, it separates the direct waves and the multiples and simplifies the wavefields. This makes time-lapse monitoring easier to implement.
In the following, we first introduce the Brady DASV data (Section 2) and the deconvolution interferometry method (Section 3.1). We focus on analyzing the ringing signals on the shallower part of the borehole as their energies are consistent. In Section 3.2, we show the deconvolved wavefields and how they can be explained by our proposed models. Detailed parameters for modeling are given in Appendix A. In Section 3.3, we analyze the velocity variations of the signal versus measured depth, observation time, temperature, and 56A-1 (P sensor) Figure 1. Locations of the target boreholes in the PoroTomo experiment. The survey was at the Brady geothermal field in Nevada, USA (black cross in the inset). The red star is the borehole with DASV and DTS (well 56-1). The green dot is the borehole with the pressure (P) sensor (well 56A-1). The blue triangles are locations of the vibroseis shots. The gray lines are the DASH cable on the surface. We use DASV, DTS, and pressure data in this study.
pressure. In Section 3.4, we apply a normal-mode analysis to the vibration modes of the waves. We also show the signals we extracted from waves propagating in the formation in Appendix B. For these signals, we cannot get precise velocity changes due to poor coupling. But we note that if coupling was better, similar techniques can be applied to these signals for analyzing temporal changes in the formation.

Data
We focus on the DASV, the DTS temperature, and the pressure data from the PoroTomo project [22,23]. Figure 1 shows the location of the wells relative to the entire DASH array and vibroseis shots on the surface. The DASV and DTS fibers are co-located in well 56-1 (the red star) that is about 380 m deep. The DASV fiber is single-mode and the DTS fibre is multi-mode. Both fibers are high temperature acrylate-coated that are tested to be resilient up to 150 • C. For resistance, the fibers are protected by stainless steel double tubing. The DASV system has 384 channels with channel spacing of approximately 1 m. The gauge length is 10 m. The sampling rate is 1000/s. The unit of the DAS raw data is radian/millisecond per gauge length. The total DASV data size is 981 GB stored in SEG-Y format. The DTS system has channel interval of 0.126 m and the sampling interval is 62 s. The pressure sensor (P sensor) is located at a nearby well 56A-1 (the green dot). The pressure sensor is at an elevation corresponding to channel 219 of the DASV system (i.e., measured depth = 219 m). The sampling interval of the pressure sensor is 60 s. The two wells are around 100 m away from each other. Their wellheads are at about the same sea level (1230 m). Patterson [25] suggested the two wells are hydraulically connected based on simultaneous responses between the DTS and the pressure sensor. Hence, we assume the pressure measurements can represent the co-located pressure changes with the DASV and DTS at measured depth of 219 m. Figure 2 shows the evolution of the pressure, temperature, and an overview of the DASV DC values and its Root-Mean-Square (RMS) amplitudes. We focus on the eight days (Mar [18][19][20][21][22][23][24][25][26] where the DASV was actively recording. Initially during Mar 18, the pressure drops drastically due to resuming operation after a shutdown period (yellow to blue shade in Figure 2a). Then, the pressure increases slowly due to increasing injection,   until resuming to normal operation on Mar 24 (blue to green shade). The sudden pressure rise at the end of Mar 25 is due to a plant shutdown [36]. The temperature increases with depth with a heat deficit below 320 m due to geothermal explorations (Figure 2b; Miller et al. 15). The lower temperature in early Mar 18 is due to cool water treatments before cable installation. Figure 2c,2d show the DASV data contains many disturbances under these changing pressure and temperature conditions. Patterson et al. [24] and Miller et al. [15] investigated these events. Here, we analyze the changes of extracted waves in the deconvolved wavefields.

Review of deconvolution interferometry
We use deconvolution interferometry to extract coherent waves from the data. The receiver used for deconvolution is a "virtual source". The deconvolution operation modifies the boundary conditions of the wavefield depending on the virtual source [27][28][29]. These boundary conditions include coupling, attenuation, and damping at the boundaries that obscure the pure response of the system. By examining the wavefields that satisfy different boundary conditions, we can potentially separate these unwanted effects. For time-lapse monitoring, this allows us to track the pure response of the structure. Snieder and Safak [30] and Nakata et al. [29] used this deconvolution method to retrieve the vibration modes of the building with receivers deployed along the building floors. Nakata and Snieder [34] monitored monthly and annual changes in shear-wave velocity using near-surface and borehole sensors. Sawazaki et al. [32], Yamada et al. [33], Nakata and Snieder [37] and Bonilla et al. [35] analyzed the near-surface velocity changes during earthquake strong ground motions. Here, we use it to analyze the reverberations that are commonly observed for DAS in a vertical borehole. We also show its potential for time-lapse borehole condition and reservoir monitoring.
The deconvolved wavefield D in the frequency domain is [29] where z is the depth of each channel, z a is the depth of the virtual source channel, ω is the angular frequency, and * denotes the complex conjugate. The deconvolution operation in the frequency domain is the division of the data recorded at each depth (U z (ω)) by the data recorded by the receiver that is used as the virtual source (U z a (ω)). The instability in Equation 1 comes from the division, and we stabilize it with a water level ε = 0.5% that scales with the average power spectrum (⟨|U z a | 2 ⟩ in Equation 2). Given a virtual source channel, we calculate the deconvolved wavefield using Equation 2 for the entire 1D array and stack the resulting wavefields over a time span to improve the signal to noise ratio. Figure 3 summarizes the workflow for calculating the deconvolved wavefield.
We use two sets of time windows to calculate the deconvolved wavefields. In Section 3.2, we use 30 minutes time window, 50% time overlap (time step=15 minutes), and then stack the deconvolved wavefields over 3 hours (time span=3 hours) to enhance the signal to noise ratio. In Section 3.3, we use 1 minute time window, 50% time overlap (time step=30 seconds), and then stack them over 1 hour (time span= 1 hours). Since the deconvolution is conducted in the frequency domain (Equation 2), we demean, detrend, and taper (10% on both sides) the raw data at each window before Fourier transform. For simplicity, we omit the "stacked" term and call the final retrieved wavefield as "deconvolved wavefield" for the rest of this paper. Figure 4 shows the deconvolved wavefields in the upper half of the borehole between 0-200 m (the top panels in subfigures a-f). We obtain strong reverberating signals that bounce between 10 and 165 m. They are only present when the virtual source is within the same depth interval. When we put the virtual source below 200 m, the reverberations almost disappear. This suggests that these waves are restricted in this depth interval. Figure 4a-4c are wavefields at the same time window but deconvolved with different virtual sources marked by the red dash lines. They show distinct differences. The waves in the deconvolved wavefield are coherent energies assuming they are excited at the virtual source.

Deconvolved wavefields
To explain the observed deconvolved wavefields, we use a simple string model and derive its mathematical notation. Figure 5 shows the sketch of the model (model 1). This string model has two reflectors (R 1 and R 2 ) on the top and the bottom as boundaries and two sources (S 1 and S 2 ) at those boundaries. This model can also represent the case when we have sources that are further away from the end points outside of this receiver line [34]. Hence, one should consider S 1 and S 2 as the incoming waves from the top and the bottom to the system. The wavefield of a single source can be expressed by a sum of a power series as shown in Nakata et al. [29]. Expanding from this, the wavefield of two sources with configuration in Figure 5 are superposition of their individual wavefields. That is   Figure  4. The model has a string with a line of receivers on it (blue line) bounded by two reflectors at z = 0 (R 1 ) and z = H (R 2 ). The two sources are located at z = 0 (S 1 ) and z = H (S 2 ) (magenta stars).
where z is depth, ω is the angular frequency, i = the imaginary number, k is the wave number, and H is the length of the structure, γ is the attenuation factor where γ = 1 2Q [38], Q is the quality factor, S 1 and S 2 denote the spectrum of the two sources terms and R 1 and R 2 are the reflection coefficients of the top and bottom reflectors, respectively. In the nominator, e z(ik−γ|k|) and R 2 e (2H−z)(ik−γ|k|) are the direct wave and the first reflection for S 1 , while e (H−z)(ik−γ|k|) and R 1 e (H+z)(ik−γ|k|) are those for S 2 . Their amplitudes are scaled by the attenuation terms that involve γ. The R 1 R 2 e 2H(ik−γ|k|) term in the denominator is the common ratio in the power series representing higher-order reverberations between two reflectors.
We simulate the deconvolved wavefields using Equations 2 and 3 and compare them with the observed deconvolved wavefields (Figures 4a-4f). After a series of parameter tests shown in Appendix A.1, we set all sources terms to be mutually uncorrelated with their cross-correlation coefficient cc = 0.01. This choice is because correlated source would generate simultaneous direct waves from the virtual source, as shown in Appendix A.1, which we do not observe in the wavefield. Other parameters used are Q = 500, ω/k = 4600 m/s and ε = 0.0001%, R 1 = R 2 = 0.9 for 4a-4d, and R 1 = R 2 = 0.5 for 4e,4f. These choices are based on the low attenuation across depth, apparent velocity of the signal, and high reflectivity at the boundaries in the observed data. Figure 4 shows that we can reproduce the wavefields using model 1. In Figure 4a-4c, the same time window is examined using three different virtual sources (the red lines). The dominant waves exhibit symmetry between causal and acausal times (i.e., left and right to the blue lines) regardless of the virtual source. To achieve this symmetry regardless of the virtual source, the two sources in model 1 need to have comparable amplitudes (Appendix A.2). Figure 4d-4f show that the model can also reproduce the three special cases observed in the data. In Figure 4d, the multiples are much weaker than in Figure 4a. We reproduce it by letting the source S 1 at the depth of the virtual source, be much weaker than S 2 , the source at the opposite end of the interval. In Figure 4e,4f, the dominant waves are asymmetric with only causal waves. We reproduce this by minimizing the amplitude of one of the sources and using the main source as the virtual source. Hence, we can explain these special cases with unequal amplitudes of S 1 and S 2 . In Appendix A.2, we analyze the effect of varying relative source amplitudes. Some observed deconvolved wavefields suggest a more complicated model ( Figure  6). The observed wavefield in the top panel of Figure 6a shows a reflector at near 90-100 m. We reproduce this wavefield using model 2 shown in Figure 6b. In model 2, we add an additional source S 1a co-located with S 1 at z = 0 (the dark blue star). This additional source generates waves that propagate between z = 0 and z = H/2 (the dark blue line). A reflector R 3 at z = H/2 acts as a lower boundary for this wave. The high RMS amplitude near 90-100 m in Figure 2d supports this model.

Time-lapse changes of wave velocities
In this section, we analyze the velocity evolution using the extracted wave. The deconvolved wavefields are calculated with 1 minute time window, 50% overlap, and stacked over 1 hour. We calculate deconvolved wavefields with the virtual source at 180 m and measure the arrival times by picking the peaks of the up-going direct wave within 70-120 m. The signals are the most consistent over the eight days between this depth range. We calculate the velocities for a channel by dividing the measured travel length (between the source channel and the target channel) by the picked arrival time. In Figure 7, we plot the estimated velocities against measured depth, observation time, temperature, and pressure. Each gray dot is a velocity measurement at a channel. In general, the velocity of this signal is at around 3600-5000 m/s. This velocity range is much higher than that of the local formation (V p =1000-2500 m/s; Parker et al. 39,Thurber et al. 40). They are closer to the compressional velocity of steel (5000-5250 m/s; Haynes 41). The waves are likely propagate in the stainless steel DAS cable jacket or the steel well-casing as Miller et al. [15] suggested. Figure 7a shows the velocities across measured depth of 70-120 m. The velocities show a slight decreasing trend of -6.6 m/s per meter, which reflects the negative temperaturevelocity dependency in Figure 7c since the temperature increases with depth at this depth range (Figure 2b). The velocity variations (the width of the blue shade) are larger near 72 m and 100 m. This larger variation potentially indicate poor coupling of the DAS cable, or it can be related to the complicated structure and additional source observed in Figure 6.  of velocity during Mar 18-19 is likely associated with disturbances in the borehole. This disturbance is caused by depressurization boiling due to the initial pressure drop related to increasing operation intensity [25] (Figure 2a). During this time, the DAS data also have a high DC level (Figure 2c). Figure 7c shows the velocity decreases with increasing temperature with a slope of -17.1 m/s/ • C. This temperature sensitivity is much higher than that measured in the lab for pure steel material (-0.5 m/s/ • C; Mott 42; Droney et al. 43). We have two possible explanations for this. If the waves propagate in the DAS cable jacket, then this higher sensitivity might suggest the cable, or the fiber inside being subjected to the high temperature. We note that the DAS fibre is rated to 150 • C while the highest temperature in the borehole is beyond 160 • C (Figure 2b). On the other hand, if the waves propagate in the well-casing, it suggests the casing might have higher sensitivity to temperature. In Figure 7d, we do not observe obvious relation between velocity and pressure due to lacking samples at higher pressure.

Normal-mode analysis
The deconvolved wavefield of a vibrating 1D structure can be written as the summation of normal modes [29,30]. This is observed in our results. Figure 8 shows the amplitude spectrum of one of the deconvolved wavefields we used for time-lapse analysis. The normal modes of the signal are clearly decomposed from 10 Hz to over 200 Hz. The frequency interval between different modes is about 18 Hz and consistent over all modes as expected.
The system has closed boundaries on both ends. The top boundary is due to the free surface that behaves as closed boundary for P-wave multiples. The bottom boundary is because of the deconvolution modifying the boundary condition to clamped boundary (a delta function) at the virtual source [27]. For this system, the wavelength of mode m is [44] where H is the length of the system. Hence, the phase velocity for mode m is where f m is the mode frequency. We estimate f m and H in the hourly stacked amplitude spectrum at 6 or 7 am on each day. This is the time with relatively high signal to noise ratio. We focusing on the 2nd (∼38 Hz), the 3rd (∼55 Hz), and the 4th (∼71 Hz) modes, since these three modes are the most significant. We pick f m at the peak amplitude of each mode. We estimate H by picking the starting and ending depths of the mode. Then, we calculate the phase velocity using Equation 5.
The temporal variations of mode frequency, system length, and phase velocity are shown in Figure 9. The mode frequencies increase slightly with time while the system length decreases with time ( Figure 9a). The velocities estimated using higher modes are lower, suggesting a negative frequency-velocity dependency. Hence, the velocities are dispersive [45]. In general, the three modes show similar trends: the velocities increase for the first few days before Mar 20, and then, they continuously decrease until the end of the analysis period.

Discussions
We extract coherent waves in the borehole DASV data using deconvolution seismic interferometry (Figure 4). The extracted waves are the ringing of the DAS cable and the well-casing based on the velocity (Figure 7 and 9b). They are caused by poor coupling between the cable to the well, or between the well and the formation. By using different virtual sources, we examine the wavefields that satisfy different boundary conditions. A simple model with two sources and two reflectors (model 1) can explain the deconvolved wavefields ( Figure 4). Some wavefields exhibit more complexity and suggest a more sophisticated model (model 2; Figure 6). We use numerical simulations to qualitatively reproduce the direct waves and the multiples in the deconvolved wavefields. In model 1, the reflectors are associated with free surface and potential casing defects [15,25]. In model 2, the added sources and the wave trapped in the upper half of the system (the dark blue line in Figure 6b m of the well [15]. In fact, the actual conditions might be even more complicated and we note that the solutions are not unique. Therefore, this model may not be appropriate if one wants to analyze the inverted full-wavefield. But this model might be the simplest that can explain our extracted waves well.
One important feature of the wavefields in Figure 4a-4c is the symmetry. According to Nakata and Snieder [31], to have symmetry between the causal and acausal times for all virtual sources in this model, we must have more than one source. In Appendix A.2 and A.3, we find the asymmetry is produced by uneven amplitudes of the sources, or uneven reflectivities of the reflectors. The effect of the former on asymmetry is more dominant than the later. We reproduce the symmetry in simulations by having two sources with comparable amplitudes and reflectors with identical reflectivities each other.
The main energy sources in this system are borehole processes, surface operations, and traffic noises. The relative amplitudes of these sources change over time, resulting in different observed cases of deconvolved wavefields in Figures 4. This noise variation can generate the asymmetry discussed above. The borehole processes include depressurization boiling and fluid exchange activities at potential casing defects [15,25]. These processes are the most intense during the initial pressure drop. Hence, the deconvolved wavefield shows strong upgoing waves during this time (Figures 4e). The surface operations include site activities and vibroseismic experiments that were conducted 10 hours a day. During these vibroseis experiments, we observe strong down going waves (Figures 4f). In addition, the interstate highway on the north-western side of the survey region provides traffic noises as a general energy source of the extracted wave [46].
We separate the direct waves and the multiples of the ringing signals by using the base of this system as the virtual source. We track the velocity variations of the direct waves. The dense spatial sampling of DAS allows us to observe the trend of velocity variations along depth, time, temperature, and pressure within a 50 m interval (Figure 7). The velocity variations provide potential information for coupling conditions along the cable and parameters that these ringing signals are sensitive to. In Figure 7a, the depth with large velocity variations might suggest poor coupling, or presence of energy source or complex structure. The later is also supported by model 2 which explains occasional variation of the wavefields ( Figure 6). As discussed in Section 3.3, the temporal correlation between the high velocity and the time of large borehole disturbances suggests that the ringings are sensitive to these disturbances (Figure 7b). Finally, in Figure 7c, the decreasing velocity with temperature indicates the DAS cable or the casing potentially being sensitive to high temperature. Hence, by monitoring the extracted waves, we can gain information of the medium that the waves propagate in.
The velocities estimated by picking arrival times on the propagating-wave (Figures 7) are slightly slower than the normal-mode method (Figure 9). This is because the normalmode analysis is done in the lower frequency modes that have higher velocities (Figure 9b) whereas the propagating waves contain all frequencies. The frequency-dependent velocities from the normal-mode analysis are potentially useful to obtain attenuation and structures at different distances from the well. However, in this case, since the coupling (either between the DAS cable to the casing or between the casing to the formation) was poor, the dispersion relation is less sensitive to the structure. Instead, the negative frequency-velocity relation might be caused by the casing and fluid in the borehole, but we need a further experiment to understand the dispersion of the waves.
We note that this deconvolution method can be useful for monitoring changes in the reservoir. In Appendix B, we show the signals we extracted on the lower portion of the DASV cable (below 200 m). We are able to obtain signals during some of the vibroseis experiments. However, the poor signal to noise ratio prevents us from analyzing the Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 29 December 2021 doi:10.20944/preprints202105.0014.v4 time-lapse changes with good precision. If the coupling was better, the signal to noise ratio would be improved and we could have signals outside of the vibroseismic experiment times. Then, we can apply the similar time-lapse velocity analysis on the obtained signal and infer for reservoir properties changes.

Conclusion
We use deconvolution seismic interferometry to analyze the reverberations in the distributed acoustic sensing (DAS) in the borehole. This method is useful for understanding complicated wavefields. The waves we extracted on the shallower part of the borehole is the cable and casing ringing. By examining the wavefield at different boundary conditions, we can qualitatively interpret the system using a simple 1D string model. An important observation is the symmetry of the wavefield. The keys to explain our observations are the source correlations, relative source amplitudes, and reflection coefficients in this system. The deconvolution method also allows us to separate the direct waves and the multiples and track the velocity changes of the direct waves over time. The velocity experiences a rise during the initial pressure drop that was associated with increasing operation intensity. The velocity decreases with increasing temperature and depth. The velocity sensitivity to temperature is higher in our results than that for pure steel measured in the lab. This suggests the DAS cable or the well-casing potentially being affected by the high temperature. The technique proposed here can be applied to many different borehole DAS applications. These applications include diagnosis of the condition of the casing structure and monitoring changes of reservoir properties. For the later, we need better coupling than simply friction in the vertical borehole to obtain more energy from the formation. In this section, we simulate the deconvolved wavefields with varying parameters in model 1. We base on the analysis here to choose the parameters in Section 3.2.
Appendix A.1. The effect of source correlation We vary the correlation coefficient between S 1 and S 2 from 0.01 (uncorrelated), 0.5 (partially correlated), to 0.99 (highly correlated). To generate synthetic data with certain degree of correlation, we first generate random, normalized data time-series and put them in rows to form matrix A. We build a covariance matrix R with the desired correlation coefficient (cc) on the non-diagonals and 1 on the diagonals. Then, we use the Cholesky decomposition to calculate matrix C such that CC T = R. Multiplying A with C gives a new matrix where the cc between each row is as desired. Figure A1 shows the simulation results when the cc equals 0.01, 0.5, and 0.99. When the two sources are not correlated ( Figure A1a), only the virtual source emits waves from time zero. When the two sources are correlated ( Figures A1b,A1c), the correlated source emits another set of waves in addition to that from the virtual source. The higher the correlation, the larger the amplitudes of those simultaneous direct waves. Since we do not observe this simultaneous direct waves in the data, we set cc = 0.01 in the simulations in We vary the relative amplitude of the two source (|S 1 |/|S 2 |) from 0.1, 1, to 10. When |S 1 |/|S 2 | = 1 (Figure A2b), the relative amplitudes on the causal and acausal axes are wellmatched regardless of the depth of the virtual source. The wavefields are symmetric. When |S 1 |/|S 2 | = 0.1 ( Figure A2a), for channels above the virtual source (the red dashed lines), the waves at causal times have larger amplitude, whereas for channels below the virtual source, the waves at acausal times have larger amplitudes. Vice versa, when |S 1 |/|S 2 | = 10 ( Figure A2c), the patterns reverse.
When one of the sources is dominant (Figures A2a,A2c), the deconvolved wavefields approach the one-source cases. This is predicted by the equations. Based on Equations 1 and 3, the deconvolved wavefield using virtual source at z a (0 ≤ z a ≤ H) can be written as = (e (z−z a )(ik−γ|k|) + R 2 e (2H−z−z a )(ik−γ|k|) )+ Similarly, in the case when |S 2 |/|S 1 | ≈ 0 and z a = 0, Equation A2 becomes the infinite series of Equation 9 in Nakata et al. [29] that has similar form. The ik terms in the exponents in Equation A3 and Equation 9 in Nakata et al. [29] are all positive. Hence, the wavefield is asymmetric.

Appendix A.3. The effect of reflection coefficients
We vary the reflection coefficient on the top boundary (R 1 ) from 0.01, 0.5, to 0.99. The effect is not obvious in the mathematical notation but observable in the simulated wavefields ( Figure A3). When R 1 gets larger, for channels above the virtual source, the acausal waves are enhanced; whereas for channels below the virtual source, the causal waves are enhanced.
This phenomenon of a larger R 1 is the opposite of the effect of larger S 1 . That is, based on A.2, if |S 1 | increases relative to |S 2 |, we expect the causal waves being enhanced for channel above the virtual source while the acausal waves being enhanced for channel below the virtual source. Hence, the relative amplitudes between the causal and acausal axes can be affected by both the relative source amplitude and the reflection coefficients. However, the influence of the reflection coefficient on the symmetry is subtle. We set R 1 = R 2 = 0.9 in the simulations in

Appendix B. Deconvolved wavefields at the lower part below 200 m
We extract coherent waves from the formations. But these waves are only present during the times when the vibroseis shots were close to the DASV well. Figure B1 shows three sets of waves we extract between depth range 165-300 m. This deconvolved wavefield is calculated using 30 minutes time window, 50% overlap, and stacked over 3 hours.
The first two signals ( Figure B1a-B1b) travel downward with the apparent velocities of 2100 m/s (green dashed lines) and 1100 m/s (pink dashed line). The V p /V s ratio is 1.91. This is consistent for shallow formations in Brady that consists of volcanic sediments, limestone, lacustrine sediments, and geothermal features such as carbonate tufa [47]. The measured velocities are close to previously estimated local velocities (V p = 2300 m/s and V s = 1200 m/s; Parker et al. 39;Matzel et al. 46). The slower apparent velocities might be due to incident angles. Potentially, we could estimate the time-lapse changes by measuring the relative velocity changes of these waves. However, the poor coupling condition prevents us from getting more scattering energy.
The third signal ( Figure B1c) has apparent velocity of 1400 m/s (the yellow dash line) and propagates upward. It is weaker than the first two signals. The source of this signal can be a reflection from nearby faults or bedding planes [26,47]. This the most likely case that we can have an up-going wave here. However, we cannot identify the reflection point due to limited amount of good-quality data.