Modeling Subsurface Explosions Recorded on a Distributed Fiber Optic Sensor

Fiber optic distributed acoustic sensors (DAS) are becoming a widely used tool for seismic sensing. Here we examine recordings of two subsurface chemical explosions, DAG‐1 and DAG‐3, each of which was about one metric ton (TNT equivalent), that were recorded from a helical fiber installed in two boreholes 80 m away from the source location. Several clear phases including the initial P wave, a weak S wave, and a surface reflected P wave are observed on the helical DAS data. We estimate a velocity model using arrival times measured from the fiber. The DAS waveform data were compared with colocated accelerometers at specific depths in both frequency and time domains. The spectra of the DAS data matched spectra estimated from the accelerometer records. Comparisons of observed waveform shape between the accelerometer records and the fiber measurements (strain‐rate) show reasonable agreement except for the data near the event depth. The DAS data and the accelerometer agreed in relative amplitudes but we had difficulties in matching absolute amplitudes, possibly due to errors in metadata. Synthetic strain‐rate waveforms were calculated using a 2D wavenumber algorithm and matched the waveform shape and relative amplitudes. In general, DAS is effective at recording strong ground motions at high spatial density. Comparison of the synthetic seismograms with observed data indicate that the waveforms are not consistent with a pure isotropic explosion source and that the observed S waves originate from very near the source region.

