Introduction to Interferometry of Fiber Optic Strain Measurements

E.R. Martin∗1,2, N.J. Lindsey, J.B. Ajo-Franklin, B. Biondi 1 Institute of Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA 2 Earth and Environmental Sciences Area, Lawrence Berkeley National Laboratory, Berkeley, CA, USA 3 Department of Earth and Planetary Sciences, University of California Berkeley, Berkeley, CA, USA 4 Geophysics Department, Stanford University, Stanford, CA, USA *Can be contacted at ermartin@sep.stanford.edu


INTRODUCTION
The cost of long-term seismic studies using dense arrays have historically been too expensive to maintain for many applications. By performing interferometry on ambient noise recorded by geophones, researchers can extract coherent signals mimicking surface wave Green's functions (Wapenaar 2010), (Lin 2008). This allows scientists to avoid the cost of a source crew, but long-term maintenance of a dense geophone array is too costly for some investigations. Underground fiber optic distributed acoustic sensing (DAS) arrays are an attractive alternative to geophone arrays due to ease of installation (potentially even using existing telecommunications infrastructure), durability, and the requirement for a single power source for thousands of sensors at meter-scale spacing. However, distributed acoustic sensing measures average axial strain rate along a subset of fiber, while geophones measure particle velocity in the vertical and sometimes both horizontal directions. This transition from a vector quantity (particle velocity) to a tensor quantity (strain rate) leads to significant changes to the angular sensitivity of the sensors, particularly for waves with particle motion primarily orthogonal to their direction of propagation.
Over the past few years, researchers have been able to extract coherent signals from ambient noise recorded by collinear channels at multiple surface DAS arrays installed in different ways (Dou 2017), (Lancelle 2016), (Martin 2015), (Martin 2016), (Zeng 2017A), (Zeng 2017B). Because these channels are already in a radial-radial alignment, the noise correlation functions are interpreted primarily as estimates of Rayleigh wave Green's functions. These results showed promise for recovering near-surface structure directly beneath fiber optic cables, but did not provide information about the subsurface between cables.
More recently, we have observed coherent signals in noise correlation functions between fiber optic lines at both a trenched array (Lindsey 2017A) and an array using fiber optics in existing telecommunications conduits (Martin & Biondi 2017), (Martin 2017), but we were previously unable to interpret these signals for use in subsurface imaging between fiber optic lines. In this paper, we develop the theory of ambient noise interferometry for fiber channels at arbitrary orientations and distances in the hope of greatly expanding the coverage of ambient noise interferometry for near surface imaging in include regions between fiber optic lines.
In addition to ambient noise interferometry studies, distributed acoustic sensing has become an increasingly popular option for microseismicity monitoring during energy production (Webster et al. 2013), (Karrenbach 2017), and is being investigated for earthquake monitoring at larger scales (Lindsey 2017B). Detection of small seismic events often relies on cross-correlations of recordings, so understanding the effect of directional sensitivity on both directly recorded signals and cross-correlations of signals could have implications for quantifying potential biases in statistics describing the distribution of these events estimated with a given array geometry.
In this paper we review the differences between theoretical geophone and DAS measurements of plane waves: both body waves and surface waves. We further analyze their cross-correlations for individual plane wave events, and for point sources. In fact, after doing this analysis, it was brought to our attention that a very different type of strain measurement device was studied nearly a century ago (Benioff 1935), so Figures 7 and8 of Benioff (1935) actually show up in our analysis. Using the theoretical point source responses and cross-correlations, we perform a simple thought experiment from (Wapenaar 2010) demonstrating how ambient noise interferometry works for random point sources recorded by both both horizontal geophones and DAS in radial-radial and transverse-transverse configurations. Building on the intuition from the thought experiment, we simulate noise correlation functions from cross-correlation of signals recorded by two parallel fiber lines and two orthogonal fiber lines. As in any ambient noise interferometry experiment, the source distribution matters, but biases caused by uneven source distributions for geophones and DAS are different.

DAS Measurement Process
Let a wave be described by particle velocityu(x, y, z, t) = (ux(x, y, z, t),uy(x, y, z, t),uz(x, y, z, t)) then the axial strain rate in the x-direction at a point is ∂ux ∂x | (x,y,z,t) . Let's say more generally we want to observe this wave field using a horizontal fiber measuring an axial strain rate in the (cos(θ), sin(θ), 0) direction. As seen in Appendix A, by applying tensor rotation matrices, the axial strain rate observed in the (cos(θ), sin(θ), 0) direction at that point iṡ where we use C θ and S θ as short notation for cos(θ) and sin(θ), respectively, andσ θ and all derivatives ofu are evaluated at (x, y, z, t). The measurement that would be detected by a fiber optic channel at that same θ orientation centered on (x, y, z) is the average axial strain rate over a gauge length of straight fiber: The gauge length g over which the interrogator unit averages the axial strain is in general different from the channel spacing. The channel spacing is the distance between the starting point of new segments. So, for instance if the channel spacing were 1 meter with a gauge length of 10 meters, the we would have a channel that observes the average axial strain along the fiber between 1 and 11 meters from the interrogator, the next channel would return the average axial strain between distances 2 and 12 meters, etc... The ability of the user (person acquiring data) to vary gauge length and channel spacing varies depending on manufacturer and age of the model.
In this paper we assume simplified geophone measurements,u at a given point with true amplitude response at all frequencies. We make two simplifying assumptions on DAS systems to keep our results more general across different DAS implementations: (1) It is possible some DAS implementations may use a different weighted averaging over the gauge length, and it has been shown that the signal may vary depending on pulse shape and gauge length (Bona 2017). We simply pretend that the axial strain rate along all points of the gauge are weighted equally. (2) We assume that the average along the gauge length of the axial strain rate of the medium is observed. In reality there is a scaling coefficient provided by each manufacturer, and the Lamé parameters of the fiber and jacket as well as the friction between the fiber, jacket and ground affect these measurements. More details on this effect can be found in (Kuvshinov 2016). To first order, we expect the math in this paper to explain the most significant changes observed when switching from geophone to DAS measurements, but more accurate modeling of any given data set could account for these implementation-specific details provided by the DAS interrogator unit manufacturer.

