Evaluating Seismic Beamforming Capabilities of Distributed Acoustic Sensing Arrays

The versatility and cost-efficiency of fibre-optic Distributed Acoustic Sensing (DAS) technologies facilitate geophysical monitoring in environments that were previously inaccessible for instrumentation. Moreover, the spatio-temporal data density permitted by DAS naturally appeals to seismic array processing techniques, such as beamforming for source location. However, the measurement principle of DAS is inherently different from that of conventional seismometers, providing measurements of ground strain rather than ground motion, and so the suitability of traditional seismological methods requires 5 in-depth evaluation. In this study, we evaluate the performance of a DAS array in the task of seismic beamforming, in comparison with a co-located nodal seismometer array. We find that, even though the nodal array achieves excellent performance in localising a regional ML 4.3 earthquake, the DAS array exhibits poor waveform coherence and consequently produces inadequate beamforming results that are dominated by the signatures of shallow scattered waves. We demonstrate that this behaviour is likely inherent to the DAS measurement principle, and so new strategies need to be adopted to tailor array processing tech10 niques to this emerging measurement technology. One strategy demonstrated here is to convert the DAS strain rates to particle velocities by spatial integration using the nodal seismometer recordings as a reference, which dramatically improves waveform coherence and beamforming performance, and warrants new types of “hybrid" array design that combine dense DAS arrays with sparse seismometer arrays.


Introduction
Dense seismometer arrays play a central role in understanding various geological phenomena, including earthquake rupture behaviour (Kiser and Ishii, 2017;Meng et al., 2011), micro-seismicity (Inbal et al., 2016), fault zone structure (Zigone et al., 2019), and deep crustal and mantle geology (Jiang et al., 2018;Lin et al., 2013). Moreover, seismic arrays also serve civil 15 protection purposes through monitoring nuclear test ban treaty violations (Ringdal and Husebye, 1982), monitoring volcano deformation and activity (Inza et al., 2011;Nakamichi et al., 2013), and potentially issuing earthquake early warnings (Meng et al., 2014). While the benefits of seismic arrays are evident, the deployment and maintenance of these arrays is (logistically) costly, and consequently they are often deployed as part of temporary campaigns.
The recent emergence of fibre-optic Distributed Acoustic Sensing (DAS;Hartog, 2017;Zhan, 2020) has opened up a plethora 20 of possibilities and applications in seismic-and transient deformation monitoring. Fibre-optic cables are relatively inexpensive, require little to no maintenance, and can be deployed in environments that were previously impractical for or inaccessible to traditional seismometers, such as urban environments Fang et al., 2020), glaciers and permafrost regions Walter et al., 2020), deep boreholes (Cole et al., 2018;Lellouch et al., 2019), and in lakes and submarine environments (Lindsey et al., 2019;Sladen et al., 2019) -see also Zhan (2020) for a concise review of applications of the PoroTomo nodal seismometer array (black dots) and the fibre-optic cable (red line); c-e) Ground motions recorded on the N/E/Zcomponents of nodal seismometer #97, which is marked in panel b by a green dot; f) Strain rates recorded by a channel co-located with nodal seismometer #97. All waveforms are filtered in a 0.5-10 Hz pass band. the data from the nodal seismometer array, but did not attempt to make a comparison with DAS data. In the present study, we retrieve and interpret the same data set as was analysed by Wang et al. (2018), and so we build upon the conclusions drawn 60 from this previous study.