The borehole DAS measurements provided an unprecedented snapshot of the seismic wavefield from the explosion and record both P and S waves. Previous analysis of the wavefield has revealed new observations about the nature of spall created by the wavefield (Abbott, 2019). The goal of this paper is to conduct preliminary modeling and understand the strengths and limitations of the sensors and how the measurements can be used to constrain more sophisticated analysis focused on the S wave, as well as to provide a basis for future experiments and deployments. An understanding of the fiber response, and the ability to model the observed signals, will greatly enhance the utility of this data set as well as future datasets and provide constraints and insight into the S wave generation.
One challenge is understanding the precise response of these fiber sensors with respect to frequency and amplitude of the seismic wavefield. The response of standard seismic sensors and associated digitizers can be characterized in a variety of ways but is typically accurate within 5% of the actual response within a specific bandwidth for calibrated observatory grade instruments (Ringler et al., 2012). That is, after application of the instrument correction, the measured ground motion should be very close to the actual ground motion within a specific bandwidth. The underlying technology and basic measurement of fiber sensors differs significantly from seismic sensors currently commonly used in monitoring and uncertainty exists in the transfer function needed to convert the data from fiber sensors to ground units. Previous quantitative studies comparing DAS and geophones in boreholes suggest that DAS provides accurate estimates of seismic travel times and amplitude for weak motion when compared with geophones (Booth et al., 2020;Correa, 2018;Willis et al., 2016). Here, we evaluate the DAS response for strong motion by comparison with synthetics and co-located accelerometers and geophones.
Fiber optic sensors are often referred to as fiber optic acoustic sensors (FOAS) or distributed acoustic sensing (DAS). As the sensors respond to both scalar acoustic and vector seismic waves, referring to the sensors as "acoustic" is somewhat misleading, but, following current convention, we will refer to the sensors as DAS in this paper. In the following section, we provide an overview of the basics of current DAS Rayleigh-scattering technology, but we note that distributed sensor technology is an area of active and rapidly progressing research and we focus on topics directly relevant to the data set rather than attempt a comprehensive review of all related technologies.
Optical fibers possess a long history as acoustic and hydroacoustic sensors (Budiansky et al., 1979;Hughes & Jarzynski, 1980;Kirkendall & Dandridge, 2004), strain meters (DeWolf et al., 2015), and point fiber Bragg sensors (Laudati et al., 2007). Fiber-based distributed acoustic sensors differ in that they measure the signal at all points along the fiber using the fiber itself as a sensor (Hartog, 2000). This is enabled by taking advantage of inherent (or engineered) heterogeneities in the fiber that cause Rayleigh backscattering of an input signal (Hartog, 2000). Other types of scattering such as Raman or Brillouin also occur in optical fibers and are useful for distributed temperature or low-frequency strain measurements (Miah & Potter, 2017;Soga & Luo, 2018), but Rayleigh scattering is preferred for seismic sensing as the frequency of the back-scattered 10.1029/2021JB022690 3 of 22 optical pulse remains unaltered (unlike Brillouin or Raman scattering), and it is possible to measure at the frequencies necessary to measure seismic signals.
An optical fiber consists of a core surrounded by cladding, both of which are composed of silica glass but with slightly different refractive indices. Propagation may be singlemode or multimode, depending on the core diameter and frequency of the light (Miah & Potter, 2017). Single mode usually shows the least attenuation (roughly 0.2 db/km at 1,550 nm wavelength) and allows longer ranges than multimode fibers. Each fiber is coated, typically with polymer, for mechanical protection. The fiber, or fibers, are then bundled within a cable. The design of the cable depends on the intended use but the cable generally composed of polymer reinforced with steel and possibly with gel filling for strain relief during installation. It is possible to wind the fiber in a helix around the cable to vary the angular response (Kuvshinov, 2016) and these helical cables are generally thicker than standard cables.
Rayleigh distributed seismic sensing measures the returned phase of the back-scattered light from a laser pulses sent along the fiber (Masoudi & Newson, 2016). The scatterers are small variations in refractive index in the fiber that occur as part of the manufacturing process or by a deliberate enhancement. When a seismic wavefield crosses a fiber, the induced strain changes both the refractive index and the relative position of the scattering centers, which in turn alters the intensity and phase of the back-scattered light. Unfortunately, the change in intensity is not proportional to the strain amplitude, which makes precise reconstruction of the strain using only intensity impossible. However, the relative phase change is linearly related to the strain amplitude (Hartog, 2000). Hence, precise seismological measurements require an interrogator capable of measuring differential phase of the signals. The differential phase measurement is made using interferometric techniques and several implementations are possible (Hartog, 2000). Depending on the implementation, the interrogator may be sensitive to strain or strain-rate in the fiber. The measurement of the returned signal is made by a photodetector and the signal is typically converted to in-phase and quadrature components, demodulated, converted to phase, and then unwrapped (Nishiguchi, 2016).
The strain in the fiber is induced by two effects: direct strain coupling of the fiber to the surrounding solid media and perpendicular stress on the fiber which creates additional fiber changes due to Poisson-induced elongation on the fiber, fiber coating, and cable (Miah & Potter, 2017). In water, for example, where shear coupling is minimal, the optical fiber responds to pressure changes and effectively act as a hydrophone (Hughes & Jarzynski, 1980). In solid media with a well-coupled cable, direct transfer of strain from the surrounding media to the cable is the predominant mechanism (Reinsch et al., 2017). As noted above, photo-elastic effects change the refractive index of the fiber as well (Kuvshinov, 2016) attributes 80% of the phase change for a well-coupled fiber to a change in length with the remainder due to photo-elastic effects. A key distinction is that pressure-induced fiber elongation should be independent of fiber orientation while strain coupling will vary with fiber orientation. Temperature variations will also cause variations at low frequencies (e.g., <0.1 Hz) but may be mitigated by using two fibers with different coefficients of expansion or measured with a distributed temperature sensor.
The transfer efficiency of strain coupling for most seismic studies with strongly embedded cables appears to be close to 100% (Reinsch et al., 2017) at the relatively low frequencies encountered in seismic studies. Fiber-based laser strain-meters match adjacent vacuum strain-meters very well (within 5%) up to periods of several days (DeWolf et al., 2015) and that provides confidence that fiber-based sensors can measure ground motion with high fidelity. In other settings such as cables inside casing in boreholes or in conduits, as is common for telecom fibers used as sensors, the transfer efficiency is not well quantified and may vary spatially. Changes in material properties between the surrounding media and the cable will also affect the amplitude of the transmitted shear strain. These can be modeled in a variety of ways including semi-analytical methods (Hughes & Jarzynski, 1980) and wavenumber methods (Wang et al., 2018).
Noise in DAS arises from several effects, both in the fiber and in the interrogator. One source is destructive interference (fading) along the fiber, which reduces the amplitude of the already weak signal (Gabai & Eyal, 2016). To improve signal-to-noise, the back-scattered returns are averaged along a segment of fiber, referred to as the gauge length, which typically ranges in length from 1 to 10 m . This length is adjustable, and longer gauge lengths permit longer sensing distances at the expense of spatial resolution. The gauge length usually exceeds the channel spacing, which is the spatial distance between adjacent sampling points, each of which yields a separate time series. As the channel spacing is usually less than, or equal to, the gauge length, measurements at adjacent channels are not completely independent. Specially engineered fibers can be fabricated to enhance signal response, which allows greater sensitivity and shorter gauge lengths. For the signals in this study, which are generated by chemical explosives detonated 80 m from the fiber, low signal amplitudes are not a significant problem but the shorter gauge lengths aids in spatial and dynamic resolution.
The electronics inside the interrogator box are also a source of noise (Kirkendall & Dandridge, 2004). Laser drift, for example, may cause amplitude variations for a specific pulse at all channels, simultaneously. For some interrogator designs, vibration of the interrogation box will affect the reference fiber coil and contaminate the desired signal. Postprocessing steps such as common mode filtering can be implemented to alleviate these issues. Another potential source of noise is due to errors in phase unwrapping. Typically, if implemented, this is performed in realtime with fast 1D unwrapping algorithms. Phase unwrapping errors would differ on each channel and tend to reduce the measured signal amplitude and underestimate the actual signal. It is possible that phase unwrapping errors might cause abrupt discontinuities between channels but the exact expression depends on the unwrapping algorithm. Unwrapping algorithms implemented as part of a high-speed real-time data acquisition in an interrogator may differ from standard unwrapping techniques and may include filtering as part of the implementation. As far as we know, the only way to definitely identify phase errors would be to have access to the original output of the interrogator photodetector prior to unwrapping. This would permit testing of alternate unwrapping schemes or error detection (Ghiglia & Pritt, 1998).
An operational source of error is that the exact spatial location of each channel is uncertain. While the distance that the signal travels along the fiber is known, the correspondence of the fiber distance with the exact spatial locations can be difficult to ascertain. For surface fibers, tap tests by striking the ground near the fiber can be conducted which allows identification of a specific channel location (within ∼1 m). Borehole fibers are more difficult and require colocated borehole sensors (or known subsurface sources) combined with tap tests at the surface. Channel locations are then interpolated between the known points.
In terms of response, the basic concern is the amplitude and frequency response. In principle, amplitude is straightforward. The change in phase of the backscattered laser pulse is proportional to the change in fiber length, and therefore proportional to the strain amplitude of the seismic wave. For strain rate sensors, the phase is proportional to the change in strain with respect to time. This proportionality factor depends on wavelength and the photo elastic refraction index for silica glass, both of which are known. The exact numerical factor depends on interrogator and fiber design but an example of the nominal response for a generic interrogator and standard fiber (non-engineered) is provided in SEAFOM (2018). This sensitivity factor will differ, perhaps significantly, for engineered fibers that are designed to enhance the returned signal.
Calibration of the nominal fiber response is typically performed either by inducing a known strain (SEA-FOM, 2018) or by colocation of a seismic sensor(s) (Lindsey et al., 2020). Neither approach is optimal. The induced strains generated on a lab bench differ in temporal and spatial characteristics from those produced by a transient seismic wave. Co-located seismic sensors require conversion between strain/strain-rate over the gauge length at a specific orientation and velocity/acceleration at a point. As a result, the transfer function between DAS and point sensors is complicated although the amplitude/frequency response of the DAS itself to ground motion is reasonably simple (Jousset et al., 2018). We are unaware of accurate numerical frequency/amplitude representations such as pole-zero representations for specific interrogators, as is standard for seismic sensors. As mentioned above, simple scalar conversion constants based on fiber parameters and interrogator settings are typically available, which convert the digital counts to strain or strain-rate but are independent of frequency.
The high-frequency response depends mainly on the sample rate and the effective gauge length. For the DAG series, where high amplitude signals were expected, high sample rates were tested to try to ensure high fidelity of the signal. In turn, this limited the possible length of the interrogated fiber, as each measurement requires sufficient time for a laser pulse to travel the length of the fiber (and back). In any case, the original sample rates (10's of kHz) are much higher than standard seismic sensors and the data are typically down sampled to reduce data volume. While details are proprietary, DAS interrogators may not include anti-aliasing filter prior to digitization, at least not in the exact fashion as the digitizers used in seismic sensors. 10.1029/2021JB022690

of 22
More importantly, the gauge length, which is the distance over which the strain field is measured, acts as a spatial filter and signals with apparent wavelengths at or less than the gauge length will be aliased (Figure 1). Even apparent wavelengths up to five times the gauge length will be decreased in amplitude Lomnitz, 1997). As the apparent wavelength depends on both wave speed and angle, this means that the response varies with the direction of propagation of the seismic wave. For DAG with an effective fiber gauge length of 1 m, wavelengths at this length occur at ∼1,000 Hz for P and 600 Hz for S in the slowest materials near the surface, with essentially no reduction in amplitude expected at frequencies lower than 200 and 120 Hz respectively. As the fiber is wound in a helix, the effective gauge length with respect to the apparent wavelength of the seismic waves is shorter (Ning & Sava, 2018) and hence the frequency resolution may be slightly higher. At the other end of the frequency bandwidth, very low frequencies with periods of 1,000's of seconds (Becker & Coleman, 2019) can be measured but the sensitivity appears to depend on the implementation, system noise, and the application of various corrections to account for factors such as temperature variation.
Another factor is the directional nature of the measurement, which differs from the point sensors common in seismology. Fiber sensors measure a specific component, or combination of components of the average strain, a tensor, over a specific length. The amplitude and frequency response of a linear fiber with a 10 m gauge length is similar to that of a 10 m strain-meter (Benioff, 1935, Figure 1). For P waves, maximum amplitude sensitivity occurs when the phase is traveling parallel to the fiber and decreases by a factor of cos 2 (θ) as the angle of incidence (θ) increases (Mateeva et al., 2014). For an S wave, the response with azimuth is sin (2θ) and therefore possesses four maxima. The pronounced directional response of a linear fiber can be alleviated by winding the fiber in a helix around the cable. This increases sensitivity to P seismic waves that are not traveling parallel to the fiber (Kuvshinov, 2016) and field tests confirm this capability (Hornman, 2017;Yavuz et al., 2019). However, a helical winding reduces the sensitivity to S waves (Baird, 2020).
As described above, the fundamental signal is the strain-rate along the fiber as measured over a set of spatially distributed sections. Our goal is to model these signals and gain insight into the strengths and weaknesses of quantitative analysis of fiber seismic data at these distances and amplitudes. In addition, as described above, one of the challenges is the lack of a well-defined response for fiber sensors. Modeling may provide some insight into defining the response for the sensor and ideally, it may be possible to construct a response from first principles. We examine data from DAG-1 and DAG-3 to understand how well we can model the data using simple assumptions.