SENSITIVITY OF DAS TO FAR-FIELD SOURCES
We are interested in understanding from which directions far-field sources are emphasized when recorded by a horizontal sensor array, both directly recorded and after cross-correlation. Far away sources are generally well represented by plane waves, so we study this approximation. Depending on whether the wave detected is a P, S, Rayleigh or Love wave, different source angles are emphasized. In Appendix B following the notation of (Pujol 2003) we derive the sensitivity of a horizontal fiber oriented in the direction (C θ , S θ , 0) to surface waves (Rayleigh and Love) propagating in the direction (C φ , S φ , 0), as drawn in Figure 1. Similarly, we derive these results for body waves (P and S) propagating in the direction (C φ 1 C φ 2 , S φ 1 C φ 2 , S φ 2 ). We assume these are monochromatic plane waves traveling in a half-space at phase velocity c, wavenumber k and frequency ω = kc. For each type of wave, we use the following parameters to describe their particle displacement and velocity in Table 1.
• Rayleigh waves: We assume a homogeneous half-space. Let α and β be the velocities of the irrotational and solenoidal parts of the solution to the elastic wave equation (so uα + u β = u, × uα = 0 and · u β = 0). Then define γα = 1 − c 2 α 2 and γ β = 1 − c 2 β 2 . A and B are amplitude factors, whose ratio defines the ellipticity of the Rayleigh wave. ϕ θ (x,y,z) Figure 1. We study surface waves (indicated by the red and blue bars) that propagate in the direction (C φ , S φ , 0) as they are observed by horizontal sensors (indicated by the green line) at (x, y, z) oriented in the (C θ , S θ , 0) direction.
• Love waves: We assume a homogeneous half-space with group velocity β2 underneath a top layer of thickness H with group velocity β1. In this paper, we consider only surface seismic, so we just use the displacement in the top H meters of the subsurface. Let A and B be amplitude terms of the depth factor in the top layer (independent of the A and B used to describe Rayleigh waves), and η1 = c 2 β 2 1 − 1.
• P waves: We assume a homogeneous half space. Let A be an amplitude factor (independent of amplitudes of other wave types). All other notation is common to all plane waves.
• S waves: We assume a homogeneous half space. We split the S-wave motion into two parts: SH which is just the horizontal particle motion, and SV which occurs in the plane spanned by the z-axis and the direction of propagation. Let A be an amplitude factor (independent of amplitudes of other wave types). All other notation is common to all plane waves.
As summarized in Table 2, in Appendix B we calculate the geophone particle velocity response to each type of wave in the (C θ , S θ , 0) direction, which we denote byu θ (x, y, z, t). We also calculate the point-wise axial strain rate in the same directionσ θ (x, y, z, t) and the expected DAS signal denoted byσ θ,g (x, y, z, t). Note that in the long wavelength (k → 0) limit, a DAS channel with a finite gauge length is predicted to yield the same signal as a point-wise axial strain rate measurement.
Given a predicted measurement,u θ ,σ θ ,σ θ,g , to a particular wave, define its sensitivity as that measurement disregarding any terms only dependent on depth and any oscillatory terms. In general, foru θ ,σ θ andσ θ,g when g is short enough relative to the wavelength, the sensitivity is similar for P-waves, SV-waves and Rayleigh waves (a two-lobed sensitivity pattern). Love waves and SH-waves are similar to each other, but,u θ has a π/2 rotated two-lobed sensitivity pattern, whileσ θ , andσ θ,g (for small enough g) have four-lobed sensitivity patterns.
First, let's look at the surface waves sensitivities. We plot the sensitivity to Rayleigh and Love waves coming from different angles in Figure 2 for a geophone, a point-wise axial strain rate measurement,and a DAS channel with gauge length 10 meters (a typical setting for seismic acquisition). While the geophone is most sensitive to Rayleigh waves with φ − θ = 0, π, it is equally sensitive to Love waves with φ − θ = ±π/2 in the sense that these sensitivity patterns are just rotated versions of each other. Like the geophone, the point-wise axial strain rate measurement is also most sensitive to Rayleigh waves such that φ − θ = 0, π, but the sensitivity is more concentrated around these peak angles (i.e. is scaled by an extra kC (φ−θ) compared to the geophone). Unlike the geophone, the point-wise axial strain rate has a four-lobed sensitivity to Love waves with peak sensitivity to Love waves coming in at φ − θ = ±π/4, ±3π/4. Further, at its peak sensitivity angles to Love waves, the point-wise axial strain rate measurment only detects half the amplitude that would be detected by the same theoretical sensor to Rayleigh waves at φ − θ = 0, π. Thus, the switch to strain means Love wave sensitivity is weaker and distributed over more angles.
For wavelengths a few times longer than g the sensitivity patterns ofσ θ well-approximate the sensitivity ofσ θ,g . As seen in Figure 2, for wavelengths just slightly bigger than g the sensitivity patterns become more flattened out near peak-sensitivity angles, then for wavelengths shorter than g, additional sensitivity lobes form, consistent with the description for P-waves in (Dean 2017). This trend can be seen in more detail in Figure 3, which shows that for larger gauge lengths, these stranger behaviors (deviating from the 2-lobe and 4-lobe trends) can occur at even longer wavelengths.
The response to body waves is slightly more complicated depending on how vertically the wave is propagating, indicated by φ2. We plot the sensitivity for φ2 = 3π/8, π/4, and π/8 in Figures 4, 5, and 6, respectively. In general, the P-wave and SV-wave sensitivity patterns Table 1. Plane wave particle displacements and velocities as described in Appendix B. The oscillatory factor o RL in the surface waves o R L = e ik(ct−xC φ −yS φ ) , and for body waves the oscillatory factor is o Wave type quantity value quantity value Table 2. Plane wave particle velocity,u θ , point-wise axial strain rate,σ θ , and strain rate averaged over a gauge length,σ θ in the (C θ , S θ , 0) direction. The oscillatory factor o RL in the surface waves o RL = e ik(ct−xC φ −yS φ ) , and the oscillatory factor o P S in the body waves is Wave type quantity value Ao P S Figure 2. Each polar plot shows the amplitude response of a measurement to planar surface waves of varying azimuth and wavelength. The radius of each line represents the sensitivity of DAS with a 2 meter (left), 5 meter (2nd column), 10 meter (3rd column) and 20 meter (right) gauge length to a range of wavelengths in a 400 m/s velocity material for both Rayleigh (green) and Love (red) plane waves coming from each angle . The plots represent sensitivities for 9 Hz (top), 19 Hz (2nd row), 29 Hz (3rd row), and 39 Hz (bottom).
are more like Rayleigh wave sensitivity (two lobed), while the SH-wave sensitivity is more like Love wave sensitivity (four lobed). When φ2 is closer to π/2 that means the wave is propagating almost vertically, so most particle displacement is closer to vertical and there is less P-wave sensitivity for all three horizontal measurements compared to P-waves traveling at smaller φ2 values. While SV geophone response is maximized for nearly-vertically propagating waves (φ2 = π/2) and minimized for nearly horizontally propagating waves (φ2 = 0), the point-wise strain rate and DAS responses are minimized for both horizontal and vertically propagating waves and maximized when there is as much vertical as their is horizontal propagation (φ2 = π/4). On the other hand, the geophone response to SH waves is independent of how much energy is propagating vertically, but the point-wise strain rate and DAS measurements are maximized when an SH-wave is propagating horizontally (φ2 = 0).