MUSIC beamforming
Seismic beamforming is a commonly employed array processing technique for estimating the direction-of-arrival (azimuth) and slowness of the seismic waves arriving at a seismic array (Capon et al., 1967;Hutchison and Ghosh, 2017;Krüger et al., 1993). It is assumed in most beamforming applications that the signal recorded at the k-th station in the array can be represented deployed at the surface, thus θ is the azimuth of propagation of the incident wave. Throughout the study we assume a single source (N = 1), so that the frequency-domain representation of the recorded signal can be written as: where e k is the noise recorded at the k-th station, and a k = e iωτ k is the steering vector that dictates the phase shift (time delay) 70 of the signals at each station, relative to the centre of the array. The theoretical time delay τ k = −S (∆x k sin θ + ∆y k cos θ) is computed over a grid of candidate apparent slowness values S and azimuths θ, with a given station location (∆x k , ∆y k ) relative to the centre of the array. In traditional delay-and-sum beamforming, the likelihood of each candidate in the grid of S and θ is estimated as the projection of the steering vector a onto the covariance matrix C 2 , defined as: wherex denotes the complex conjugate of x. Here, the spectra and cross-spectra involved in the equation above are estimated by the multi-taper method (Thomson, 1982), following Meng et al. (2011). Note that the covariance matrix is complex, and that it is scaled by the norms of the waveforms x such that 0 ≤ |C 2 | ≤ 1. Consequently, the magnitude of C 2 is not affected by amplitude differences between x i and x j , e.g. due to spatial variations in coupling or fibre orientation, which could be represented as a station-specific factor α k (ω) multiplying the first term of the right-hand-side of Eq. (1). However, local effects 80 leading to spatial variability of waveform shape are not compensated by this normalization.
MUltiple SIgnal Classification (MUSIC) is an extension of classical beamforming approaches that acknowledges sparsity in the number of signals arriving at the array, resulting in higher-resolution estimates of the back-azimuth and slowness of the seismic waves (Goldstein and Archuleta, 1987;Meng et al., 2011;Schmidt, 1986). Instead of projecting the steering vectors onto the full covariance matrix, a pseudo-power of the signal is estimated as the reciprocal of the projection of the steering 85 vectors onto the noise-space of the covariance matrix, which is found through an eigenvalue decomposition of C 2 . For a detailed exposition of MUSIC, the reader is referred to Schmidt (1986).

Signal characteristics and coherence
Before attempting to perform beamforming on the array data, we first consider the spectral characteristics of the recorded 90 signals -see Fig. 2a and b. The velocity spectra of the three components of the wavefield recorded by the nodal array are first converted into acceleration spectra, which are proportional to the DAS strain rate spectra under the assumption of a single plane wave, with the phase velocity as the proportionality constant. An apparent phase velocity of 0.2 km s −1 for both the Pand S-phases gives a good comparison between the nodal spectra and the DAS spectra, which suggests that a common type of wave (e.g. scattered surface waves) dominates the spectra of these time windows. For reference, the median power of the noise 95 recorded by the nodal array prior to the P-wave arrivals is indicated by a grey band in Fig. 2. The durations of the noise, P-, S-wave d Figure 2. a-b) Median power spectral densities of the nodal and DAS arrays for the P-wave and S-wave. The nodal seismometer recordings are converted into acceleration spectra, which are proportional to the DAS strain rate spectra. The proportionality constant (the apparent phase velocity) is taken as 0.2 km s −1 . The noise floor is shown as a grey band; c-d) Mean power in selected frequency bins as a function of the fibre-optic cable orientation with respect to the back-azimuth of the seismic source. The theoretical sensitivities (Martin et al., 2018) are included for reference. and S-wave windows over which the spectral power was computed are all taken to be 5 s, which is long enough to include low frequency information.
Owing to the nature of the measurement principle of DAS (i.e. measuring strains rather than particle motions), the directional sensitivity of the fibre to P-and S-waves is different from nodal seismometers (Kuvshinov, 2016;Zhan, 2020). For a gauge 100 length that is much smaller than the seismic wavelength, the DAS strain rate is proportional to cos 2 θ for a P-wave or SV-wave, and sin 2θ for an SH-wave, assuming a plane wave with incidence angle θ relative to the fibre (Martin et al., 2018). These theoretical sensitivities are plotted for reference in Fig. 2c   orientation bins. As was also concluded by Wang et al. (2018) from analysing the directional dependence of the signal-to-noise ratio, no directionality of the mean power is observed. Moreover, the variability within a given frequency band exceeds one 105 order of magnitude. Wang et al. (2018) interpreted this as an effect of the heterogeneous site response, which likely exerts a first-order control on the amplitudes and directionality of the ground motions. This will be demonstrated in more detail in Section 3.4.