Data
The data set is derived from the Dry Alluvium Geology (DAG) series of the Source Physics Experiment, which consisted of four subsurface chemical explosions at depths ranging from 385 (DAG-1) to 52 (DAG-4) meters in a borehole at the Nevada National Security Site (Bonner, 2018). These chemical explosion events, which varied in size between one and 50 metric tons (TNT equivalent), were recorded by fiber cables installed in boreholes located 80 m from the source borehole ( Figure 2). The fiber cable was a custom helically wound engineered fiber designed to be sensitive to multiple components of the seismic wavefield and with higher sensitivity. The boreholes were also equipped with accelerometers at multiple depths and a single geophone at the bottom of the SW80 borehole. The boreholes are emplaced in the alluvial fill of Yucca Flat, a fault-bounded basin that is filled with eroded debris from the surrounding mountains. In general, the fill consists of mixtures of silt, sand, and gravel deposited in coalescing alluvial fans. Yucca Flat has been extensively drilled, logged, and the target of multiple geophysical surveys and the subsurface is well characterized both from a geological (Wagoner & Prothro, 2020) and seismological perspective Pitarka & Mellors, 2021).
The data set ( Figure 2) used in this study consists of data from DAG-1 and DAG-3 which were recorded on fiber installed in two boreholes situated 80 m away from the source emplacement borehole (U2ez). The two monitoring boreholes are referred to, using a sophisticated naming scheme, as E80 (80-m east of the explosion) and SW80 (80 m southwest). SW80 has a total depth of 450 m and E80 is 385 m deep. All boreholes are vertical. Both DAG-1 and DAG-3 consisted of one metric ton (TNT equivalent) of chemical explosive with similar composition, canisters, emplacement, and stemming. The source depth for DAG-1 and DAG-3 was 385 and 150 m, respectively. The series included two other events, DAG-2 (50 metric tons at 300 m) and DAG-4 (10 metric tons at 52 m), which were also recorded by the same fiber sensors. However, as DAG-2 and DAG-4 were larger, nonlinear interactions may play a significant role in the nearfield wavefield, which increases the complexity of the analysis. Therefore, in this paper we focus on DAG-1 and DAG-3, as we believe that the fiber was primarily in the elastic zone (Perret & Bass, 1975) for these events, and therefore use of a linear elastic algorithm to model the wavefield at the fiber is reasonable as a first-order approximation.
The data were recorded with a channel spacing of 1.02 m and gauge length of 2 m using an engineered fiber wound in a nominal 30° helix (fabrication angle measured with respect to the cable axis) around a cable with a diameter of roughly 5 cm. The fabrication angle is equal to 90, helical angle. The effective gauge length is roughly 1.0 m, as each meter of borehole accommodates 2.0 m of fiber. The term "engineered" refers to an optical fiber that has been custom fabricated to yield high sensitivity (up to 100× times that of standard telecom fiber) and precision (Naldrett et al., 2020). The data are strain-rate as a function of time and distance along the fiber and the delivered data is sampled at 2,000 Hz having been decimated from an original sample rate of 50,000 and 100,000 Hz for E80 and SW80, respectively. As this type of close-in deployment had not been attempted, previously, the setup and recording parameters were very much in experimental mode, especially for DAG-1, and several iterations of processed data products were produced. This led to complexities and uncertainties in the metadata.
The Silixa system is an interferometer, which does not keep track of integer multiples of full cycle strain along gauge lengths between subsequent measurements. This means that large strain rate signals, greater than 2 pi radians over a gauge length, cannot be recovered with phase unwrapping. There are two ways to increase system 10.1029/2021JB022690 7 of 22 dynamic range to ensure that the system does not saturate when measuring large strain rate signals and both were employed for DAG: decrease the gauge length and increase the sampling frequency. Decreasing the gauge length simply reduces the strain rate base length and makes each sample less sensitive to strain changes. Increasing the sampling frequency samples large strain rate signals at points more closely separated in time, and so the phase change from one time step to the next is smaller and will more likely fit within the phase limitations imposed by the interferometer. As it was unclear prior to the experiment which sample rate would be optimal, both 50,000 and 100,000 Hz sample rates were tested. Figure 2 shows DAG-1 as recorded on both wells and DAG-3 on E80 (only one of every 10 channels is shown). The data quality is good for DAG-1 but data from DAG-3 is noisier (with more spikes and missing traces). The cause is unclear but may be due to poorer coupling, perhaps due to permanent subsurface disruption caused by the large DAG-2 event. The direct P is clearly visible along with a weaker S arrival. A reflected PP phase is also evident. In general, the DAG-1 traces appear to be higher frequency than the DAG-3 counterparts. The P phase near the shot depth also appears weaker and less impulsive on DAG-3 than on DAG-1. The PP phase also appears much lower frequency than the direct P.
Along with the fiber data, a set of borehole accelerometers was deployed in each borehole for the DAG events (Cashion, 2018;Steedman, 2019). Accelerometer packages were installed at depths corresponding to the event depths (50, 150, 300, and 385 m) as well as some supplementary levels for specific events (Figure 2). The SW80 well also had an accelerometer and a geophone at a depth of 450 m. Each accelerometer installation was composed of three sets of three-component accelerometers, yielding nine total measurements. The triple redundancy was installed to compensate for possible failures. Each individual accelerometer package included two orthogonal horizontal components and one vertical, but the packages were oriented 120° apart on the horizontal plane. The data was rotated (2D) into an L (longitudinal or vertical), R (radial), and T (tangential) components after acquisition so that the L component was parallel to the borehole. 3D rotated data, which was rotated in a vertical plane as well to align the radial with source hypocenter, was also available but we use the 2D rotated data in this study, as it should be more consistent with the borehole fiber geometry. We examined the 3D rotated data as well, but the results were similar.
We compare the DAS and the accelerometer data for consistency. Initially, we looked at the raw data (DAS data in nominal strain rate and m/s/s for the accelerometers) (Figure 3). Several observations can be made. The waveforms at E80 and SW80 are similar although there is a systematic difference by a factor of two between the two boreholes (which is corrected in Figure 3 and which we believe is due to an artifact in the metadata). Near the event level for DAG-1 (385), the waveforms recorded at the two boreholes differs. The accelerometer data, on the other hand, is similar in amplitude between the two boreholes with less obvious differences in waveform at the event depth. To provide a common basis for quantitative comparison, the data were resampled to a common sample rate (1,000 Hz) and then filtered between 1 and 50 Hz. We did not include the single record from the velocity sensor in the SW80 well at the 450 level in the analysis but did confirm that it was essentially identical (after differentiation) to the accelerometer waveform at that point. The DAS data were extracted directly from the SEGY files. The comparison was conducted using correlation to examine waveform shapes and the relative scaling of the raw data using measured peak amplitudes of the initial arrival. Four sets of correlations and amplitudes were estimated: (a) cross-correlation of the DAS data between the two wells (E80 and SW80) at identical depths for both events (e.g., DAG-1 E80 channel at 100 m was cross-correlated with DAG-1 SW80 channel 100 m); (b) cross correlation between DAS and accelerometer (L) waveforms at the same depth; (c) cross-correlation of the accelerometer data between the two wells for each event (e.g., E80 sensor 1 L at 150 m with SW80 sensors 1 L at 150 m); and finally, cross-correlation of the three sets of accelerometer data at each level for each well and event (e.g., DAG-1 E80 The DAS/DAS correlation was estimated using channels spaced 10 m apart from the surface to a depth of 380 m (bottom of the E80 well). Correlation (absolute value) was high, with a mean of 0.88 (±0.09) for DAG-1 and 84% (±0.14) for DAG-3. Correlation was lowest near the event depth and on channels where the S wave, which is more pronounced on SW80, was strongest. The amplitudes between the two wells differ systematically for DAG-1, with SW80 showing higher amplitudes by an overall factor of ∼2. This is not clearly observed for DAG-3 or in the accelerometer data and we suspect that it is due to an error in postprocessing but were unable to confirm. The analysis also showed a difference in polarity between the two wells for DAG-3, which is also an artifact from the data acquisition.
We compared the DAS data with accelerometer data (L) in the same well and at the same depth. Correlation was as high as 0.99 for specific pairs but decreased for sensors at the event depth. Also, some specific accelerometer packages showed considerable variation even within sensors in the same package. The overall mean was 0.88 (±0.12) for DAG-1 and 0.86 (±0.11) for DAG-3. For DAG-1, the polarity was reversed between the DAS and the acceleration at 450 m, below the depth of the event. This is expected, as DAS data, which measures strain rate, is insensitive to changes in direction of wave motion but the accelerometers do distinguish between upward and downward motion. This switch in polarities occurs for DAG-3 as well.
We note that the accelerometer data set was extensively examined by Steedman (2019) and our analysis is consistent with those results, although the analysis presented here is not as comprehensive. In at least one case, (SW80 450), the amplitude varies by a factor of almost 10, which suggests an issue with the instrumentation. We use the correlations as a guide for reliability of the data of the sensor packages in the following analysis. Accelerometer data from a sensor package with inconsistent data were not used in the comparison with the DAS time series.
To evaluate the data further, we estimated the spectra for each DAS channel and then smoothed by averaging every 10 channels. A multitaper algorithm (Thomson, 1982) as implemented by (Krischer, 2016;Prieto et al., 2009) was used. The results are shown in Figure 4. We see a distinct shift toward lower frequencies from DAG-1 to DAG-3. Peak frequencies for DAG-1 lie between 50 and 100 Hz, while DAG-3 are almost uniformly less than 50 Hz. A check of the acceleration records yielded similar spectral results (shown in red on Figure 4) as the fiber, with the exception of the DAG-3 accelerometer at 380 m for SW80. Previous work (Hornman, 2017) also reported good correspondence between the spectra measured by helically wound fiber and geophones. As strain-rate on a linear fiber is related to acceleration by a scalar factor (apparent phase velocity), consistency in spectra between the two measurements is expected.

Travel-Times and Velocity Model
The first step in analysis of the DAS data was to measure phase travel times and estimate a velocity model based on the fiber data. The velocity model is useful for assessment of reliability of fiber for travel times, comparison with previous velocity models and for computing synthetics. Abbott et al. (2019) used the DAG data set to invert P arrival time using data from all four DAG events. One challenge is that the velocity model appears to vary over the experiment sequence, as the explosions disrupt the subsurface sufficiently to create variations in the seismic velocities. DAG-3, which was equivalent to 50 tons of TNT, causes a clear decrease in seismic velocities centered 10.1029/2021JB022690 10 of 22 at the depth of the event . We decided to redo the velocity model analysis using a simple approach to ensure clear understanding of the data set and consistency with the modeling.
As noted above, the mapping of channel locations to depth is a key element in analysis of borehole fiber and we used the surface and a colocated accelerometer to match channel locations to depth. While a wide variety of approaches are possible to estimate velocity from a vertical profile (Lellouch, Yuan, Spica, et al., 2019), the technique we used to estimate the velocity model was simple. P phase arrival times were picked using an automatic . Normalized spectra of DAG-1 and DAG-3 DAS data for wells SW80 as estimated from DAS (black) and selected accelerometers (red). Each DAS spectra is an average of the spectra from 10 channels. Yellow stars denote the position of DAG-1 and DAG-3. The location of DAG-2, a large (50 metric ton TNT equivalent) event between DAG-1 and DAG-3 is also shown.
picker and the results reviewed manually. We then conducted a grid search to match the observed P travel times. The velocity model was parameterized using two variables: surface velocity and a constant velocity gradient. A WKBJ algorithm used to calculate travel times from the source to the fiber channel locations (Stockwell, 2008). The parametrization assumes that the velocity increases with depth, which is reasonable given the Yucca Flat velocity structure (Toney et al., 2019;Wagoner & Prothro, 2020). The advantage of a grid search is that it provides a robust indication of the error. For DAG-1, estimates of the velocity model were independently made for each well (E80 and SW80). The difference between the two wells was within the variance of the combined results. Therefore, we use a common velocity model constrained by data from both wells. The predicted travel times match the observed ones well for DAG-1 and the results were consistent with previous studies and a sonic log of the U2EZ borehole ( Figure 5).
As a check, we tested an alternate algorithm that used a staged Monte Carlo approach. This algorithm included a specific number of depth points with randomly varying velocities. Depths between the specific points were interpolated using a cubic spline. The first stage defined the vertical profile with three points, which then increased to five and seven points. Each stage uses the best fit from the preceding stage as a starting point. The initial model was a constant velocity (set at 1,750 m/s). This algorithm, unlike the preceding gradient grid search, allowed for the possibility of low velocity zones. Results for DAG-1 were similar to the gradient grid search although the results demonstrated that the velocities at the top and bottom of the data set are poorly constrained by these algorithms. For DAG-3, both the gradient algorithm and the random algorithm showed indications of a decrease in velocity with respect to DAG-1. DAG-3 data also showed more variation in the estimated velocities between the E80 and the SW80 borehole ( Figure 5), which matches conclusions reached by (Abbott, 2019).
To match the S waves, we estimated the Vp/Vs ratio using a Wadati diagram with a known origin time (Kisslinger & Engdahl, 1973). For DAG-1, although some differences were notes between SW80 and E80, the analysis yielded a value of 1.75, which roughly matches Vp/Vs for the same location obtained by Toney et al., 2019 based on active source surface wave dispersion data. For DAG-3, the results from the Wadati-based analysis was poorer, possibly indicating changes with depth and yielded a Vs/Vp of 1.80. We suspect that the disruption caused by DAG-1 and DAG-2 contributes to the observed increase in Vs/Vp. However, when tested with synthetics, the velocity models underestimated the expected PP travel times. Tests using the PP time as an additional constraint resulted in significantly poorer fit to the initial P arrival times. We suspect that this misfit is the result of several factors: (a) previous studies (Wagoner & Prothro, 2020) indicate lower velocities near (within ∼30 m) of the surface; (b) the accelerations at the surface cause spall, in which the upper surface layers physically separate from the underlying subsurface at the initial arrival; and (c) systematic delay in picking the PP arrival time due to high background noise from coda. Spall in the upper layer leads to a decrease in velocity and higher attenuation, which would explain both the delayed arrival and the decrease in frequency content. We suspect that spall on the upper near-surface layer may be the primary cause. The DAG fiber data set deserves a more comprehensive study to understand the variations in the velocity model induced by the subsurface events, but for the purposes of this study, the fit to the observed travel times is adequate to explore modeling of the fiber data, which is our primary objective.

Fiber and Accelerometer Waveforms
We compared the fiber waveforms with the colocated accelerometers in SW80 and E80. In theory, strain-rate measured on a linear DAS fiber can be directly converted into acceleration. For a plane wave in a homogenous media where the seismic wavelength is much longer than the gauge length:  is strain-rate, and E p is apparent phase slowness (Lior et al., 2021). In the case of DAG, the velocities are not homogenous, and the wavefront is not planar, but tests using acceleration and strain-rate synthetics indicated that the above relationship is a reasonable approximation for DAG. The major problems are due to estimation of the apparent velocity, which varies with time along the waveform with each phase (P, S, and PP), and effects of the helical fiber. Initially, we compare waveform shape and apply an empirical scale factor to match amplitude for ease of comparison. To avoid comparing with accelerometer data of uncertain quality, we only used data from accelerometer packages where all three sets of accelerometer sensors agreed with each other. We use the L component of the accelerometer for comparison, as it displayed the highest correlation with the DAS data at these depths. Figure 6 shows a good match in waveform shape between the DAS and the L accelerometer at depths away from the source depth for both SW80 and E80 for DAG-1 and DAG-3. The DAG-1 source was at a depth 385 m, so all measurements are located at distances of 85 and 235 m above the explosion. Given the geometry, and the relatively strong velocity gradient, it is not too surprising that the signal from the vertical fiber matches the L accelerometer component most closely. The DAG-3 accelerometer data (L) at the 300 or 385 levels also matches the DAS Figure 6. Comparison of Dry Alluvium Geology (DAS) fiber data amplitudes (black) with accelerometer (red) at selected receiver locations and components not within 80 m of the source depth. The match between the DAS and the accelerometers is good at these locations. The DAS data have been multiplied by a factor listed in Table 1. data reasonably well. The visual comparison matches the cross-correlation estimates. The similarity between accelerometers and DAS was higher for DAG-1 than DAG-3. We suspect that the large DAG-2 event may have affected coupling of both the fiber and possibly the accelerometer packages. Figure 7 shows the seismograms at source level (385 m for DAG-1 and 150 m for DAG-3). Here the match between accelerometer data and DAS data is not as clear and we show all three accelerometer components for clarity. For DAG-1, the DAS data does not match any accelerometer component particularly well and the amplitudes at 0.1 s are as large as the initial amplitudes, unlike the accelerometer data. For DAG-3, the highest amplitudes occur at approximately the same arrival time as the accelerometer data but the waveforms differ. Some differences are expected. In theory, the helical wound DAS should possess some sensitivity to the pure broadside radial arrival but the response should be reduced in amplitude. Similarly, the helical fiber should be sensitive to the L and T components as well, but with different response due to the asymmetry in fiber geometry in the L direction as compared with the T and R directions. Finally, the fiber measures strain-rate, which is proportional to acceleration but may differ in polarity depending on the apparent phase velocity. In general, the variations in amplitudes observed for DAG-1 on the DAS are surprising.
One possibility is that the high amplitudes in the radial direction simply exceed the capability of the interrogator to respond. If the signal change is fast enough (more than 2π radians per sample), the wrapped phase will be aliased and in general, the aliased signal should be smaller in amplitude that then actual signal (SEAFOM, 2018). Unlike standard seismic instruments, the fiber signal will not clip but will typically display a ragged signal with lower amplitude. Phase saturation errors should also cause differences between adjacent DAS traces even if the Figure 7. Comparison of Dry Alluvium Geology (DAS) fiber data amplitudes (black) with accelerometer (red) at depths near the source depth. The DAS data has been multiplied by a factor scaled to the expected apparent velocity with respect to the vertical fiber and by an empirical factor (Table 1).
channel separation is less than the gauge length, and examination of adjacent traces shows some differences. These differences appear more distinct on E80 than SW80.
Next, we examine the amplitudes of the DAS data. Due to approximations required for helical fiber, as well as differences caused by uncertainties in the estimates of the acoustic impedance of the cable, grout, and surrounding rock, we chose to estimate empirical scaling factors for the DAS data based on the accelerometer data, which we assume is correct. These empirical factors were applied in Figures 6 and 7.
This empirical scale factors (Table 1) were estimated in two ways: (a) a simple ratio between the maximum absolute amplitudes of the vertical acceleration and the DAS and (b) the best-fit scale factor to match the two waveform amplitudes over a time window centered on the P arrival. In general, the two different methods provided similar results except where the two waveforms (DAS and accelerometer) showed significant differences in shape. A complication is that possible errors in metadata during data acquisition adds further confusion, as we believe there is a systematic factor of two between E80 and SW80 for the DAG-1 data

Empirical and Theoretical Scaling Between DAS and Accelerometers
In theory, these empirical values can be used to test various models of the helical fiber response and estimated parameters (Baird, 2020;Kuvshinov, 2016;Ning & Sava, 2018). Direct comparison with the accelerometer data requires projection of the strain tensor components measured by the fiber on the vector components measured by the accelerometer and the appropriate apparent phase velocities, which will differ in each direction. For a linear fiber, this process is reasonably straightforward and conversion of strain-rate amplitude to acceleration amplitudes for a linear fiber requires multiplication by the apparent phase velocity (Wang et al., 2018). For helical fiber, the transformation between DAS and a point sensor is more complex, especially when the seismogram contains different phases (P, S, and PP), each with differing phase velocities.
We attempted to directly compare the DAS strain-rate and the accelerometer waveform amplitudes. We focus on the DAG-1 seismograms shown in Figure 6 that match closely in waveform shape and at distance where the wavefront is traveling at a relatively small incident angle to the fiber axis. At receiver locations away from the source depth, we assume that the fiber response can be approximately by a linear fiber in the Z direction and the correction will be proportional to the apparent phase velocity and the fiber sensitivity. The vertical phase velocity at each depth was estimated using three approaches: (a) observed travel times, (b) estimates from the velocity model, and (c) estimates from the ratios of acceleration and strain-rate derived from synthetic waveforms. Phase velocities estimated from the observed travel times and velocity model are highly sensitive to minor differences in the phase arrivals, especially at the event depth. Estimates from the synthetics were slower at the event depth due to finite frequency content of the synthetics. Unfortunately, the expected theoretical corrections did not match the observed empirical correction used to generate Figures 6 and 7 and listed in Table 1. This may be due to incomplete modeling of the response of the helical fiber, errors in the phase velocity estimates, and inaccuracies in the existing metadata and fiber parameters especially with respect to the fiber sensitivity factor. We were unable to resolve the differences despite several iterations in reprocessing of the original data set.

Synthetics
We attempt to model the DAS signal directly. The goals are to test our understanding of the fiber response and to constrain the source. We use two approaches to modeling: a 1D wavenumber code (Herrmann, 2013) and a 3D finite difference code (SW4) (Sjögreen & Petersson, 2012). The advantage of the wavenumber code is that it is fast and can generate high frequencies. The disadvantage is that only displacement, velocity, or acceleration time series can be created. It is possible to difference velocity or acceleration time series to match single components of the strain tensor measured by DAS (Wang et al., 2018), but generating multiple components is more complicated, al- Note. The DAS traces shown in Figures 6 and 7 were multiplied by the empirical scale factor derived from the ratio between the DAS and the vertical acceleration records. Asterisks indicate event level and the L or the R denote accelerometer channel.

Table 1
Measured Parameters From Dry Alluvium Geology (DAS) and Accelerometers for Selected Traces though possible. The finite difference code will generate the strain tensor as well as velocities and displacement, and can handle 3D variations in velocity and anisotropy, but high frequencies are computationally intensive. We restrict the synthetics to frequencies below 100 Hz to avoid complications due to the gauge length. Assuming a minimum seismic velocity of 800 m/s (S waves near the surface) and a 100 Hz signal, the wavelength is 8 m, which is much longer than the effective gauge length. The spectra shown in Figure 2 also indicate that the signals are mostly below 100 Hz.
To approximate the strain in the longitudinal direction using the wavenumber code, we calculate wavenumber synthetics on a 2 m vertical spacing and difference the Z component to estimate the strain-rate waveform. We chose 2 m as it was smaller than the shortest expected wavelength but longer than the effective gauge length. The velocity model derived from the travel time grid search was used and Vp/Vs from the estimates based on the Wadati diagrams (1.75 and 1.80, respectively for DAG-1 and DAG-3). We also created acceleration synthetics to compare directly with the accelerometers. As discussed earlier, the ratio between the strain-rate synthetics and the acceleration synthetics can be used to estimate the effective phase velocity at each point. For the finite difference code, we generate the strain-rate synthetics on a 10 m spacing (for computational reasons) and at frequencies lower than 100 Hz. We use a Gaussian pulse as a source. A comparison of the wavenumber and finite difference synthetics, after compensating for differences in the source frequency, showed comparable results.
Our initial tests use a point source explosion (isotropic) mechanism. This is an approximation of the actual source as a simple point source will not replicate the spectral characteristics of an explosion source (Mueller & Murphy, 1971;Walter & Ford, 2018). However, as this paper primarily focusses on understanding the response of the fiber, we want to investigate how well simple models can match the data before attempting more complex source representations. A secondary issue is that S waves are clearly present in the observed data, which will not be created by a pure explosion source synthetic. To evaluate the possibility of a non-isotropic component to the source we add an arbitrary double-couple strike-slip component. The purpose is to investigate whether the observed S waves can be explained by a source component in terms of timing and approximate amplitude.
As before with the accelerometer data, comparison of the synthetics with the data is not completely straightforward. We compare the observed data with the Z component calculated using a linear approximation to the fiber, as at distances away from the immediate sources region, we expect that most of the seismic energy will be traveling in the Z direction. We begin ( Figure 8) using a subset of the waveforms shown in Figures 6 and 7. The lower part of Figure 8 compares wavenumber synthetics using a pure explosion source with the fiber data. Amplitudes are adjusted using a common empirical scaling for all traces to match. The initial P phase is matched but the pure explosion source synthetics, as expected, fail to generate any S wave energy. We tested adding a 20% double-couple (DC) component using an arbitrary vertical strike-slip mechanism. This creates an S arrival similar to that observed in the data (Figure 8, upper three traces).
Figures 9 and 10 compares the synthetics with the observed data for the entire wavefield for both boreholes (E80 and SW80) using a pure explosion and with a double-couple component added. The synthetics replicate the observed data fairly well except for the difference in timing and amplitude of the reflected P wave and that the waveforms at event level are not well matched, as was noted previously. As the timing of the S waves matches on the synthetics and the data it appears to indicate that they originate from the near-source region rather than a conversion at a velocity interface between the source and the fiber. This demonstrates that S waves energy can be generated at the source region of an explosive source.

Results and Conclusions
This analysis yields a number of results relating to both the performance of the DAS fiber and to the DAG tests. We see that the travel times measured by the fiber are consistent with the available accelerometer data and that the travel times can interpreted to constrain a velocity model but do not uniquely define the velocities. The most notable difficulty was in the time fitting of the surface reflected P wave, which may be caused by anomalous transient spall effects at the surface (or possibly errors in the velocity model). We also corroborate a significant decrease in the seismic wave velocities centered around a depth of 150 m that may have been caused by the large DAG-2 event, as has been noted previously .

10.1029/2021JB022690
17 of 22 Estimates of the signal spectral content between the fiber DAS and acceleration data are similar. As the DAS data should be related to the accelerometer data by a scalar constant, this similarity is expected but is useful to confirm. The DAS data provide a much denser spatial coverage and a systematic difference in frequency content is observed between DAG-1 and DAG-3 event. The cause of this shift is uncertain, but may be related to the decrease in seismic velocities and possible changes in attenuation apparently caused by the large DAG-2 event.
Alternatively, it may be caused by a difference at the source due to changes in coupling between DAG-1 and DAG-3. The data set may be a useful target for a systematic study of attenuation.
The DAS data and the accelerometers show similar waveforms except for the data recorded at the source depth. It was important to carefully sort the accelerometer data to ensure that the anomalous traces that are likely due to instrumental artifacts are eliminated. At most depths, the fiber waveforms matched the longitudinal component of the accelerometer data. The large DAG-2 event may have impacted the later recording of the DAG-3 explosion, possibly by reducing the coupling between the cable and the surrounding alluvium.
Although the waveforms matched closely, we were unable to directly match amplitudes between the DAS strainrate and accelerometer data after applying the expected corrections based on the apparent phase velocity and fiber geometry. In particular, we had difficulty in matching the fiber and accelerometer data at the source depth with accelerometer data for both the waveform shape and absolute amplitude. This may be due to rapid changes in strain rate that exceeded the capabilities of the phase unwrapping algorithms in the interrogator. We also suspect inaccuracies in the metadata, which we attempted to resolve but have not yet been successful. As a check of our understanding of the fiber response and of the source characteristics, we calculate simple synthetics of the DAS data. Two approaches were tested: a wavenumber algorithm and a 3D finite difference algorithm. We approximate the response of the helical fiber with a simple approach that ignores the effect of the cable and treats the fiber as linear. Given the uncertainties in measured amplitudes, we focused on waveform shape rather than absolute amplitudes. For this case, where we believe a 1D velocity model is an adequate representation of the velocity structure, we find that the wavenumber algorithm is a useful approach due to the speed and capability to generate high frequencies. We test both a pure explosion and explosion with partial double-couple added; the source with 20% double couple matches both the observed P wave and the S wave timing. This indicates that the S wave energy emanates from very near the source region. The relative S wave amplitudes at the two boreholes differ, which indicates either a variation in source properties with azimuth or possibly in the lateral velocity structure. While the approximation of the helix as linear is very simple, it seems to be effective at distances away from the event depth.
As construction of strain synthetics has been shown to be fairly straightforward, we believe that this is a preferred approach to modeling DAS data, as opposed to the more roundabout and imprecise method of converting DAS data to velocity (or acceleration) data and then modeling. Care should be taken to evaluate to investigate the possibility of errors in phase unwrapping, as signals that exceed the dynamic range will not display clipping as appears on standard seismic sensors, but will typically show lower amplitude incoherent signals. In summary, we find that the DAS provides high-quality recording of strong ground motions with good waveform fidelity and excellent spatial density. Recovery of high amplitudes requires high sample rates and short gauge lengths. While it was difficult to match exact absolute amplitudes, we suspect that this was due to inaccuracies in metadata and does not reflect fundamental physical limits. Relative amplitudes for the DAS for each data collect appeared internally consistent. The high-resolution spatial data was valuable for resolving the velocity structure. We expect that the DAS data could be useful for source modeling as well, both for the moment tensor and for other parameters, such as stress drop and attenuation.
We conclude that DAS is clearly useful for measuring nearfield explosion data and may be useful for other sources of seismic strong motion such as earthquakes. A DAS fiber deployed in the nearfield of an earthquake would be effective and could provide previously unavailable data.
Future efforts should be expended in understanding the optimal interrogator parameters, the exact response of the helical fiber and to improve modeling capability. The combination of different fiber geometries such as both linear and helical fiber would provide useful constraints. These efforts should emphasize field recordings in carefully controlled settings with calibrated ancillary sensors and combined with modeling.