SENSITIVITY OF DAS CROSS-CORRELATIONS TO PLANE WAVE SOURCES
When attempting to use weak signals, as in ambient noise interferometry or some microseismicity detection methods, we often crosscorrelate two time series signals, s(xA, yA, zA, t) and s(xB, yB, zB, t), recorded concurrently at different sensors to bring their joint signal, C(τ ) = 1 2T T −T s(xA, yA, zA, t)s * (xB, yB, zB, t + τ )dt, above the noise level. In particular, this is used to more accurately detect and locate microseismic events using an array, and is also used to extract surface wave Green's function approximations from ambient noise recorded by surface arrays. While the extensive sensor coverage of DAS is clearly an advantage in microseismicity detection, could the use of certain channel geometries relative to events lead to different biases in the estimated statistical distribution of events? When performing ambient noise interferometry in the presence of an ideal noise field, does the extracted signal actually approximate the same arrival times as the true Green's functions?
The answers to these questions begin by studying the differences between cross-correlations of pairs of geophones, point-wise axial strain rate measurements, and DAS channels responding to plane waves. When using 3C geophones, we can rotate any pair of sensors into radial and transverse horizontal components that clearly yield Rayleigh and Love wave components (Lin 2008), but as we only observe one component of the strain tensor with fiber, we cannot do the same with a horizontal fiber array. Thus, we calculate cross-correlations between Figure 3. The radius of each line represents the sensitivity of DAS with a 2 meter (left), 5 meter (2nd column), 10 meter (3rd column) and 20 meter (right) gauge length to a range of wavelengths in a 400 m/s velocity material for both Rayleigh (green) and Love (red) plane waves coming from each angle . The plots represent sensitivities for 9 Hz (top), 19 Hz (2nd row), 29 Hz (3rd row), and 39 Hz (bottom). a sensor at (x1, y1, z1) oriented in the (C θ 1 , S θ 1 , 0) direction and a second sensor at (x2, y2, z2) oriented in the (C θ 2 , S θ 2 , 0) direction. These cross-correlations are summarized in Table 3, as derived in Appendix B.
Say we have two sensors at two surface locations x1 and x2 such that x 1 −x 2 x 1 −x 2 = (1, 0, 0). A plane surface wave coming in from angle φ = 0 or π would hit x1 and x2 at what appears to be the true velocity, but a plane wave coming in from angle φ = ±π/2 would arrive at x1 and x2 at the same time at an infinitely fast velocity, and any angle in between will have some fast biased apparent velocity, so ideally, we want whichever wave type we are trying to detect, to be emphasized close to φ = 0, π.
Ignoring the oscillatory terms and exponentially decaying depth-dependent surface wave terms, we have plotted cross-correlation sensitivities for x1 and x2 if they were geophones, point-wise strain rates, and DAS responding to surface waves in Figures 7 and 8, which show their radial-radial (θ1 = θ2 = 0) and transverse-transverse (θ1 = θ2 = π/2) cross-correlation sensitivities respectively. These are the cross-correlations we use for extracting Rayleigh waves and Love waves, respectively from geophones, and as expected, the radial-radial cross-correlation in Figure 7 is very sensitive to Rayleigh waves at φ = 0, π, and the transverse-transverse cross-correlation in Figure 8 is very sensitive to Love waves at φ = 0, π. However, when we look at the point-wise strain rate diagrams (the same as DAS when the gauge length is chosen appropriately), the radial-radial correlations in Figure 7 are very sensitive to Rayleigh waves at φ = 0, π, but the transverse-transverse correlations in Figure 8 are not very sensitive to any waves at φ = 0, π, and in fact are very sensitive to Rayleigh waves that would yield a very fast (even infinite) velocity.
So how to we detect any Love waves from cross-correlations? The transverse-transverse cross-correlations are sensitive to Love waves coming from φ = ±π/4, ±3π/4, just not as sensitive as they are to Rayleigh waves near φ = ±π/2 which would appear to have an extremely fast velocity. If we have reason to believe that Love waves dominate the wavefield, then transverse-transverse cross-correlations will result in an apparent velocity that is c/ cos(Φ), where c is the true Love wave velocity at that frequency, and Φ is the peak Love wave cross-correlation sensitivity angle within the frequency band of interest (Φ = π/4 for longer wavelengths relative to the gauge length).