Beamforming results of the nodal and DAS arrays
To set a baseline, we first beamform the P-and S-waves recorded by the nodal array for each component separately. To visualise 110 the coherence of the wavefield in each direction, we select a 10-second time window starting a few seconds before the P-arrival.
The waveforms are then ordered by distance from the earthquake epicentre, band-pass filtered in the 0.5-1 Hz range, and each trace is scaled by its standard deviation -see Fig. 3. In particular the vertical waveforms exhibit very strong coherence across the entire array. Among the horizontal components, the N-component is more coherent, consistent with a source that is oriented almost directly south of the array (with a back-azimuth of 157°from the centre of the array). Similarly, the S-waves (not shown 115 here) exhibit strong coherence particularly in the E-direction, followed by the N-and Z-directions.  with radial grid increments of 2 km s −1 (indicated by the red numbers). The frequency bands and components are indicated above each column and beside each row, respectively. The true back-azimuth of the source is indicated by a red star in each panel.
The P-wave beamforming results for the nodal array (Fig. 4) show a well-resolved source in the southeast, with an azimuth close to the true back-azimuth of 157°. As expected from the waveform coherence, this source is most stable and well-resolved for the vertical component, with an apparent propagation velocity between 4 and 6 km s −1 . The beams formed from the Ncomponent also indicate a southeast direction-of-arrival, but with a less well-resolved apparent velocity. The beams formed 120 from the E-component suggest weak, poorly resolved sources in the south, west, and east, which are likely scattered P-waves.
The sources indicated by beamforming of the S-wave (Fig. 5) are even better resolved than the P-wave sources, particularly in the E-component, with an apparent propagation velocity between 2 and 4 km s −1 . For an assumed true P-and S-wave speeds of 2.1 and 1.3 km s −1 , respectively (crudely estimated from Feigl and the PoroTomo Team, 2018), the inferred apparent velocities would correspond with an inclination of the direction-of-arrival of 65°(consistent with the ratio of vertical to horizontal 125 amplitudes of the nodal P-waveforms).
In strong contrast to the nodal array, the P-and S-waveforms recorded by the DAS array show a low degree of coherence ( Fig. 6). While individual cable segments may exhibit some internal coherence (analysed further in Section 3.4), this coherence does not persist across the array. There are several factors that contribute to the incoherence of the recorded signals. Firstly, for a horizontal cable, the DAS strain is a combination of gradients of the two horizontal components of the wavefield, which may 130 lead to unfavourable interference. Secondly, the amplitudes of the DAS recordings depend strongly on the coupling of the fibreoptic cable to the ground (Wang et al., 2018) and on the angle of incidence of the incoming plane wave (Martin et al., 2018), so that various segments at different locations and with different fibre orientations experience variable signal-to-noise ratios.
Thirdly, depending on the orientation of the fibre, polarity flips are anticipated (Fang et al., 2020). Lastly, spatial gradients of the particle velocity (i.e. strain rates) are highly sensitive to local heterogeneities (Singh et al., 2020), and their amplitudes are 135 inversely proportional to the apparent phase velocity so that slow waves (often scattered waves) are amplified.
When we nonetheless continue to perform beamforming on the DAS array recordings, we obtain highly variable results ( Fig. 7). At the lower frequencies (below 2 Hz), we find a diffuse spread of pseudo-power over a range of potential source azimuths and apparent velocities. By contrast, at the higher frequencies, we find several well-resolved sources pointing in the southeast and east directions, but with very low apparent velocities (less than 2 km s −1 ). These apparent velocities between 1 140 and 2 km s −1 are consistent with the inferred S-wave speeds at depths of a few hundred metres, suggesting a shallow "source" (most likely a seismic scatterer). As mentioned above, slow phase velocities amplify the recorded DAS strain rates, and so it is not unexpected that, in the absence of strong coherent direct arrivals, slow, scattered waves dominate the beamforming solutions. Also recall that through the definition of the covariance matrix, the absolute amplitude of the recorded signals S-wave Figure 6. P-and S-waveforms recorded by the DAS array, ordered by distance from the earthquake epicentre, filtered in a 0.5-2 Hz pass band, and scaled by the standard deviation of each trace. Time is relative to the start of the recordings, and is synchronised with the waveforms shown in Fig. 1. For clarity, only every 20 th channel waveform is plotted. is irrelevant, which partially addresses the issues of (potential) directionality and ground coupling. However, the wavefield 145 is composed of multiple waves (e.g. direct and scattered arrivals) and their relative phase amplitudes may vary with fibre orientation, which still affects the pseudo-power.

