the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Evaluating seismic beamforming capabilities of distributed acoustic sensing arrays
Martijn P. A. van den Ende
JeanPaul Ampuero
The versatility and cost efficiency of fibreoptic distributed acoustic sensing (DAS) technologies facilitate geophysical monitoring in environments that were previously inaccessible for instrumentation. Moreover, the spatiotemporal 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 indepth evaluation. In this study, we evaluate the performance of a DAS array in the task of seismic beamforming, in comparison with a colocated nodal seismometer array. We find that, even though the nodal array achieves excellent performance in localising a regional M_{L} 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 techniques 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.
Dense seismometer arrays play a central role in understanding various geological phenomena, including earthquake rupture behaviour (Kiser and Ishii, 2017; Meng et al., 2011), microseismicity (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 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 fibreoptic distributed acoustic sensing (DAS; Hartog, 2017; Zhan, 2020) has opened up a plethora of possibilities and applications in seismic and transient deformation monitoring. Fibreoptic 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 (Dou et al., 2017; Fang et al., 2020), glaciers and permafrost regions (AjoFranklin et al., 2017; 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 in geosciences. DAS thus has an enormous potential to complement or replace seismometer arrays (Jousset et al., 2018). However, the measurement principles of DAS are inherently different from those of conventional seismometers (Zhan, 2020), which presents new challenges in interpreting DAS data. Traditional array processing techniques, such as seismic beamforming, need to be reevaluated for the application to DAS data.
Even though several studies already reported first results of applying seismic beamforming to linear and Lshaped DAS arrays (Fang et al., 2020; Lindsey et al., 2017, 2019), the potential of DAS in beamforming requires further exploration. In this study, we directly compare beamforming results of data from a nodal seismometer array and from a colocated optical fibre cable at the Brady Hot Springs geothermal site, Nevada, USA (Feigl and the PoroTomo Team, 2018). Specifically, we analyse the recordings of the March 2016 M_{L} 4.3 Hawthorne earthquake, which occurred 150 km south of the Brady Hot Springs site and was well captured by both the nodal and DAS arrays. The comparison suggests that the beamforming of the DASrecorded waveforms is severely hampered by shallow seismic scattering and by spatial variations in phase velocities, to which DAS measurements are highly sensitive. This is consistent with previous theoretical findings that ground motion gradients (or strains) are more severely affected by smallscale heterogeneities than ground motions themselves. To remedy this, we propose a method in which we convert the DAS strain rates to particle velocities, yielding DAS beamforming results that are on par with those of the nodal array. We conclude by putting these observations in a broader context of beamforming capabilities of DAS arrays of larger aperture, and their application in seismic source monitoring and earthquake early warning.
2.1 The PoroTomo experiment
The Poroelastic Tomography (PoroTomo) project is a hydrogeological experiment conducted in March 2016 (phase II) at a geothermal site near Brady Hot Springs, Nevada, USA (Feigl and the PoroTomo Team, 2018) – see Fig. 1. For the purpose of highresolution monitoring of changes in rockmechanical properties during operation of the enhanced geothermal system, an array of 238 Fairfield Nodal ZLand 3C seismometers was deployed over an area spanning 1500 by 500 m, as well as several fibreoptic cables for distributed acoustic sensing and distributed temperature sensing. These fibreoptic cables were laid out horizontally in a trench of 8700 m in total length and 0.5 m in depth, and vertically in a borehole down to 400 m. The gauge length was taken to be 10 m, which was supersampled to give a channel spacing of 1 m (i.e. one strain rate measurement was made every 1 m). The geothermal reservoir of Tertiary volcanic rocks is overlain by a thick alluvium of several hundreds of metres in thickness (Jolie et al., 2015). The nearsurface velocity structure of the site has been inferred from the analysis of highfrequency vibroseis sweeps and from noise correlation functions (Feigl and the PoroTomo Team, 2018), showing strong variations over distances of tens of metres.
During the experiment, on 21 March 2016 at 07:37:10 UTC, an M_{L} 4.3 strikeslip earthquake occurred 150 km SSE of the geothermal site at a depth of 9.9 km. The ground motion data of this event recorded by both the nodal and DAS arrays are available at the National Geothermal Data Repository (Feigl, 2016a, b). For convenience, we downsample the nodal and DAS data to 100 Hz, which is largely sufficient for the frequency bands selected in our analysis (up to 5 Hz). The nodal data are corrected for the geophone instrument response (damping factor of 0.7 at a 5 Hz corner frequency). An indepth analysis of the ground motions in terms of frequency content, spatial variability of signaltonoise ratios, and the comparison between DAS and nodal seismometers was performed by Wang et al. (2018). These authors also performed a preliminary beamforming analysis using 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 dataset as was analysed by Wang et al. (2018), and so we build upon the conclusions drawn from this previous study.
2.2 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 kth station in the array can be represented by a superposition of N plane waves, each carrying a signal s and impinging on the array at an angle θ. We consider arrays 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 frequencydomain representation of the recorded signal can be written as follows:
where e_{k} is the noise recorded at the kth station, and ${\mathit{a}}_{k}={e}^{i\mathit{\omega}{\mathit{\tau}}_{k}}$ is the steering vector that dictates the phase shift (time delay) of the signals at each station, relative to the centre of the array. The theoretical time delay ${\mathit{\tau}}_{k}=S\left(\mathrm{\Delta}{x}_{k}\mathrm{sin}\mathit{\theta}+\mathrm{\Delta}{y}_{k}\mathrm{cos}\mathit{\theta}\right)$ 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 delayandsum 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 follows:
where $\stackrel{\mathrm{\u203e}}{x}$ denotes the complex conjugate of x. Here, the spectra and crossspectra involved in the equation above are estimated by the multitaper 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 $\mathrm{0}\le \left{\mathbf{C}}^{\mathrm{2}}\right\le \mathrm{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 stationspecific factor α_{k}(ω) multiplying the first term of the righthand side of Eq. (1). However, local effects leading to spatial variability of waveform shape are not compensated by this normalisation.
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 higherresolution estimates of the backazimuth 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 pseudopower of the signal is estimated as the reciprocal of the projection of the steering vectors onto the noise space of the covariance matrix, which is found through an eigenvalue decomposition of C^{2}. The procedure of estimating C^{2} is as described above, and so the sole difference between MUSIC and traditional beamforming lies in the projection of the steering vectors onto the noise space (and taking the reciprocal), rather than projecting onto the full space of C^{2}. For a detailed exposition of MUSIC, the reader is referred to Schmidt (1986).
3.1 Signal characteristics and coherence
Before attempting to perform beamforming on the array data, we first consider the spectral characteristics of the recorded 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 P and 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 recorded by the nodal array prior to the Pwave arrivals is indicated by a grey line in Fig. 2. The durations of the noise and P and Swave windows over which the spectral power was computed are all taken to be 5 s, which is long enough to include lowfrequency 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 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 and d, alongside the mean power measured within selected fibre orientation bins. As was also concluded by Wang et al. (2018) from analysing the directional dependence of the signaltonoise ratio, no directionality of the mean power is observed. Moreover, the variability within a given frequency band exceeds one order of magnitude. Wang et al. (2018) interpreted this as an effect of the heterogeneous site response, which likely exerts a firstorder control on the amplitudes and directionality of the ground motions. This will be demonstrated in more detail in Sect. 3.4.
3.2 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. We take a time window from 2 s before to 8 s after the first arrival of each respective phase (i.e. 10 s in total). To visualise the coherence of the wavefield in each direction, we select a 10 s time window starting near the P arrival. The waveforms are then ordered by distance from the earthquake epicentre and bandpass 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 backazimuth of 157^{∘} from the centre of the array). Similarly, the S waves (not shown here) exhibit strong coherence particularly in the E direction, followed by the N and Z directions.
The Pwave beamforming results using all the nodes in the nodal array show a wellresolved source in the southeast, with an azimuth close to the true backazimuth of 157^{∘} (Fig. 4). As expected from the waveform coherence, this source is most stable and wellresolved for the vertical component, with an apparent propagation velocity between 4 and 6 km s^{−1}. Only in the 0.5–1 Hz frequency band does the beamforming of this component lead to a relatively poorly resolved location, which may be due to the influence of the corner frequency (typically around 1–3 Hz for an M_{L} 4.3 event; see for example Scholz, 2019). The beams formed from the N component also indicate a southeast direction of arrival, but with a less wellresolved apparent velocity. The beams formed 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 Pwave sources, particularly in the E component, with an apparent propagation velocity between 2 and 4 km s^{−1}. For 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 amplitudes of the nodal P waveforms measured across the nodal array).
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 Sect. 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 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 signaltonoise ratios. Thirdly, depending on the orientation of the fibre, Swave polarity flips are anticipated (Fang et al., 2020). These polarity flips are due to the projection of the particle motion onto the fibre, leading to contraction in some segments and extension in others (Lindsey et al., 2017). 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 inversely proportional to the apparent phase velocity so that slow waves (often scattered waves) are amplified.
When we nonetheless continue to perform beamforming on all DAS array recordings (8621 channels in total), we obtain highly variable results (Fig. 7) for the same window length and frequency range as was used for the nodal array. At the lower frequencies (below 2 Hz), we find a diffuse spread of pseudopower over a range of potential source azimuths and apparent velocities. By contrast, at the higher frequencies, we find several wellresolved 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 and 2 km s^{−1} are consistent with the inferred Swave 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 is irrelevant, which partially addresses the issues of (potential) directionality and ground coupling. However, the wavefield 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 pseudopower.
3.3 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 principle: the onedimensional strain rate measurement that aggregates multiple components of the particle velocity field. A DAS measurement provides only one 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 seismometers separated by a small distance L can be accurately converted into the average longitudinal strain rate $\dot{\mathit{\epsilon}}$ between the two nodes (expressed here at their midpoint x):
where $\dot{u}$ 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 L whose end points are colocated with the nodes, if the cable is straight and has a spatially uniform coupling between the two nodes. A similar relation holds when the distance between the nodes is a multiple of the DAS gauge length (see Wang et al., 2018, their Eq. 5).
Given the density of the nodal array, in which most nodes are positioned less than 100 m from their 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 horizontal (N and E) components of velocity onto the “virtual” DAS fibre orientation θ (relative to east) as $\dot{u}={\dot{u}}^{N}\mathrm{sin}\mathit{\theta}+{\dot{u}}^{E}\mathrm{cos}\mathit{\theta}$. Substitution of $\dot{u}$ into Eq. (3) yields the mean strain rate along the simulated DAS fibre in 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 coordinate point (x,y).
The resulting simulated strain rate P waveforms are shown in Fig. 8b, for selected segments with an orientation within $\pm \mathrm{10}{}^{\circ}$ from the event backazimuth (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 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 is done mathematically, like was done here, or physically, like in a DAS fibre). Also, since only segments were selected that are nearly parallel to each other, the lack of waveform coherence cannot be attributed to 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 it does not persist throughout the array. Owing to the directional sensitivity of the (simulated) DAS measurement, combining segments of different orientations may affect the beamforming results. We repeated the beamforming on selected segments with an orientation $\pm \mathrm{10}{}^{\circ}$ from the event backazimuth (see Fig. 8a). When only these subparallel segments are selected, an ambiguous source arises in the west with high apparent velocity, which is inconsistent with the backazimuth 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 DASspecific 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.
3.4 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 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 Pwave speed from the first Pwave 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 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 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 tradeoff between backazimuth and apparent velocity for linear segments. This tradeoff is most clearly seen when plotting the beam pseudopower in a space spanned by the azimuth and apparent slowness (reciprocal apparent velocity), in which the beam pseudopower 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 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 individual 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 definition of the MUSIC spectrum yields the following (Rieken and Fuhrmann, 2004):
with G_{m} representing the projection onto the noise subspace of the covariance matrix ${\mathbf{C}}_{m}^{\mathrm{2}}$ (as defined in Eq. 2) and a_{m} the steering vector of the mth 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 pseudospectra of each subarray. When applying this rationale to the subarray beams presented in Fig. 11, we obtain a combined beam which shows a wellresolved source with a backazimuth 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 Pwave 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 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 quasilinear segment of the cable over the 0.5–1 Hz frequency band and select those segments with c^{2}>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 segments 1, 2, 3, and 5 suggest a direction of arrival between the east and the south with a maximum apparent velocity of 2–3 km s^{−1}, the other two segments (4 and 6) 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 pseudopower of each segment through Eq. (5), we obtain an apparent source with a backazimuth pointing northeast, and an apparent velocity of 0.7 km s^{−1} (Fig. 12c). At first, this seems counterintuitive, 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 implicitly assumes that the wavefields at all selected subarrays are dominated by the same sets of waves: maximising the MUSIC pseudospectra 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 pseudopower in slowness space (or at the origin in velocity space). This then leads to the diffuse spread of pseudopower seen in Fig. 12c with the maximum pseudopower at high slowness values.
When we exclude the two segments with an apparent source north of the array (so that four segments remain), we obtain a maximum beam power at an azimuth of 122^{∘} and apparent velocity of 2 km s^{−1}. While this result is in better agreement with the true direction of arrival of the seismic source (157^{∘} with an apparent Pwave 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 of conflicting directions of arrival that lead to artificial results like in Fig. 12c.
3.5 Converting DAS strain rates to velocities
As demonstrated in Sect. 3.3, one can convert the nodal seismometer data into strain rate using Eq. (4). In doing so, much of the waveform coherence is lost, and slow phases (e.g. surface waves) are emphasised in the beamforming results. Following the same reasoning, the converse is expected to be true: to alleviate the effects of heterogeneities and slow surface waves as discussed in the previous sections, one could integrate the DAS strain rate data along the cable to obtain particle velocities. As will be shown in the following section, performing this conversion greatly improves the DAS beamforming results. However, great care must be taken in performing this analysis, and so we will start with a detailed description of the adopted procedure.
We begin by taking Eq. (3) and expressing the particle velocity at a point x as follows:
where $\dot{u}$ is the particle velocity in the direction parallel to the fibre (defined as positive when pointing from the reference point x_{ref} to x), Δx is the distance between x_{ref} and x, and $\dot{\mathit{\epsilon}}(x;{x}_{\mathrm{ref}})$ is the average strain rate between x and x_{ref} (taken positive in extension). We emphasise that $\dot{\mathit{\epsilon}}$ is not a finitedifference approximation of the strain at the midpoint of the segment, but the exact average of the strain rate along the segment, and that DAS measures these strain rates averaged over onegaugelengthlong segments. Moreover, no assumptions are made regarding the apparent propagation speed of the signals to convert strain rate into particle motion (as is typically done in other methods; e.g. Lior et al., 2021; Zhu et al., 2021, this issue). Therefore, in the case where a seismometer is colocated with the DAS cable, we can take the horizontal component of the recorded particle velocity in the direction of the fibre as a starting point and compute the particle velocity at the next DAS channel through integration of the recorded strain rate. This procedure can then be extended to obtain the particle velocity at each subsequent gauge length along a straight fibre by using the following equation for a segment of length nL (L being the gauge length):
Naturally this approach will accumulate integration errors as the strain rate is integrated farther and farther away from the seismometer owing to departures from the assumptions of a straight cable and uniform coupling, as well as instrument noise. We can evaluate this to some extent by considering DAS segments that have two colocated nodes, integrating from one node to the other and comparing the converted DAS particle velocity with the terminal node. It is critical to note that DAS records only one component of strain, namely the longitudinal strain associated with the local direction of the cable. Integrating the DAS recordings will therefore only yield one component of the particle motion parallel to the fibre, which cannot be decomposed into perpendicular components (like those recorded by the nodal seismometers). Nevertheless, the singledirection measurement of particle velocity is likely more coherent than the singledirection measurement of strain rate.
The integration along the cable is systematically performed in the direction from low channel numbers to high channel numbers. In the current dataset, this corresponds to the direction running from the DAS interrogator to the opposite end of the cable. Let $\dot{\mathit{u}}\left({\mathit{x}}_{i}\right)$ denote the threecomponent ground velocity field at location x_{i} along the cable, with subscript i denoting the channel index. The longitudinal (alongcable) strain rate recorded by DAS on the channel corresponds to the gaugelong segment $\left[{\mathit{x}}_{i},{\mathit{x}}_{i+\mathrm{1}}\right]$ is $\dot{\mathit{\epsilon}}({\mathit{x}}_{i+\mathrm{1}},{\mathit{x}}_{i})=\frac{\mathrm{1}}{L}\left(\dot{\mathit{u}}\left({\mathit{x}}_{i+\mathrm{1}}\right)\dot{\mathit{u}}\left({\mathit{x}}_{i}\right)\right)\cdot {\mathit{n}}_{i}$, where ${\mathit{n}}_{i}=\frac{\mathrm{1}}{L}\left({\mathit{x}}_{i+\mathrm{1}}{\mathit{x}}_{i}\right)$ is the unit vector pointing from x_{i} to x_{i+1}. In the local reference frame associated with the cable, $\dot{u}$ in Eq. (8) is defined as the projection of the particle velocity in the direction of increasing channel numbers along the cable: $\dot{u}=\dot{\mathit{u}}\cdot \mathit{n}$. Since the cable curves and is laid out in a zigzag pattern, even parallel cable segments may have unit vectors pointing in opposite directions. To bring all velocities to a common reference frame to be used for beamforming, after evaluating Eq. (8) for each segment we multiply $\dot{u}$ by sgn[e⋅n], where e is a reference unit vector taken here as the unit vector pointing east. If this sign correction is not performed, the integrated waveforms will change polarity depending on the orientation of the fibre, and coherence will be broken as a result.
Lastly, in the PoroTomo experiment the recordings of strain rate were made every 1 m, i.e. the gauge length of 10 m was oversampled by a factor of 10. However, owing to the averaging effect of the gauge these additional samples are not independent of one another, and so only one independent measurement of the strain field is made every gauge length. We therefore perform the integration over nonoverlapping gauges, i.e. over channels that are one gauge length away from each other, as indicated in Eq. (8).
Keeping these notations in mind, we perform the conversion of DAS strain rates to particle velocities as described above. We identify the nodal seismometers that are at most 1 gauge length (10 m) away from the nearest DAS channel. We then select linear portions of the DAS array that have a node at the start and at the end of the segment. For each segment we compute the horizontal component of the wavefield (parallel to the segment) as recorded by the starting node and perform the integration along the fibre in the direction of the terminal node (accounting for the polarity as described above). For the purpose of beamforming, we select the DAS segments that are within $\pm \mathrm{10}{}^{\circ}$ of east (see Fig. 13a), which should exhibit a strong sensitivity to the horizontal component of the S wave. To demonstrate that the proposed integration method is accurate, we also select the longest quasilinear segment in the DAS array (cyan segment in Fig. 13a) for inspection of the waveforms. As can be seen in Fig. 13b and c, the particle velocity obtained from the integration of the DAS strain rates up to the last channel is practically identical to the waveforms recorded by the terminal nodal seismometer, demonstrating the accuracy of the method. The comparison is shown in Fig. 13 up to 2 Hz, but the accuracy persists up to 5 Hz. This segment is 230 m in length, and consequently the wavefield changes substantially from the start to the end of the segment (particularly for the S wave; Fig. 13c). Nonetheless, the difference between the converted DAS and the nodal data is minimal, which is encouraging for the integration of DAS strain rates over even longer distances (potentially up to several kilometres).
As expected, converting the DAS strain rates to particle velocities restores the waveform coherence throughout the array (compare Fig. 13d and e). Particularly in the 0.5–1 Hz frequency band the Swave coherence is very strong, and the converted DAS data compare onetoone with the nodal seismometers at the start and end of each selected segment. This similarity and waveform coherence clearly manifests itself in the beamforming (Fig. 13f–h): while the original DAS strain rate waveforms yield ambiguous and complicated beamforming results, the converted DAS and nodal Swave beamforming results are practically identical, suggesting a wellresolved backazimuth close to the true backazimuth of the seismic source and an apparent velocity of 3 km s^{−1}, consistent with previous results (Fig. 5). Hence, we can conclude from this exercise that converting the DAS strain rate to particle velocity alleviates all of the issues associated with the measurement principle pointed out in the previous sections.
4.1 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 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:

For P and S waves, DAS strain rate measurements are most sensitive to strains parallel to the direction of the fibre (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, freesurface 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.

As pointed out by Daley et al. (2016), the relation between the strain ε induced by a plane wave and the particle velocity $\dot{u}$ in the alongfibre direction is $\mathit{\epsilon}=\pm \frac{\mathrm{1}}{c}\dot{u}$, 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 subhorizontallytravelling waves such as those generated by shallow scattering.

From the numerical simulations and theoretical analysis of Singh et al. (2020), it is immediately clear that gradients in the particle velocity field are highly sensitive to heterogeneities. It has also been observed in realworld DAS experiments (e.g. Jousset et al., 2018; Lindsey et al., 2019) that shallow subsurface features such as faults clearly manifest themselves in the DAS records, further attesting to the sensitivity of DAS strain (rate) measurements 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 Figs. 2 and 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 “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 dataset.

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 destructive interference, even though this was not clearly seen in our analysis in Sect. 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 largescale heterogeneities, lowfrequency signals may be less affected by seismic scattering. Coincidentally, the DAS beamforming results in the lowestfrequency band (0.5–1 Hz; Fig. 7) do suggest a source with an azimuth that corresponds with the true backazimuth, although the spread in the beam pseudopower is diffuse. This may be a consequence of the signaltonoise ratio in the lowest frequency range, since the source exhibits low spectral power in this range (Fig. 2). Nevertheless, DAS has a flat frequency response in strain even at very low frequencies (Lindsey et al., 2020; Paitz et al., 2020), so that further investigations of DAS beamforming at low frequencies are warranted.
4.2 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 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. Fibreoptical cables deployed on more homogeneous bedrock terranes may not suffer as much from highamplitude 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 (Dou et al., 2017; Fang et al., 2020) and 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 (1500 by 500 m), which limits the resolution of beamforming methods (being proportional to the span of the array). Sparse Lshaped and quasilinear array configurations provide a much larger array span, at the expense of increased source azimuth ambiguity inherent to linear arrays. As was done in Sect. 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 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 segmentspecific 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 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 fibreoptic cables of several tens of kilometres in length has been demonstrated to be feasible (Lindsey et al., 2019; Sladen et al., 2019), these assumptions may break down for local and regional seismic sources. Moreover, for earthquake early warning purposes for example, fault zones may be instrumented with fibreoptic cables running alongstrike 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 backprojection exercise (Kiser and Ishii, 2017; Zhu and Stensrud, 2019). 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.
4.3 Hybrid array design
As shown in Sect. 3.5, converting the DAS strain rates to particle velocities through integration dramatically increases the coherence across the array and closes the gap between DAS and seismometer arrays in terms of beamforming performance. While at least one seismometer needs to be colocated with the DAS array, this approach motivates the design of “hybrid” arrays that comprise both DAS and seismometers. We here propose three such hybrid arrays (Fig. 14) that optimally employ the available seismic stations while reducing the length of DAS cable used in the array:

A minimal array configuration that records the two independent horizontal components of the particle velocity field features a single seismometer at the centre of a number4shaped array (Fig. 14a). The strain rates recorded by DAS are integrated in the directions away from the central seismometer to give onecomponent measurements of the wavefield in the east–west or north–south direction. The diagonal segment that connects the extremities of the array is not directly connected to a reference station, and so it cannot be integrated to yield particle velocities in this direction. Nonetheless, this connecting segment can be identified and ignored in the analysis, focusing solely on the orthogonal segments connected to the central seismometer. The total length of cable used in this configuration is $\left(\mathrm{4}+\sqrt{\mathrm{2}}\right)L$, where L is the segment length from the centre to one of the corners of the array.

Since DAS measurements have a distinct directional sensitivity to various types of ground motions (Martin et al., 2018), two independent directions of integration may not be sufficient for the analysis of seismic sources of arbitrary backazimuth. To improve the likelihood of having a DAS segment that is optimally oriented for a given backazimuth, a single DAS cable can be deployed to form multiple segments at fixed orientation increments (Fig. 14b). This umbrellashaped array is a generalisation of the 4shaped array shown in Fig. 14a, and it could feature an arbitrary number of radial segments. The example shown in Fig. 14b features six independent segment orientations. The total DAS cable length for an umbrella array with n independent segments (i.e. n=2 in Fig. 14a and n=6 in Fig. 14b) is approximately $\mathrm{2}L\left(n+\mathit{\pi}\left[n\mathrm{1}\right]/\mathrm{2}n\right)$.

In the array configurations using only a single seismometer, integration along each DAS segment yields at most one component of the particle velocity field for that given segment. When multiple seismometers are available, the DAS cable can be deployed in a gridlike fashion with the seismometers positioned along the diagonal of the grid – see Fig. 14c. At the locations where two DAS segments intersect at right angles, the two independent horizontal components of the particle velocity field can be inferred. While other configurations of seismometer placement are possible (e.g. along the boundaries of the checkerboard array), placing the seismometers along a diagonal yields the optimal number of DAS cable intersections per seismometer (the total number of DAS intersections is given as n^{2}−n, with n being the number of seismometers). Correspondingly, the total length of DAS cable required for this array design is 2Δx(n^{2}−1), Δx being the grid spacing.
Owing to attenuation of the scattered light intensity with increasing cable length, standard commercial fibreoptic cables permit a maximum sensing distance of roughly 50 km. Taking this as the maximum length L of the cable available for the array, and taking a grid spacing of 200 m, a checkerboard grid can be constructed with $n=\mathrm{floor}\left(\sqrt{\frac{L}{\mathrm{2}\mathrm{\Delta}x}+\mathrm{1}}\right)=\mathrm{11}$ seismometers. Consequently, this array will have a total span of $\left(n\mathrm{1}\right)\mathrm{\Delta}x=\mathrm{2}$ km, with n^{2}=121 twocomponent measurements every 200 m. For a grid spacing of Δx=20 m (ignoring the limitation of the gauge length for the moment), the corresponding number of twocomponent measurements and total array span are 35^{2}=1225 and 700 m, respectively. For comparison, this array span and density are similar to the San Jacinto Fault array that was deployed by Roux et al. (2016), which featured 1108 singlecomponent geophones distributed over a 650×700 m^{2} area. But instead of requiring 1108 individual stations, the hybrid array only requires 35 seismic stations, which dramatically reduces the hardware costs for such an operation. From this calculation it is apparent that there is a tradeoff between the extent of the array (which roughly scales as $\sqrt{\mathrm{\Delta}x}$) and the spatial density of the measurements (which roughly scales as $\mathrm{1}/\mathrm{\Delta}x$). The optimal spacing of the measurements must therefore be chosen appropriately for the application.
The complexity of deployment of a fibreoptic array is higher than for a temporary nodal array due to the need for trenching. However, since permanent seismometers are typically installed in shallow holes in the ground, the complexity and challenges of deploying a hybrid DASseismometer array is likely similar to those of a permanent seismometer array. This renders hybrid array designs with strategically deployed seismometers and DAS cable segments a costefficient alternative to permanent seismometer arrays, facilitating many applications including microseismicity monitoring, ambient noise tomography, and seismic source characterisation. One particular scenario that we wish to highlight here is that of rapid aftershock monitoring: the deployment of seismic stations following a large earthquake is time consuming and logistically challenging, and most often the temporary deployments are not completed in time to capture the earliest stage of the aftershock sequence. In contrast, fibreoptic cables can be deployed permanently in the 4shape or umbrella array configurations around permanent seismic stations in earthquakeprone regions, and, as soon as an earthquake of interest occurs, an interrogator unit may be readily connected to record the very earliest stages of the aftershock sequence.
This study considered the potential of fibreoptic distributed acoustic sensing (DAS) arrays for the purpose of seismic beamforming. 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 colocated 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 the nodal array produces an extremely wellresolved source location that is consistent with the true backazimuth 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 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. Fortunately, all of the above issues can be alleviated by converting the DAS strain rate measurements into particle velocities, taking the recordings of a nodal seismometer as a reference point and incrementally integrating the strain rates along a linear DAS segment. This approach warrants the design of “hybrid” arrays comprising both DAS and nodal seismometer arrays.
Python scripts that reproduce the results and figures in this paper are available at https://doi.org/10.6084/m9.figshare.12899288 (van den Ende and Ampuero, 2020).
MPAvdE conceptualised the study and performed the analyses. JPA supervised MPAvdE and contributed to methodology and interpretation. Both authors discussed and prepared the contents of the paper.
The authors declare that they have no conflict of interest.
This article is part of the special issue “Fibreoptic sensing in Earth sciences”. It is not associated with a conference.
The authors thank Lingsen Meng for sharing his MUSIC backprojection codes, from which the beamforming codes used in this study were derived. The nodal seismometer data were processed using ObsPy (Beyreuther et al., 2010), and generic data manipulations were performed with NumPy (Harris et al., 2020). Data visualisation was done using Matplotlib (Hunter, 2007). The authors are grateful for the constructive comments of the handling editor Philippe Jousset, two anonymous reviewers, and Eileen Martin.
This work was supported by the French government through the UCA^{JEDI} Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR15IDEX01.
This paper was edited by Philippe Jousset and reviewed by two anonymous referees.
AjoFranklin, J., Dou, S., Daley, T., Freifeld, B., Robertson, M., Ulrich, C., Wood, T., Eckblaw, I., Lindsey, N., Martin, E., and Wagner, A.: TimeLapse Surface Wave Monitoring of Permafrost Thaw Using Distributed Acoustic Sensing and a Permanent Automated Seismic Source, in: 2017 SEG International Exposition and Annual Meeting, Society of Exploration Geophysicists, SEG, 6093 pp., 2017. a
Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., and Wassermann, J.: ObsPy: A Python Toolbox for Seismology, Seismol. Res. Lett., 81, 530–533, https://doi.org/10.1785/gssrl.81.3.530, 2010. a
Capon, J., Greenfield, R., and Kolker, R.: Multidimensional MaximumLikelihood Processing of a Large Aperture Seismic Array, P. IEEE, 55, 192–211, https://doi.org/10.1109/PROC.1967.5439, 1967. a
Cole, S., Karrenbach, M., Kahn, D., Rich, J., Silver, K., and Langton, D.: Source Parameter Estimation from DAS Microseismic Data, in: 2018 SEG International Exposition and Annual Meeting, Society of Exploration Geophysicists, SEG, 5520 pp., 2018. a
Daley, T. M., Miller, D. E., Dodds, K., Cook, P., and Freifeld, B. M.: Field Testing of Modular Borehole Monitoring with Simultaneous Distributed Acoustic Sensing and Geophone Vertical Seismic Profiles at Citronelle, Alabama, Geophys. Prospect., 64, 1318–1334, https://doi.org/10.1111/13652478.12324, 2016. a
Dou, S., Lindsey, N., Wagner, A. M., Daley, T. M., Freifeld, B., Robertson, M., Peterson, J., Ulrich, C., Martin, E. R., and AjoFranklin, J. B.: Distributed Acoustic Sensing for Seismic Monitoring of The Near Surface: A TrafficNoise Interferometry Case Study, Sci. Rep.UK, 7, 11620, https://doi.org/10.1038/s41598017119864, 2017. a, b
Fang, G., Li, Y. E., Zhao, Y., and Martin, E. R.: Urban NearSurface Seismic Monitoring Using Distributed Acoustic Sensing, Geophys. Res. Lett., 47, e2019GL086115, https://doi.org/10.1029/2019GL086115, 2020. a, b, c, d
Feigl, K.: Brady's Geothermal Field DAS Earthquake Data [data set], University of Wisconsin, USA, https://doi.org/10.15121/1334285, 2016a. a
Feigl, K.: Brady's Geothermal Field Nodal Seismometer Earthquake Data [data set], University of Wisconsin, USA, https://doi.org/10.15121/1334284, 2016b. a
Feigl, K. L. and the PoroTomo Team: Overview and Preliminary Results from the PoroTomo Project at Brady Hot Springs, Nevada: Poroelastic Tomography by Adjoint Inverse Modeling of Data from Seismology, Geodesy, and Hydrology, in: 43rd Workshop on Geothermal Reservoir Engineering, Stanford University, Stanford, USA, 13–15 February 2017, 1715 pp., 2018. a, b, c, d, e
Goldstein, P. and Archuleta, R. J.: Array Analysis of Seismic Signals, Geophys. Res. Lett., 14, 13–16, https://doi.org/10.1029/GL014i001p00013, 1987. a
Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., GérardMarchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.: Array Programming with NumPy, Nature, 585, 357–362, https://doi.org/10.1038/s4158602026492, 2020. a
Hartog, A. H.: An Introduction to Distributed Optical Fibre Sensors, edn. 1, CRC Press, Boca Raton, Florida, USA, 472 pp., https://doi.org/10.1201/9781315119014, 2017. a
Hunter, J. D.: Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng., 9, 90–95, https://doi.org/10.1109/MCSE.2007.55, 2007. a
Hutchison, A. A. and Ghosh, A.: Ambient Tectonic Tremor in the San Jacinto Fault, near the Anza Gap, Detected by Multiple Mini Seismic Arrays, B. Seismol. Soc. Am., 107, 1985–1993, https://doi.org/10.1785/0120160385, 2017. a, b
Inbal, A., Ampuero, J. P., and Clayton, R. W.: Localized Seismic Deformation in the Upper Mantle Revealed by Dense Seismic Arrays, Science, 354, 88–92, https://doi.org/10.1126/science.aaf1370, 2016. a
Inza, L. A., Mars, J. I., Métaxian, J. P., O'Brien, G. S., and Macedo, O.: SeismoVolcano Source Localization with Triaxial BroadBand Seismic Array, Geophys. J. Int., 187, 371–384, https://doi.org/10.1111/j.1365246X.2011.05148.x, 2011. a
Jiang, C., Schmandt, B., Ward, K. M., Lin, F.C., and Worthington, L. L.: Upper Mantle Seismic Structure of Alaska From Rayleigh and S Wave Tomography, Geophys. Res. Lett., 45, 10350–10359, https://doi.org/10.1029/2018GL079406, 2018. a
Jolie, E., Moeck, I., and Faulds, J. E.: Quantitative StructuralGeological Exploration of FaultControlled Geothermal Systems – A Case Study from the Basin and RangeProvince, Nevada (USA), Geothermics, 54, 54–67, https://doi.org/10.1016/j.geothermics.2014.10.003, 2015. a
Jousset, P., Reinsch, T., Ryberg, T., Blanck, H., Clarke, A., Aghayev, R., Hersir, G. P., Henninges, J., Weber, M., and Krawczyk, C. M.: Dynamic Strain Determination Using FibreOptic Cables Allows Imaging of Seismological and Structural Features, Nat. Commun., 9, 2509, https://doi.org/10.1038/s4146701804860y, 2018. a, b
Kiser, E. and Ishii, M.: BackProjection Imaging of Earthquakes, Annu. Rev. Earth Pl. Sc., 45, 271–299, https://doi.org/10.1146/annurevearth063016015801, 2017. a, b
Krüger, F., Weber, M., Scherbaum, F., and Schlittenhardt, J.: Double Beam Analysis of Anomalies in the CoreMantle Boundary Region, Geophys. Res. Lett., 20, 1475–1478, https://doi.org/10.1029/93GL01311, 1993. a
Kuvshinov, B. N.: Interaction of Helically Wound FibreOptic Cables with Plane Seismic Waves, Geophys. Prospect., 64, 671–688, https://doi.org/10.1111/13652478.12303, 2016. a, b
Lellouch, A., Yuan, S., Spica, Z., Biondi, B., and Ellsworth, W. L.: Seismic Velocity Estimation Using Passive Downhole Distributed Acoustic Sensing Records: Examples From the San Andreas Fault Observatory at Depth, J. Geophys. Res.Sol. Ea., 124, 6931–6948, https://doi.org/10.1029/2019JB017533, 2019. a
Lin, F.C., Li, D., Clayton, R. W., and Hollis, D.: HighResolution 3D Shallow Crustal Structure in Long Beach, California: Application of Ambient Noise Tomography on a Dense Seismic Array, Geophysics, 78, 45–56, https://doi.org/10.1190/geo20120453.1, 2013. a
Lindsey, N. J., Martin, E. R., Dreger, D. S., Freifeld, B., Cole, S., James, S. R., Biondi, B. L., and AjoFranklin, J. B.: FiberOptic Network Observations of Earthquake Wavefields, Geophys. Res. Lett., 44, 11792–11799, https://doi.org/10.1002/2017GL075722, 2017. a, b
Lindsey, N. J., Dawe, T. C., and AjoFranklin, J. B.: Illuminating Seafloor Faults and Ocean Dynamics with Dark Fiber Distributed Acoustic Sensing, Science, 366, 1103–1107, https://doi.org/10.1126/science.aay5881, 2019. a, b, c, d, e, f
Lindsey, N. J., Rademacher, H., and AjoFranklin, J. B.: On the Broadband Instrument Response of FiberOptic DAS Arrays, J. Geophys. Res.Sol. Ea., 125, e2019JB018145, https://doi.org/10.1029/2019JB018145, 2020. a
Lior, I., Sladen, A., Rivet, D., Ampuero, J.P., Hello, Y., Becerril, C., Martins, H. F., Lamare, P., Jestin, C., Tsagkli, S., and Markou, C.: On the Detection Capabilities of Underwater DAS, J. Geophys. Res.Sol. Ea., 126, e2020JB020925, https://doi.org/10.1029/2020JB020925, 2021. a
Liu, H., Ma, J., Yan, W., Liu, W., Zhang, X., and Li, C.: Traffic Flow Detection Using Distributed Fiber Optic Acoustic Sensing, IEEE Access, 6, 68968–68980, https://doi.org/10.1109/ACCESS.2018.2868418, 2018. a
Martin, E. R., Lindsey, N., AjoFranklin, J., and Biondi, B.: Introduction to Interferometry of Fiber Optic Strain Measurements, EarthArXiv (preprint), https://doi.org/10.31223/osf.io/s2tjd, 2018. a, b, c, d, e
Meng, L., Inbal, A., and Ampuero, J.P.: A Window into the Complexity of the Dynamic Rupture of the 2011 Mw 9 TohokuOki Earthquake, Geophys. Res. Lett., 38, L00G07, https://doi.org/10.1029/2011GL048118, 2011. a, b, c
Meng, L., Allen, R. M., and Ampuero, J.P.: Application of Seismic Array Processing to Earthquake Early Warning, B. Seismol. Soc. Am., 104, 2553–2561, https://doi.org/10.1785/0120130277, 2014. a
Nakamichi, H., Yamanaka, Y., Terakawa, T., Horikawa, S., Okuda, T., and Yamazaki, F.: Continuous LongTerm Array Analysis of Seismic Records Observed during the 2011 Shinmoedake Eruption Activity of Kirishima Volcano, Southwest Japan, Earth Planets Space, 65, 7, https://doi.org/10.5047/eps.2013.03.002, 2013. a
Paitz, P., Edme, P., Gräff, D., Walter, F., Doetsch, J., Chalari, A., Schmelzbach, C., and Fichtner, A.: Empirical Investigations of the Instrument Response for Distributed Acoustic Sensing (DAS) across 17 Octaves, B. Seismol. Soc. Am., 111, 1–10, https://doi.org/10.1785/0120200185, 2020. a
Rieken, D. and Fuhrmann, D.: Generalizing MUSIC and MVDR for Multiple Noncoherent Arrays, IEEE T. Signal Proces., 52, 2396–2406, https://doi.org/10.1109/TSP.2004.832153, 2004. a, b, c
Ringdal, F. and Husebye, E. S.: Application of Arrays in the Detection, Location, and Identification of Seismic Events, B. Seismol. Soc. Am., 72, 201–224, 1982. a
Roux, P., Moreau, L., Lecointre, A., Hillers, G., Campillo, M., BenZion, Y., Zigone, D., and Vernon, F.: A Methodological Approach towards HighResolution Surface Wave Imaging of the San Jacinto Fault Zone Using AmbientNoise Recordings at a Spatially Dense Array, Geophys. J. Int., 206, 980–992, https://doi.org/10.1093/gji/ggw193, 2016. a
Schmidt, R.: Multiple Emitter Location and Signal Parameter Estimation, IEEE T. Antenn. Propag., 34, 276–280, https://doi.org/10.1109/TAP.1986.1143830, 1986. a, b
Scholz, C. H.: The Mechanics of Earthquakes and Faulting, edn. 3, Cambridge University Press, Cambridge, UK, https://doi.org/10.1017/9781316681473, 2019. a
Shearer, P. M.: Deep Earth Structure: Seismic Scattering in the Deep Earth, in: Treatise on Geophysics (edn. 2), edited by: Schubert, G., Elsevier, Oxford, UK, 759–787, https://doi.org/10.1016/B9780444538024.00018X, 2015. a
Singh, S., Capdeville, Y., and Igel, H.: Correcting Wavefield Gradients for the Effects of Local SmallScale Heterogeneities, Geophys. J. Int., 220, 996–1011, https://doi.org/10.1093/gji/ggz479, 2020. a, b, c
Sladen, A., Rivet, D., Ampuero, J. P., De Barros, L., Hello, Y., Calbris, G., and Lamare, P.: Distributed Sensing of Earthquakes and OceanSolid Earth Interactions on Seafloor Telecom Cables, Nat. Commun., 10, 5777, https://doi.org/10.1038/s4146701913793z, 2019. a, b, c
Stipčević, J., Kennett, B. L. N., and Tkalčić, H.: Simultaneous Use of Multiple Seismic Arrays, Geophys. J. Int., 209, 770–783, https://doi.org/10.1093/gji/ggx027, 2017. a
Thomson, D.: Spectrum Estimation and Harmonic Analysis, P. IEEE, 70, 1055–1096, https://doi.org/10.1109/PROC.1982.12433, 1982. a
van den Ende, M. and Ampuero, J.P.: Evaluating seismic beamforming capabilities of Distributed Acoustic Sensing arrays, figshare, https://doi.org/10.6084/m9.figshare.12899288 (last access: 18 April 2021), 2020.
Walter, F., Gräff, D., Lindner, F., Paitz, P., Köpfli, M., Chmiel, M., and Fichtner, A.: Distributed Acoustic Sensing of Microseismic Sources and Wave Propagation in Glaciated Terrain, Nat. Commun., 11, 2436, https://doi.org/10.1038/s41467020158246, 2020. a
Wang, H. F., Zeng, X., Miller, D. E., Fratta, D., Feigl, K. L., Thurber, C. H., and Mellors, R. J.: Ground Motion Response to an ML 4.3 Earthquake Using CoLocated Distributed Acoustic Sensing and Seismometer Arrays, Geophys. J. Int., 213, 2020–2036, https://doi.org/10.1093/gji/ggy102, 2018. a, b, c, d, e, f, g, h, i
Wiesmeyr, C., Litzenberger, M., Waser, M., Papp, A., Garn, H., Neunteufel, G., and Döller, H.: RealTime Train Tracking from Distributed Acoustic Sensing Data, Applied Sciences, 10, 448, https://doi.org/10.3390/app10020448, 2020. a
Zhan, Z.: Distributed Acoustic Sensing Turns FiberOptic Cables into Sensitive Seismic Antennas, Seismol. Res. Lett., 91, 1–15, https://doi.org/10.1785/0220190112, 2020. a, b, c, d
Zhu, T. and Stensrud, D. J.: Characterizing ThunderInduced Ground Motions Using FiberOptic Distributed Acoustic Sensing Array, J. Geophys. Res.Atmos., 124, 12810–12823, https://doi.org/10.1029/2019JD031453, 2019. a
Zhu, T., Shen, J., and Martin, E. R.: Sensing Earth and environment dynamics by telecommunication fiberoptic sensors: an urban experiment in Pennsylvania, USA, Solid Earth, 12, 219–235, https://doi.org/10.5194/se122192021, 2021. a
Zigone, D., BenZion, Y., Lehujeur, M., Campillo, M., Hillers, G., and Vernon, F. L.: Imaging Subsurface Structures in the San Jacinto Fault Zone with HighFrequency Noise Recorded by Dense Linear Arrays, Geophys. J. Int., 217, 879–893, https://doi.org/10.1093/gji/ggz069, 2019. a