THOUGHT EXPERIMENT DEMONSTRATING AMBIENT NOISE INTERFEROMETRY TRENDS
A simple thought experiment to understand why hydrophone (pressure) and vertical geophone ambient noise cross-correlations averaged over many sources yield a signal with its peak around the Green's function arrival time of one of the receivers responding to a virtual source at the other can be found in (Wapenaar 2010). Here, I recreate that thought experiment at the same scale, but study the cross-correlations of the particle velocity and point-wise axial strain response to random sources that travel either like Rayleigh waves or Love Waves (although at a faster velocity than surface waves and without dispersion). Further, I test a range of geometries encountered in the experiments examined in this thesis: a virtual source that is a fiber optic channel cross-correlated with a set of receivers on a parallel cable, and a fiber optic virtual source cross-correlated with a set of receivers on an orthogonal cable.
Here, I assume a simple model of surface wave point sources: a Ricker wavelet set off at a point that either has particle motion in the direction of propagation (Rayleigh wave) or horizontal and orthogonal to the direction of propagation (Love wave). This is described in detail in Appendix C along with the derivations of the expected response of the xand y-components of the particle velocity and axial strain rate to these point sources, which represent idealized horizontal geophones and DAS cables. These responses are summarized in Table 4 4.1 Simple Case: Radial-Radial Cross-correlations Prior work has shown in practice that it is possible to retrieve reasonable geology from surface wave inversion using cross-correlations of data recorded on collinear DAS channels (@warning Citation 'Dou2017' on page 7 undefined). This section aims to better understand why this works. Following the thought experiment of (Wapenaar 2010), imagine two receivers at x1 = (−600, 0, 0) and x2 = (600, 0, 0) surrounded by 5000 point sources randomly distributed on an annulus with inner and outer radii of 2000 and 3000 m (uniformly distributed azimuth in [0, 2π] and uniformly distributed radius in [2000,3000]). While the original experiment studied a scalar response (pressure), in this section I am interested in the radial-radial cross-correlations of geophones or DAS channels oriented following the geometries shown in Figure 9, so θ1 = θ2 = π/2. The recordings of these sources (every eighth source for visualization) on both sensors acting as geophones (particle velocity) and as DAS channels (point-wise strain rate) is shown in Figure 10. Both the geophones and the fiber optics emphasize sources with φSrc close to Figure 5. The radius of each line represents the sensitivity of geophones (left), point-wise strain measurements (middle) and DAS with a 10 meter gauge length (right) to a different wavelengths in a 400 m/s velocity material for both P (orange), SH (blue) and SV (black) plane waves coming from each horizontal angle φ 1 − θ and vertical angle φ 2 = π/4. The plots represent sensitivities for 9 Hz (top), 19 Hz (2nd row), 29 Hz (3rd row), and 39 Hz (bottom). 0 and π (i.e. sources that aren't observed with any apparent velocity bias), and in fact, the relative emphasis of these sources by the fiber channels is stronger than the relative emphasis observed by the geophone. This is also seen in the cross-correlations for each source. As in (@warning Citation 'Wapenaar2010A' on page 7 undefined), even though the initial recordings had random time lags due to their radius being chosen over a range between 2000 and 3000 meters, their relative arrival times at x1 and x2 are consistent, so the cross-correlations show a clear trend. When we average the source-wise cross-correlations we get a clear signal with a peak at ±0.6 seconds for both the geophone and DAS experiments. Because x1 and x2 are 1200 meters apart in a 2000 m/s medium, that is the arrival time we would expect for a source emitted from one receiver's location and recorded at the other receiver.
While the average of single-source cross-correlations can start to give some intuition about why ambient noise interferometry works, in reality, each window of noise contains the responses to many sources. If a finite number, N , of point sources go off during a particular time window, then the cross-correlation of two traces (denoted by u) recording these sources would include both single-source cross-correlations and cross-terms between different sources. This is apparent when the cross-correlations are written as a product of Green's functions, G, and source functions f : Only the first term (single source cross-correlations) are explained by Figure 10. Thus, (@warning Citation 'Wapenaar2010A' on page 8 Figure 6. The radius of each line represents the sensitivity of geophones (left), point-wise strain measurements (middle) and DAS with a 10 meter gauge length (right) to a different wavelengths in a 400 m/s velocity material for both P (orange), SH (blue) and SV (black) plane waves coming from each horizontal angle φ 1 − θ and vertical angle φ 2 = π/8. The plots represent sensitivities for 9 Hz (top), 19 Hz (2nd row), 29 Hz (3rd row), and 39 Hz (bottom). ER undefined) performed a slightly more realistic test case to ensure that these cross-source-terms did not add up to make coherent changes in the extracted velocity. The test is to take a long time period and record a single trace for each sensor while letting a number of random sources go off. These are recorded at both x1 and x2 (and I repeat this for both a particle velocity and a point-wise strain rate). Then a single cross-correlation is done between the long window of recording at x1 with the long window of recording at x2. The resulting singlewindow cross-correlations for both fiber and geophones responding to 1000 random Rayleigh wave point sources throughout a 40,000 second window are shown in Figure 11. Also shown are the cross-correlations of the same experiment repeated with 1000 random Rayleigh wave point sources and 1000 random Love wave point sources of the same amplitude. The signal extracted when both Rayleigh and Love waves are present is a bit noisier than when just Rayleigh waves are present, but the peak at ±0.6 seconds is easy to pick for both the fiber and geophone responses.

Transverse-Transverse Cross-correlations
While the radial-radial cross-correlations yield a clear peak at the correct ±0.6 second time lags for both geophones and DAS. Radialradial cross-correlations were already observed to be a simpler situation in that they both respond strongly to Rayleigh wave sources at φ = 0, π (the azimuths corresponding to true velocity sources). The same does not hold true for transverse-transverse cross-correlations in the setup pictured in Figure 9. To better understand transverse-transverse cross-correlations, I repeated the exercise from (@warning Citation 'Wapenaar2010A' on page 9 undefined) but with just Love wave sources, and the results can be seen in Figure 12. While the geophones emphasize Love wave sources at φSrc = 0, π, the fiber channels emphasize Love wave sources at φSrc = −π/4, π/4, 3π/4, 5π/4. This relative difference is even more pronounced in the cross-correlations of the records for each of these sources. As has already been confirmed in practice (@warning Citation 'Lin2008' on page 9 undefined), the average geophone cross-correlation has a strong peak at ±0.6 seconds which is the correct arrival time. The average fiber cross-correlation is much more spread out, which may be in part because wavelets don't cancel as cleanly away from peak sensitivity angles, but also due to the amplitude difference predicted in Table 8. As predicted, the peak of the average fiber cross-correlation is between ±0.4 and ±0.5 seconds (since 0.42 seconds = 0.6 seconds / √ 2). I also repeated the exercise from (@warning Citation 'Wapenaar2010A' on page 9 undefined) recording 1000 random Love wave point Table 3. Cross-correlations of measurements at (x 1 , y 1 , z 1 ) in the (C θ 1 , S θ 1 , 0) direction with measurements at (x 2 , y 2 , z 2 ) in the (C θ 2 , S θ 2 , 0) direction. To keep the notation short, we introduce notation for common oscillatory terms: Wave type quantity value sources spanning a 40000 second-long trace for each location in both particle velocity and strain rate measurements, then doing a single cross-correlation for each type of measurement, pictured in Figure 13. Again, because there is only a single window with cross-terms due to many sources rather than averaging over multiple windows that each have a single source, the resulting cross-correlations are noisier than in Figure 12. The geophone cross-correlation again yields correct peaks at ±0.6 seconds, but again, the fiber signal is spread out with a peak somewhere less than ±0.5 seconds. I repeated this exercise in the presence of 1000 Love wave sources and 1000 Rayleigh wave sources of equal amplitude, also shown in Figure 13. There might be some hope of recovering a Love wave signal (that can be corrected) from this receiver geometry when only in the presence of Love wave sources, but it appears that the sensitivity to very apparently fast Rayleigh wave sources near φSrc = π/2, 3π/2 dominates the signal too much to recover a Love wave signal. Table 4. Horizontal particle velocity and point-wise axial strain rate responses at x to a simplified model of Rayleigh and Love wave point sources at xs that emit a Ricker wavelet at frequency f . To keep notation short, let Figure 7. The radius of each line represents the sensitivity of radial-radial (θ 1 = θ 2 = 0) cross-correlations geophones (left), point-wise strain measurements (middle) and DAS with a 10 meter gauge length (right) to a range of wavelengths in a 400 m/s velocity material for both Rayleigh (green) and Love (red) plane waves coming from each angle. The plots represent sensitivities for 19 Hz (top), 29 Hz (2nd row), and 39 Hz (bottom). 9 Hz was not pictured because the responses are too small to see.

SIMULATED AMBIENT NOISE INTERFEROMETRY
The exercises in Section 4 give some intuition behind the signals extracted by ambient noise interferometry ofσ θ in two common setups: radial-radial and transverse-transverse, but even in simple fiber arrays made up of parallel or orthogonal lines, we wish to understand the signals extracted from a wider range of ray paths. Although truly transverse-transverse channels may not yield a detectable Love wave signal when Rayleigh wave sources are also present, it is natural to wonder whether such a signal may be possible from two parallel channels offset by angles close to π/4 (so x A − x B lines up with the direction of peak Love wave source sensitivity. Further, we have observed clear signals extracted from orthogonal and parallel channel pairs which are neither exactly radial-radial nor exactly transverse-transverse (Lindsey 2017A), (Martin & Biondi 2017), (Martin 2017), so we wish to understand the meaning of these signals so that they can be properly used in tomography.

Signals Extracted Between Two Parallel Fiber Lines
Although the transverse-transverse fiber geometry does not appear to yield a clear signal, it seems reasonable that cross-correlations of parallel channels such that x1 − x2 is in a direction close to (C π/4 , S π/4 , 0) would have a Love wave signal since this is the direction that is most sensitive to Love waves. However, these sensors also have significant sensitivity to Rayleigh waves, so it seems plausible that the cross-correlations might include both a Rayleigh wave and Love wave signal. I put this to a test similar to the single cross-correlation over a long recording window in the previous section. I again use a 30 Hz Ricker wavelet (same expressions as in previous section) with 10000 random Love wave and 10000 random Rayleigh wave sources uniformly distributed over an annulus with inner radius 2000 m and outer radius 3000 m. These are recorded during a 400,000 second record. The Rayleigh wave velocity is again 2000 m/s, but the Love wave is 20% faster: 2400 m/s. This way it will be possible to distinguish the Rayleigh and Love wave signals from each other easily, and if the Love wave signal is apparently fast it will just lead to better separation from any Rayleigh wave signal. As pictured in Figure 14, the virtual source is a point-wise axial strain measurement oriented in the (0, 1, 0) direction and sits at (425, 425, 0), and there are 7 other point-wise axial strain receivers on a parallel line, each oriented in the Figure 8. The radius of each line represents the sensitivity of transverse-transverse (θ 1 = θ 2 = π/2) cross-correlations geophones (left), point-wise strain measurements (middle) and DAS with a 10 meter gauge length (right) to a range of wavelengths in a 400 m/s velocity material for both Rayleigh (green) and Love (red) plane waves coming from each angle. The plots represent sensitivities for 19 Hz (top), 29 Hz (2nd row), and 39 Hz (bottom). 9 Hz was not pictured because the responses are too small to see.
The cross-correlations of these long records for each receiver against the virtual source are shown in Figure 14. For the receivers within a 15 o offset of the transverse-transverse orientation (i.e. the two bluest receivers) there is a lot of energy from the very apparently fast velocity Rayleigh wave sources from φSrc = ±π/2, so there is no clear arrival at the true Love and Rayleigh wave arrival times. Moving down to the third receiver, about a 25 o offset, there is still quite a bit of this fast energy, but there is also a strong peak at the true Love wave velocity and  a smaller peak (although not nearly as clear) at the true Rayleigh wave velocity. All of the receivers farther down have both clear Rayleigh and Love arrivals and less of the early arrival energy. The Rayleigh wave arrivals get increasingly strong moving down the line because the orientation between the receivers and the virtual source gets closer to the two channels being collinear (like a radial-radial setup). Figure 11. A single long radial-radial cross-correlation of synthetic geophone (black) data recorded in the presence of many Rayleigh wave point sources (Left) yields a coherent signal the correct arrival time of ±0.6 seconds, and the same holds true for the process repeated with synthetic fiber (red) data. (Right) Even when Love wave sources are present at equal amplitudes to the Rayleigh wave sources, the correct arrival time can still be picked clearly. Figure 12. Random synthetic point sources emitting Love waves were recorded via particle velocity and strain rate at x 1 and x 2 in the (0, 1, 0) direction. (Top left) The geophones respond strongly to Love waves coming from the φ Src = 0, π directions. (Top right) The fiber channels respond strongly to Love waves coming from the φ Src = −π/4, π/4, 3π/4, 5π/4 directions. (Bottom left) For each source, the geophone cross-correlation is plotted in black and the fiber cross-correlation is plotted in red. (Bottom right) The average of these source-wise cross-correlations is plotted for the geophones in black and the fiber in red. The peak geophone signal is around ±0.6 seconds, and the peak fiber signal is more spread in time with a peak in the ±0.4 to ±0.5 second range.

Virtual Source Perpendicular to Receiver Line
At corners of arrays where one fiber line is orthogonal to another, it would be ideal to be able to use the ray paths between those two lines. To better understand this situation, I test a virtual source on one line at (425, 425, 0) oriented in the (1, 0, 0) direction, and seven equally spaced receivers oriented in the (0, 1, 0) direction between (−425, −800, 0) and (−425, 425, 0), as pictured in Figure 15. The same velocities and same configuration of Rayleigh and Love wave point sources were used as in the previous section.  The cross-correlation results are seen in Figure 15. No clear signal is visible in the cross-correlation with the top receiver that is directly across from the virtual source, but for all other offsets, both the Rayleigh and Love wave signals are clearly visible. While the relative amplitudes between the Rayleigh and Love wave signal peaks varied significantly with offset in the parallel lines setup, the relative amplitudes of the Rayleigh and Love wave peaks stay consistent over offsets in this orthogonal lines setup. Overall, it appears that the orthogonal lines setup is actually more reliable for simultaneously extracting both Rayleigh and Love wave signals than the parallel lines setup, but the parallel lines setup can be used for a limited range of offsets.

ACKNOWLEDGMENTS
We would like to acknowledge useful discussions about the content of this paper with Jason Chang (Stanford), Bob Clapp (Stanford), Nori Nakata ( However, strain rate is a tensor quantity, so any strain rate measurements observed in a rotated coordinate system actually behave differently than particle velocity. Here, we show how that difference crops up for axial strain rate measurements. As a reminder, using the notation that = (∂x, ∂y, ∂z) T (so is a column vector), and writing out the definition of our strain rate tensor, we can see how it is described in terms of andu: Now let's look atΣ , the strain tensor as observed in that same rotated coordinate system asu . First, note that if we take derivatives in a rotated (denoted by ') coordinate system, = (∂ x , ∂ y , ∂ z ) T = R . We start out with the definition of the strain rate in the rotated coordinate system, similar to Equation A.2, then plug in our rotated expressions for andu : so we see that the strain rate tensor in the rotated coordinate system is just the strain rate tensor in the original coordinate system with a rotation matrix and its transpose applied to either side. Expanding out Equation A.3, we see how each term of the rotated strain rate tensor can be expressed from the matrix multiplication: The particular term of interest isσ x x because it is the point-wise axial strain rate that would be detected by a fiber oriented in the direction (C θ , S θ , 0) if the fiber interrogator unit could return an infinitely small gauge length measurement. Expanding out the top-left term of Equation A.4, we see that the axial strain rate can be expressed in terms of strains observed the original x − y coordinate system as: As shorthand for this, we use the notationσ θ :=σ x x in the rest of this paper.

APPENDIX B: DERIVING VELOCITY AND STRAIN MEASUREMENTS OF PLANE WAVES AND INTERFEROMETRY OF THOSE MEASUREMENTS
In this appendix, we investigate the response of several types of idealized sensors to plane waves of multiple types: Rayleigh, Love, P, SH and SV. We are interested in the response of an idealized geophone (that measures exact particle velocity response at all frequencies), a point-wise axial strain rate measurement in some horizontal direction, and an average axial strain rate measurement over a line segment of length g (the gauge length). For each type of wave, we start with a displacement expression from (Pujol 2003), although we generalize our starting displacement expressions so that the plane wave propagation direction is not restricted to the x − z plane. Thus, we can assume any fixed horizontal sensor orientation, θ, then consider the response of that sensor to signals approaching the sensor from different angles (φ for surface waves, or φ1, φ2 for body waves). Starting with those displacements, we derive the following responses for each type of wave: (i) We take a time derivative for 3D particle velocitiesu = (ux,uy,uz), summarized in Table 1. (ii) By applying a single rotation matrix tou we can calculate the particle velocityu θ as detected by a geophone in the in the (cos(θ), sin(θ), 0) direction. These results are summarized in Table 2 (iii) We calculate horizontal spatial derivatives ofu then plug them into Equation A.5 to calculate the pointwise axial strain rate in the (cos(θ), sin(θ), 0) direction,σ θ . These results are summarized in Table 2.
(iv) We averageσ θ over a horizontal line segment of length g, the gauge length, yieldingσ θ,g , which mimics a DAS signal in the same location in the direction (cos(θ), sin(θ), 0). These results are summarized in Table 2.
(v) We verify that lim k→0σ θ,ġ σ θ = 1, i.e. at long wavelengths, the average strain rate over a fixed gauge length acts like a point-wise axial strain rate measurement in the same direction.

B1 Rayleigh waves
(i) We take a time derivative for 3D particle velocitiesu = (ux,uy,uz) By applying a single rotation matrix tou we can calculate the particle velocityu θ as detected by a geophone in the in the (cos(θ), sin(θ), 0) direction.
Now we can plug these in for the point-wise axial strain in the θ direction: (iv) We averageσ θ over a horizontal line segment of length g, the gauge length, yieldingσ θ,g , which mimics a DAS signal in the same location in the direction (cos(θ), sin(θ), 0).
(v) We verify that lim k→0σ θ,ġ σ θ = 1, i.e. at long wavelengths, the average strain rate over a fixed gauge length acts like a point-wise axial strain rate measurement in the same direction. This requires the use of L'Hospital's rule: (vi) We calculate the cross-correlation of two horizontal particle velocity measurements,u θ 1 (x1, y1, z) andu θ 2 (x2, y2, z), simultaneously recording the same plane wave.

2T
T −Tσ = 1, verifying that the behavior of interferometry of DAS measurements is the same as interferometry of point-wise axial strain rate measurements at long wavelengths. This requires two applications of L'Hospital's rule for indeterminate forms: B2 Love waves (i) We take a time derivative for 3D particle velocitiesu = (ux,uy,uz) (ii) By applying a single rotation matrix tou we can calculate the particle velocityu θ as detected by a geophone in the in the (cos(θ), sin(θ), 0) direction.
Now we can plug these in for the point-wise axial strain in the θ direction: (iv) We averageσ θ over a horizontal line segment of length g, the gauge length, yieldingσ θ,g , which mimics a DAS signal in the same location in the direction (cos(θ), sin(θ), 0).
(v) We verify that lim k→0σ θ,ġ σ θ = 1, i.e. at long wavelengths, the average strain rate over a fixed gauge length acts like a point-wise axial strain rate measurement in the same direction.

B3 P waves
(i) We take a time derivative for 3D particle velocitiesu = (ux,uy,uz): By applying a single rotation matrix tou we can calculate the particle velocityu θ as detected by a geophone in the in the (cos(θ), sin(θ), 0) direction.
(iii) We calculate horizontal spatial derivatives ofu then plug them into Equation A.5 to calculate the pointwise axial strain rate in the (cos(θ), sin(θ), 0) direction,σ θ . First, let's calculate the point-wise axial strain of a horizontal fiber oriented in the (C θ , S θ , 0) direction at the point (x, y.z). The first steps are to calculate our spatial derivatives of the particle velocity, keeping in mind thatux Now we can plug these in for the point-wise axial strain in the θ direction: (iv) We averageσ θ over a horizontal line segment of length g, the gauge length, yieldingσ θ,g , which mimics a DAS signal in the same location in the direction (cos(θ), sin(θ), 0).
(v) We verify that lim k→0σ θ,ġ σ θ = 1, i.e. at long wavelengths, the average strain rate over a fixed gauge length acts like a point-wise axial strain rate measurement in the same direction.

B4 SV waves
(i) We take a time derivative for 3D particle velocitiesu = (ux,uy,uz) (ii) By applying a single rotation matrix tou we can calculate the particle velocityu θ as detected by a geophone in the in the (cos(θ), sin(θ), 0) direction.
We calculate horizontal spatial derivatives ofu then plug them into Equation A.5 to calculate the pointwise axial strain rate in the (cos(θ), sin(θ), 0) direction,σ θ . First, let's calculate the point-wise axial strain of a horizontal fiber oriented in the (C θ , S θ , 0) direction at the point (x, y.z). The first steps are to calculate our spatial derivatives of the particle velocity, keeping in mind thatux Now we can plug these in for the point-wise axial strain in the θ direction: (iv) We averageσ θ over a horizontal line segment of length g, the gauge length, yieldingσ θ,g , which mimics a DAS signal in the same location in the direction (cos(θ), sin(θ), 0). σ θ,g (x, y, z, t) = 1 g (v) We verify that lim k→0σ θ,ġ σ θ = 1, i.e. at long wavelengths, the average strain rate over a fixed gauge length acts like a point-wise axial strain rate measurement in the same direction.

B5 SH waves
(i) We take a time derivative for 3D particle velocitiesu = (ux,uy,uz) By applying a single rotation matrix tou we can calculate the particle velocityu θ as detected by a geophone in the in the (cos(θ), sin(θ), 0) direction.
We calculate horizontal spatial derivatives ofu then plug them into Equation A.5 to calculate the pointwise axial strain rate in the (cos(θ), sin(θ), 0) direction,σ θ . First, let's calculate the point-wise axial strain of a horizontal fiber oriented in the (C θ , S θ , 0) direction at the point (x, y.z). The first steps are to calculate our spatial derivatives of the particle velocity, keeping in mind thatux = ickAS φ 1 oP S anduy = −ickAC φ 1 oP S , where oP S = e ik(ct−xC φ 1 C φ 2 −yS φ 1 C φ 2 −zS φ 2 ) : Now we can plug these in for the point-wise axial strain in the θ direction: (iv) We averageσ θ over a horizontal line segment of length g, the gauge length, yieldingσ θ,g , which mimics a DAS signal in the same location in the direction (cos(θ), sin(θ), 0). σ θ,g (x, y, z, t) = 1 g (v) We verify that lim k→0σ θ,ġ σ θ = 1, i.e. at long wavelengths, the average strain rate over a fixed gauge length acts like a point-wise axial strain rate measurement in the same direction.

APPENDIX C: DERIVATIONS OF RESPONSE OF DAS AND GEOPHONES TO POINT SOURCES
Here, we derive analytical expressions for idealized geophones (particle velocityu θ ) and DAS channels (point-wise axial strain rateσ θ ) responding to a surface point source at xs that emits a wave that spreads in a circle along the Earth's surface. The wave mimics either a Rayleigh wave with particle motion in the direction of propagation and vertically, or a Love wave with particle motion orthogonal to the direction of propagation. Like real surface waves, we assume a 1/R die-off of amplitudes as the waves propagate away from their source locations, but unlike real surface waves, we assume no dispersion for the sake of analyzing a simple analytical expression.

C1 Love waves
We assume a simple model of displacement u at any surface point x = (x, y, 0) responding to a surface Love wave point source at xs = (xs, ys, 0): a frequency f Ricker wavelet with particle motion in the direction tangent to the circular spreading surface that decays as 1/ x − xs 2. The equation describing this displacement is: where c is the speed of propagation and T = (x − xs) × (0, 0, 1), soT = T/ T is the horizontal unit vector orthogonal to the direction of propagation, so T = (y −ys, xs −x, 0)/ x−xs 2. Taking the time derivative of the displacement results in the particle velocitẏ u: so the measurement detected by a geophone's x-component would bė and the measurement detected by a geophone's y-component would bė To calculate the point-wise strain axial strain rate detected by a fiber oriented in the i direction, we can use the chain rule with x − xs as an intermediate variable, ∂ux(x,t;xs) ∂x = ∂ux(x,t;xs) where the second term can be calculated as ∂ x−xs ∂x = ∂x (x − xs) 2 + (y − ys) 2 = x−xs √ (x−xs) 2 +(y−ys) 2 = x−xs x−xs and the first term can be calculated using the product rule (defining the notation R := x − xs ) ∂ux(x, t; xs) so we can put the two terms together to get the i direction axial strain rate: ∂ux(x, t; xs) ∂x = exp −π 2 f 2 (t − R c ) 2 (y − ys)(x − xs) R 3 To calculate the point-wise strain axial strain rate detected by a fiber oriented in the y-direction, we again use the chain rule with x−xs as an intermediate variable where the second term can be calculated as ∂ x−xs ∂y = ∂y (x − xs) 2 + (y − ys) 2 = y−ys √ (x−xs) 2 +(y−ys) 2 = y−ys x−xs and the first term can be calculated using the product rule (defining the notation R := x − xs ) ∂uy(x, t; xs) Then when we combine the two chain-rule terms to get the axial strain measurement in the y-direction, we get ∂uy(x, t; xs) ∂y = (xs − x)(y − ys) C2 Rayleigh waves Assume the horizontal particle displacement, u, at x to a point source at xs will be a Ricker wavelet that decays as 1/distance: where c is the speed of propagation, f is the peak frequency of the Ricker wavelet, and R = x − xs, soR = R/ R is the horizontal unit vector orthogonal to the direction of propagation. While a Rayleigh wave also has a vertical displacement component, we ignore that term because it would not mix with the horizontal components as observed either through a horizontal strain rate or velocity measurement. Taking the time derivative of the displacement, we get the particle velocityu: The measurement detected by a geophone's x-component would be: ux(x, t; xs) = x − xs x − xs 2 exp −π 2 f 2 (t − x − xs c ) 2 −6π 2 f 2 (t − x − xs c ) + 4π 4 f 4 (t − x − xs c ) 3 (C.9) and its y-component is: uy(x, t; xs) = y − ys x − xs 2 exp −π 2 f 2 (t − x − xs c ) 2 −6π 2 f 2 (t − x − xs c ) + 4π 4 f 4 (t − x − xs c ) 3 (C.10) To calculate the point-wise strain axial strain rate detected by a fiber oriented in the i direction, we can use the chain rule with x − xs as an intermediate variable ∂ux(x,t;xs) ∂x = ∂ux(x,t;xs) ∂ x−xs ∂ x−xs ∂x so the second term can be calculated as x − xs x − xs (C.11) and the first term can be calculated using the product rule (defining the notation R := x − xs ) ∂ux(x, t; xs) When we combine the two derivative terms, we get point-wise axial strain rate in the i direction as: To calculate the point-wise strain axial strain rate detected by a fiber oriented in the y-direction, we can use the chain rule with x − xs as an intermediate variable and the first term can be calculated using the product rule (defining the notation R := x − xs ) ∂uy(x, t; xs) and when we combine these two terms we get the yy-component of the strain rate: ∂uy(x, t; xs) ∂y = exp −π 2 f 2 (t − R c ) 2 (y − ys) 2 R 3 ) 4 (C.14)