Simulating DAS recordings from the nodal array
In the previous section, we pointed out several potential reasons for the lack of waveform coherence and the inconsistent beamforming results. In the following section, we will explore one of these factors that is inherent to the DAS measurement 150 principle: the one-dimensional strain rate measurement that aggregates multiple components of the particle velocity field. A 9 https://doi.org/10.5194/se-2020-157 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. DAS measurement provides only on component of strain, the longitudinal strain along the direction of the fibre. This limitation can be mitigated by making the measurements along helically-wound cables (Kuvshinov, 2016), but since such cable designs were not deployed during the PoroTomo experiment, one has to resort to alternative approaches.
As has been clearly demonstrated by Wang et al. (2018), the particle velocity measurements recorded on two nodal seis-155 mometers separated by a small distance 2L can be accurately converted into the average longitudinal strain rateε between the two nodes (expressed here at their midpoint x): https://doi.org/10.5194/se-2020-157 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
whereu is the particle velocity in the direction parallel to the positional difference vector between the two nodes. This average strain rate is equal (or proportional) to the DAS strain measured along a gauge length 2L whose end points are co-located with 160 the nodes, if the cable is straight and has a spatially uniform coupling between the two nodes. Given the density of the nodal array, in which most nodes are positioned less than 100 m from its nearest neighbour, we can use this relationship to simulate the response of a DAS array to the strain field induced by the Hawthorne earthquake, and test the effect of superimposing multiple independent components on the beamforming performance. To this end, we triangulate the node coordinates (Fig. 8a), yielding pairs of stations that define each edge in the mesh. For node pairs separated by a distance less than 80 m, we rotate the 165 horizontal (N and E) components of velocity onto the 'virtual' DAS fibre orientation θ asu =u N sin θ +u E cos θ. Substitution ofu into Eq.
(3) yields the strain rate along the simulated DAS fibre at the midpoint between each pair of stations (414 in total): where the subscripts A and B indicate the two seismometers between which the strain rate is calculated, each located at a 170 coordinate point (x, y).
The resulting simulated strain rate P-waveforms are shown in Fig. 8b, for selected segments with an orientation within ±10°from the event back-azimuth (red segments in Fig. 8a). Even though individually the N-and E-components recorded by the nodal stations exhibit some coherence across the array (see Fig. 3), the horizontal strain rates, involving differences of velocities on two horizontal components, are not coherent. Moreover, if the strain rate (the gradient of the particle velocity 175 field) is calculated only on the basis of the strongly coherent N-or Z-components ( Fig. 8c and d, respectively), then the coherence that is seen in the particle velocity measurements (Fig. 3) is almost completely lost. This indicates that it is not just the superposition of two components that causes destructive interference, but that this is caused by the gradient operator itself (regardless whether this operation be done mathematically, like was done here, or physically, like in a DAS fibre). Also, since only segments were selected that are near-parallel to each other, the lack of waveform coherence cannot be attributed to 180 directionality effects.
When we perform the beamforming on the P-waveforms recorded by the virtual DAS array (all segments; Fig. 9), the only sources that stand out are those with very slow apparent phase velocity (< 2 km s −1 ) and azimuths that vary from west to east. Since the overall waveform coherence across the array is low, these sources likely result from subregions in the array that locally exhibit moderate coherence, but which does not persist throughout the array. Owing to the directional sensitivity of the 185 (simulated) DAS measurement, combining segments of different orientations may affect the beamforming results. We repeated the beamforming on selected segments with an orientation ±10°from the event back-azimuth (see Fig. 8a). When only these sub-parallel segments are selected, an ambiguous source arises in the west with high apparent velocity, which is inconsistent with the back-azimuth and apparent phase velocities of the seismic source. This exercise demonstrates that the measurement principle of DAS challenges beamforming methods that are traditionally applied to particle velocities rather than to strain rates.
Since we derived the strain rates directly from the nodal seismometers, the lack of coherence and beam resolution seen in the DAS data (Figs. 6 and 7) cannot be attributed to DAS-specific technicalities like coupling of the DAS cable with the ground or phase unwrapping artefacts, since the nodal seismometers and their derived data do not suffer from this. Parallel segments (v < 2 km/s) Figure 9. Beamforming results of the P-waveforms recorded by the virtual DAS array, derived from the E-and N-components of the nodal seismometers (i.e. Eq. (4)). The panel representation is the same as in Fig. 4. Each second row of panels is a magnification of the upper row, highlighting the sources with low apparent phase velocity (< 2 km s −1 ). The upper two rows are the results of beamforming all the segments, while the lower rows are the results of beamforming only sub-parallel segments (indicated in red in Fig. 8a).

Selective beamforming of the DAS array
Even though the DAS array as a whole does not exhibit strong waveform coherence, there are short individual segments that 195 do exhibit excellent coherence locally. Examples of this can be found in Wang et al. (2018) (their Fig. 14), who selected three segments to estimate the apparent P-wave speed from the first P-wave arrivals recorded by the DAS fibre. The apparent phase velocities obtained from this analysis ranged from 1.124 to 1.452 km s −1 , which are much lower than the apparent velocities obtained from the nodal array beamforming (between 4 and 6 km s −1 ), and are suggestive of a shallow, scattered source.
Since these segments exhibit strong waveform coherence (Fig. 10), we can attempt to form a stable beam by selecting only the 200 channels associated with these segments.  Beamforming on multiple small, linear segments comes with additional challenges. First of all, even though the waveforms may be coherent within one segment, they are not necessarily coherent between segments, which affects the beamforming performance when the waveforms of all the segments are combined. On the other hand, if beamforming is performed on each segment individually, ambiguity arises from the linear geometry of the segments. A plane wave travelling parallel to the length 205 of the cable at a high speed will induce the same phase shift as a plane wave impinging on the cable with a larger angle but at a lower speed. In other words, there exists a perfect trade-off between back-azimuth and apparent velocity for linear segments.
This trade-off is most clearly seen when plotting the beam pseudo-power in a space spanned by the azimuth and apparent slowness (reciprocal apparent velocity), in which the beam pseudo-power will appear as a straight line perpendicular to the orientation of the cable (see Fig. 11; see also Fig. 2e in Lindsey et al. (2019)). This ambiguity can be resolved by combining 210 the beamforming results of multiple segments: the source azimuth and slowness that produces the phase shifts consistent with all segments is the one where the linear bands of beam power intersect. As laid out by Rieken and Fuhrmann (2004), the intersection of the invididual signal spaces of all subarrays can be found by the method of projection onto convex sets. From this analysis, it follows that the signal space of M subarrays combined furnish a block matrix, which upon substitution into the with G m representing the projection onto the noise-subspace of the covariance matrix C 2 m (as defined in Eq. (2)) and a m the steering vector of the m-th subarray, and the superscript H denoting the conjugate transpose. In other words, the analysis of the intersections of signal subspaces of the subarrays naturally leads to a harmonic mean of the MUSIC pseudo-spectra of each subarray. When applying this rationale to the subarray beams presented in Fig. 11, we obtain a combined beam which shows a 220 well-resolved source with a back-azimuth towards the east and an apparent slowness of 1.2 s km −1 (0.83 km s −1 ). Consistent with our previous interpretations, the azimuth and phase velocity of this source are incompatible with direct P-wave arrivals from the Hawthorne earthquake; it is therefore interpreted as a shallow scattered wave.
Since the internal waveform coherence of each segment varies across the array (and also with the selected frequency band), manually selecting a few segments for the beamforming may introduce a bias. However, if internal waveform coherence is the 225 main selection criterion, whether or not to include individual segments in the beamforming analysis can be determined on the basis of the L 2 -norm of the covariance matrix, i.e.: where N is the number of channels in a given segment for which the covariance matrix C ij is computed. We compute c 2 for the P-waveforms of each quasi-linear segment of the cable over the 0.5-1 Hz frequency band, and select those segments with c 2 > 230 0.9 (six in total; Fig. 12a). Since these segments are all linear, we obtain ambiguous results in terms of the azimuth and apparent velocity (which trade-off with one another), but within this ambiguity the sources are well resolved (Fig. 12b). However, while most of the selected segments suggest a direction-of-arrival between the east and the south with a maximum apparent velocity of 2-3 km s −1 , there are two segments that suggest a direction-of-arrival from the north with a maximum apparent velocity at around 1 km s −1 , which most likely signify the predominance of locally scattered waves. When we combine the beam pseudo-235 power of each segment through Eq. (5), we obtain an apparent source with a back-azimuth pointing northeast, and an apparent velocity of 0.7 km s −1 (Fig. 12c). At first, this seems counter-intuitive, as none of the segments seem to suggest a source northeast of the array, while the combined result does. Moreover, the combined apparent velocity is substantially less than the maximum apparent velocities of all the segments individually. However, these results are a direct consequence of the harmonic averaging procedure adopted here, which implictly assumes that the wavefields at all selected subarrays are dominated by the 240 same sets of waves: maximising the MUSIC pseudo-spectra at the intersection of the signal spaces of all the segments becomes problematic for sources that are diametrically opposite of each other (e.g. north-south), as the intersection of the signal spaces occurs only at infinity when representing the pseudo-power in slowness space (or at the origin in velocity space). This then leads to the diffuse spread of pseudo-power seen in Fig. 12c with the maximum pseudo-power at high slowness values.
When we exclude the two segments with an apparent source north of the array, we obtain a maximum beam power at an 245 azimuth of 122°and apparent velocity of 2.0 km s −1 . While this result is in better agreement with the true direction of arrival of the seismic source (157°with an apparent P-wave speed of 4-6 km s −1 , as inferred from the nodal array), the discrepancy is substantial. Particularly the low apparent velocity inferred from the DAS segments is suggestive of the arrival of scattered waves rather than direct waves. Regardless, this result is only obtained after manual quality control and selection of desired  segments. For automated segment selection and beamforming one must be cautious for conflicting directions-of-arrival that 250 lead to artificial results like in Fig. 12c.

Importance of site heterogeneities and seismic wave scattering
The purpose of beamforming is to locate the origin of the incoming signals, which requires that the signals be direct arrivals.
As the seismic equivalent of acoustic echoes, scattered waves will follow a trajectory different from the direct arrivals, exhibit 255 a different direction-of-arrival, and so appear to originate from a different source location. Throughout this work, we have shown that the DAS array recordings at the Brady Hot Springs natural laboratory are strongly affected by scattered waves and local heterogeneities. However, the nodal seismometer array seems much less affected, and displays strong signatures of the Hawthorne earthquake's direct arrivals. To summarise, the likely reasons for this are as follows: 1. For P-and S-waves, DAS strain rate measurements are most sensitive to strains parallel to the direction of the fibre 260 (Martin et al., 2018). Direct arrivals originating from a distant source arrive at the array at a steep inclination, and so their projection onto the direction of a horizontal fibre is comparatively small (for P-waves, the DAS sensitivity exhibits a cos 2 decay with inclination). By contrast, free-surface topography and shallow subsurface heterogeneities may cause scattered seismic waves to arrive at the array at a shallow inclination, so that these are greatly amplified in the DAS measurements.

265
2. As pointed out by Daley et al. (2016), the relation between the strain ε induced by a plane wave and the particle velocitẏ u in the along-fibre direction is ε = ± 1 cu , with c being the apparent phase velocity of the medium. This implies that apparent fast waves (i.e. arriving at steep inclinations) are damped with respect to sub-horizontally-travelling waves such as those generated by shallow scattering.
3. From the numerical simulations and theoretical analysis of Singh et al. (2020), it is immediately clear that gradients in 270 the particle velocity field are highly sensitive to heterogeneities. Since the subsurface beneath the Brady Hot Springs site is strongly heterogeneous (Feigl and the PoroTomo Team, 2018), spatial variations in the phase velocity are expected to exert a strong influence on the DAS measurements. The lack of observable directional sensitivity of the DAS array (our Fig. 2 and Fig. 22 of Wang et al. (2018)) further attest to this. Singh et al. (2020) proposed that a correction term (the "coupling tensor") may be inverted for, so that the amplitudes of the recordings may be adjusted to represent the 275 "true" velocity gradients. Since the true particle velocity field is known (being recorded by the nodal array), this offers an interesting perspective for future analysis of the Brady Hot Springs data set.
4. Lastly, even when considering the directional sensitivity of DAS, the locally recorded strain rates appear to be highly incoherent, while the particle velocity measurements themselves exhibit very strong coherence (compare Fig. 3 with  Fig. 8). The superposition of multiple orthogonal components of the particle velocity field may lead to additional de-280 structive interference, even though this was not clearly seen in our analysis in Section 3.3 (see Fig. 8b).
As summarised by Shearer (2015), the effects of phase velocity heterogeneities are most pronounced when the size of the heterogeneity is similar to the seismic wavelength. Owing to a scarcity of large-scale heterogeneities, low frequency signals may be less affected by seismic scattering. Coincidentally, the DAS beamforming results in the lowest frequency band (0.5-1 Hz; Fig. 7) do suggest a source with an azimuth that corresponds with the true back-azimuth, although the spread in the 285 beam pseudo-power is diffuse. This may be a consequence of the signal-to-noise ratio in the lowest frequency range, since the source exhibits low spectral power in this range (Fig. 2). Notwithstanding, DAS has a flat frequency response even at very low frequencies (Lindsey et al., 2020), so that further investigations of DAS beamforming at low frequencies are warranted.

Implications for beamforming on sparse and dense DAS arrays
The analysis of the PoroTomo experiment has revealed some limitations of conventional beamforming methods applied to DAS 290 array data, particularly in relation to scattered waves and heterogeneities. However, the Brady Hot Springs geothermal field may be considered a particularly unfavourable setting for DAS seismic beamforming owing to the complexity of the subsurface.
Fibre-optical cables deployed on more homogeneous bedrock terranes may not suffer as much from high-amplitude shallow scattering. On the other hand, one of the main promises of DAS is its versatility in deployment conditions, with interesting deployment targets including "heterogeneous" environments such as urban areas Fang et al., 2020) and 295 submarine basins (Lindsey et al., 2019;Sladen et al., 2019). For many civil monitoring applications, such as traffic density monitoring (Liu et al., 2018) and vehicle tracking (Wiesmeyr et al., 2020), some of the issues pointed out in the previous section do not apply, as the signals of interest arrive at the DAS fibre at a shallow (or zero) inclination. However, for the purpose of localising deep or distant sources, the inclination sensitivity of DAS starts to become directly relevant.
A second unfavourable aspect of the PoroTomo dataset is that the DAS array is deployed within a relatively small region 300 (1500 by 500 m), which limits the resolution of beamforming methods (being proportional to the span of the array). Sparse L-shaped and quasi-linear array configurations provide a much larger array span, at the expense of increased source azimuth ambiguity inherent to linear arrays. As was done in Section 3.4, multiple segments of variable orientations can be combined to resolve this ambiguity. Moreover, the same procedure can be adopted to extract the signals carried by direct arrivals: in the case of the Brady Hot Springs geothermal site, the entire array receives seismic energy from (potentially) multiple nearby 305 scattering sites, obscuring the direct arrivals. By contrast, sparse arrays that extend over long distances may receive seismic energy from different scattering sites along the trace of the cable. By selecting and combining the beam power of several segments following the harmonic averaging method proposed by Rieken and Fuhrmann (2004), sources that are common to all segments are amplified with respect to segment-specific sources (i.e. scatterers), provided that the direct arrivals exhibit a sufficiently large footprint in all segments. Large sparse arrays may therefore be more suitable for seismic beamforming than 310 compact dense arrays.
One key assumption that underlies beamforming is that the signal is carried by a plane wave. This assumption implicitly requires that the source be distant compared to the extent of the array. Moreover, the phase velocity is assumed to be uniform across the array. Both these assumptions are embedded in the definition of the steering vectors through the time delay τ (see Eq. (1)). Since DAS on fibre-optic cables of several tens of kilometres in length has been demonstrated to be feasible (Lindsey 315 et al., 2019;Sladen et al., 2019), these assumptions may break down for local and regional seismic sources. Moreover, for e.g. earthquake early warning purposes, fault zones may be instrumented with fibre-optic cables running along-strike for many tens of kilometres, so that the finite extent of the rupture and rupture complexity will prevent the application of beamforming methods. Fortunately, the plane wave and uniform velocity assumptions can be relaxed by directly computing the travel time between a candidate source location and a location along the cable, turning the beamforming problem into a back-projection 320 exercise (Kiser and Ishii, 2017). Alternatively, the azimuth estimates derived from beamforming on individual segments can be combined to produce a triangulated source location (e.g. Hutchison and Ghosh, 2017;Stipčević et al., 2017). In both cases, the source localisation results will greatly benefit from the large lateral extent of the DAS arrays.

Conclusions
This study considered the potential of fibre-optic Distributed Acoustic Sensing (DAS) arrays for the purpose of seismic beam-325 forming. This was done by performing beamforming on the ground motions generated by the March 2016 M L 4.3 Hawthorne earthquake, as recorded by a DAS array co-located with a dense nodal seismometer array at the Brady Hot Springs geothermal field. Comparing the waveforms recorded by DAS with those recorded by the nodal seismometers, we find that the strong waveform coherence of the nodal array is absent in the DAS array. Since the quality of the beamforming results depends strongly on waveform coherence, the DAS array is unable to produce a robust source azimuth and apparent velocity, whereas 330 the nodal array produces an extremely well-resolved source location that is consistent with the true back-azimuth of the earthquake epicentre. Instead, beamforming on the DAS array reveals source locations that likely correspond with shallow seismic scattering sites. We attribute the lack of DAS waveform coherence to the DAS measurement principle, which inherently leads to diminished sensitivity of DAS recordings to the direct arrivals, supposedly arriving at the array at a steep inclination, and amplifies scattered waves arriving at shallow inclinations. Moreover, we demonstrate that the spatial gradients of the particle 335 velocity field (i.e. strain rate) exhibit far lower coherence than do the particle velocity waveforms, which additionally impedes beamforming. Compared to other DAS arrays, this may be aggravated by the strong phase velocity heterogeneities present at the Brady Hot Springs geothermal field. Nonetheless, many of these issues may be resolved by expanding the extent of the array, and through the combination of the beamforming results of individual segments along the fibre-optic cable.
Code and data availability. Python scripts that reproduces the results and figures in this manuscript are available at