Noise-based ballistic wave passive seismic monitoring. Part 1: body waves

F. Brenguier ,1 R. Courbis,1,4 A. Mordret ,2 X. Campman,3 P. Boué,1 M. Chmiel,1,4 T. Takano,1,5 T. Lecocq,6 W. Van der Veen,7 S. Postif3 and D. Hollis4 1Univ. Grenoble Alpes, Institut des Sciences de la Terre, Grenoble, France. E-mail: florent.brenguier@univ-grenoble-alpes.fr 2Massachusetts Institute of Technology, Department of Earth, Atmospheric, and Planetary Sciences, Cambridge, US 3Shell, Amsterdam, The Netherlands 4Sisprobe, Grenoble, France 5Tohoku University, Solid Earth Physics Laboratory, Sendai, Japan 6Royal Observatory of Belgium, Brussels, Belgium 7Nederlandse Aardolie Maatschappij, Assen, The Netherlands


INTRODUCTION
Large earthquakes and volcanic eruptions result from long-lasting, steady, pressure buildup on faults and magmatic reservoirs. However, the triggering mechanisms of impending events are thought to be initiated by short timescale stress and pore pressure transients associated with tectonic and volcanic interactions (e.g. Bouchon et al. 2011;Brenguier et al. 2014;Khoshmanesh & Shirzaei 2018) and possibly environmental perturbations (e.g. Johnson et al. 2017).
Anthropogenic activities such as hydrocarbon extraction, wastewater disposal, CO 2 storage and geothermal production also induce fluid pore-pressure related deformation that can lead to the triggering of induced seismicity (Talwani et al. 2007;Ellsworth 2013;Chang & Segall 2016). Monitoring these stress and porepressure perturbations continuously in time with high spatial accuracy at depth is thus critical to foresee forthcoming catastrophic tectonic and volcanic events and to improve reservoir management. C Seismic velocities are sensitive to stress and pore-pressure perturbations. This has been described in situ for a wide range of processes like for example solid earth tides (Takano et al. 2014) and rainfall induced pore-pressure changes (Wang et al. 2017). Niu et al. (2008) used active seismic source monitoring in boreholes at 1 km depth along the San Andreas Fault and were able to observe the first clear preseismic seismic velocity-stress perturbations in the hours preceding the occurrence of nearby earthquakes. Although convincing, these observations have never been reproduced due to the costly and logistically complicated in situ experiment.
Over the last decade, passive noise-based seismic monitoring proved success for monitoring volcanoes (Sens-Schönfelder & Wegler 2006;Brenguier et al. 2008a;Donaldson et al. 2017), earthquakes (Brenguier et al. 2008b) and environmental/climate (Lecocq et al. 2017) changes. Even though some attempts were made, no one succeeded in using passive seismic monitoring to observe clear preseismic anomalies similar to those described by Niu et al. (2008). The standard coda-based technique suffers from shortcomings that limit our ability to detect localized seismic velocity perturbations at depth. This technique also referred to as coda-wave interferometry (Poupinet et al. 1984) has the advantage of being very stable thanks to its low sensitivity to noise source property changes (Colombi et al. 2014). It thus allows detecting weak temporal changes of seismic velocities as small as 0.01 per cent such as those associated with solid Earth tides for example (Mao et al. 2019). However, the counterpart of this high detection capability is that the complexity of coda-wave propagation limits our ability to precisely characterize both the type (P or S) of velocity change and accurate estimates of their spatial distribution at depth.
In this work, we propose a complementary monitoring approach that uses ballistic waves reconstructed from noise correlations. This paper focuses on ballistic body waves and the companion paper Mordret et al. 2019 focuses on using ballistic surface waves on the same data set with the same type of approach. Body waves have a specific sensitivity to seismic velocity changes at depth and are potentially less affected by near-surface environmental changes than surface waves. By using ballistic waves instead of coda waves, we can more easily model the spatial sensitivity of the temporal change observations to local velocity perturbations at depth. The drawback of using direct, ballistic body waves instead of coda waves is their strong sensitivity to noise source temporal variations (Colombi et al. 2014). We use dense seismic networks and azimuthal averaging to circumvent this issue but still need to carefully analyse the stability of noise sources for such type of analysis.
Different studies showed that body-wave extraction from noise correlations is possible at various scales (Roux et al. 2005;Draganov et al. 2009;Poli et al. 2012). Nakata et al. (2015) were able to implement the first passive 3-D P-wave velocity tomography from continuous ground motion recorded on a dense array of more than 2500 seismic sensors installed at Long Beach (CA, USA). Recently Brenguier et al. (2016) and Nakata et al. (2016) proved the temporal stability of direct virtual body waves between dense arrays on Piton de la Fournaise volcano thus opening the way for continuous, passive ballistic wave monitoring. Brenguier et al. 2019 assess that the seismic noise emitted by heavy freight train traffic can be used to reconstruct repetitive virtual sources of body waves for monitoring the shallow crust at a few kilometres depth with a specific focus on the San Andreas Fault.
In this paper, we describe the fundamental aspects of passive ballistic wave monitoring using dense arrays and further illustrate its potential by applying it to a network of 417 seismic stations in the Netherlands. We are able to measure temporal changes of apparent velocities from both direct and refracted P waves, and thus we are able to separate the response of the near surface sediments and the basement located at 700 m depth. By providing direct observations of the rock mass response to stress changes at depth, this new passive seismic approach paves the way for improved in situ earthquake, volcano and producing reservoir monitoring.

METHODS
The new approach is based on measuring temporal changes of apparent slowness of specific ballistic waves that have been reconstructed from noise correlations using dense arrays of seismic sensors (Boué et al. 2013b;Mordret et al. 2014;Nakata et al. 2015Nakata et al. , 2016. The underlying requirement is that, as for Nakata et al. (2015), the high number of seismic sensors (> 100) and thus of noise correlation receiver pairs allows for the reconstruction of a virtual shot-gather of sufficiently high quality to be able to isolate and measure the apparent velocity of ballistic waves such as direct P or S waves, refracted waves or surface waves with clear mode separation (Mordret et al. 2019).
For this purpose of properly extracting high quality ballistic waves, we gather all possible noise correlations for all receiver pairs from a dense network into a single seismic panel (propagation time versus virtual source-receiver offset) thus considering a 1-D velocity model for which the apparent slowness of a reconstructed ballistic wave measured at surface can be written as: where n is the apparent slowness and V the apparent velocity. We can estimate n as the slope of apparent arrival times t along distance x under the assumption that the apparent velocity V is uniform along distance range x. We are now interested in measuring a temporal change in apparent slowness Ct n, the index Ct being for Calendar time ( Fig. 1a): The temporal change of apparent slowness can be measured directly as the difference of slowness estimates (from a slant-stack or beam-forming analysis for example) at different calendar times (de Cacqueray et al. 2016). Here we use an approach that estimates the temporal change of apparent slowness as the slope of the linear regression of the traveltime shifts at different calendar times for each distance offset along distance x ( Fig. 1b): where Ct t are the measurements of traveltime shifts at different offsets x for two different calendar times and for a specific windowed ballistic wave. This approach shows the advantage of providing a direct estimate of the uncertainty of the estimated apparent slowness temporal change by assessing how the measured traveltime delays Ct t fit a linear regression model along distance range x. The deviation of traveltime shifts, Ct t, to a linear trend along distance x can be explained by both noise affecting the reconstructed body waves and lateral variability of seismic velocity temporal changes that we do not account for in our 1-D approximated model. From eq. (2), we can derive the relation linking the temporal change of apparent slowness and apparent velocity: By multiplying the above equation by V, the apparent uniform velocity of the studied ballistic wave, this equation leads to (Fig 1c): This later equation shows that, for small velocity perturbations (< 10 per cent) we can estimate the relative temporal change of apparent velocity, Ct V V , of a given ballistic wave by estimating the value of the slope of the linear regression of Ct t measurements along traveltime t. This approach is sketched on Fig. 1. It applies to any kind of ballistic arrivals that can be clearly identified and isolated on a virtual source gather section. The case of direct surface waves requires mode separation and is challenging because of dispersion that leads to different velocities for different periods at which Ct t values are measured (Mordret et al. 2019). As suggested in Nakata et al. (2016) and Brenguier et al. (2019), this approach can also be applied to noise-correlations between two distant arrays in order to monitor diving body waves probing magmatic reservoirs, seismic faults or producing reservoirs at a few hundreds to a few kilometres depth. It is interesting to note that in this situation of two distant arrays referred to as A and B, the estimates of temporal changes of apparent velocities can be achieved using noise-correlations between arrays A and B and B and A separately, thus providing two independent estimates of temporal velocity changes. This method can also be applied using noise correlations between an individual seismic station and a distant array. In this work we apply this approach to the monitoring of a direct and a refracted wave using noise-correlations within a single array of about 8 km wide.

DATA
The Groningen gas field located in the northeast of the Netherlands is one of Europe's largest natural gas field. The reservoir located at 3 km depth is thought to be 40 by 50 km wide and 250 m thick. Bourne et al. (2018) show that the gas production in this field led to a 15 MPa average reservoir pore-pressure depletion since 1995 which is associated with seismicity rates that increased as an exponential-like function.
We use a network of 417 short period seismic stations deployed in the Groningen area of the Netherlands (Fig. 2) from 11 February (day 42) to 12 March (day 71) 2017 over a time span of 30 d. The network forms a grid array with aperture of the order of 8 km and an average station distance of about 400 m.
We computed an averaged seismic section of vertical to vertical noise cross-correlations. The noise-correlations were stacked in time (over the 30 d of continuous data), in space (along distance bins of 50 m long) and for the causal and acausal parts following Boué et al. (2013a) and Nakata et al. (2015). Fig. 3 illustrates these binned data for two different frequency ranges (1-3 Hz and 3-12 Hz). The low frequency (1-3 Hz) section mainly shows the propagation of low-velocity Rayleigh waves and also of ballistic P waves at velocities between 1.5 and 3 km s -1 . The lower panel of Fig. 3 highlights the high frequency (3-12 Hz) P waves interpreted as (1) the direct, diving P wave with velocity of ∼1700 m s -1 (see model on the right panel) contoured by a blue box and (2) a refracted wave at the 700 m interface with apparent velocity of ∼3300 m s -1 contoured by a red box. Thanks to the high stability in time of the useful high-frequency ambient seismic noise, we are finally able to reconstruct repetitive in time seismic sections from the correlations of daily records of ambient seismic noise leading us to daily virtual shot seismic gathers.

ANALYSIS AND RESULTS
In order to measure temporal changes, we further isolate the direct and refracted waves by applying a tapered window (Figs 4a and b). For the time shift measurements, Ct t, we use a cross-spectrum approach (Clarke et al. 2011) in the frequency range 3-8 Hz (best signal to noise ratio). We illustrate our approach on Fig. 4 by plotting traveltime shifts for the direct (left-hand panel) and refracted (righthand panel) waves measured using reference body waves stacked from days 42 to 50 and current body waves stacked from days 53 to 62 for the direct and days 48 to 57 for the refracted waves. Even though the traveltime shift measurements show large fluctuations likely associated with imperfect direct and refracted waves reconstruction and noise source changes through time, they show a clear linear trend along distance, especially for the refracted wave, indicating a clear change in apparent velocity between these two time periods. By multiplying the slopes of the linear regressions of the Ct t over distance measurements by the apparent velocities for both the direct (1700 m s -1 ) and refracted (3300 m s -1 ) waves (step b to c on Fig. 1), we find velocity changes of -0.25 and +1.5 per cent for the direct and refracted waves, respectively.
In order to gain insights on the temporal evolution of velocity changes, we further average the daily seismic sections using a 10-d long, 1-d moving window. This leads us to 21, full 10-d averaged, daily seismic sections (from days 42 to 71). Following the method described above, we measure velocity changes ( Ct V V ) between each 10-d averaged seismic section and the reference section corresponding to a stack of days 42-50 for both the direct (Fig. 5a) and refracted (Fig. 5b) waves. The apparent velocity change curves shown on Fig. 5 indicate a velocity decrease of maximum -0.25 per cent for the direct wave and a velocity increase of maximum 1.5 per cent for the refracted wave. The error bars correspond to the uncertainty of the linear regressions estimates of the traveltime shifts, Ct t, over traveltime t (Fig. 1c) using a least-square approach following Brenguier et al. (2008a).

THE NOISE STATIONARITY PROBLEM
The most common source of error in passive, noise-based seismic monitoring results from the non-stationarity of noise sources. Furthermore, ballistic waves are much more sensitive to noise source variations than coda waves (Colombi et al. 2014). As a result, the main drawback of this ballistic wave based monitoring method is the potential error introduced by changes in noise sources and great care has to be taken when interpreting the results. We thus need to properly assess how noise source temporal azimuthal variations can hamper our results. Due to the minimum interstation distance of 300 m we are not able to beamform the noise in the frequency range of interest (3-8 Hz) and can only beamform the 1-4 s period range. The average over 1-month beamforming result of Fig. 2   shows that the noise source distribution is strongly anisotropic with a main spot coming from the North Sea, north of our array (Fig. 2).
To study the HF (3-8 Hz) noise azimuthal distribution, we sort all noise correlations by azimuth (Fig. 6). On this figure, we observe that the HF body waves are best reconstructed for azimuth ranging from 120 • to 240 • meaning that the useful noise comes from the north which is in agreement with the 1-4 s beamforming observations. By looking in more details though, the 1-4 s beamforming points towards a dominant source northwest of the array and the azimuthally sorted correlations rather to a dominant source northeast of the array (best section is 180-240 • ). The noise in these two frequency ranges could thus be related to ocean activity and shore break but it is also possible that these noise sources have different origins with the 1-4 s noise being related to oceanic activity and the high frequency noise possibly coming from vehicle traffic from the city of Delfzijl located about 10 km north of the study area ).
Our assumption of comparing the high-frequency (3-8 Hz) to the low-frequency (1-4 s) noise sources might not be relevant but still informative. To study the stationarity of noise sources in the 1-4 s period range we beamform the ambient seismic noise on a 10-d average basis. We observe that during the time span of our analysis, the noise source distribution in the 1-4 s period range is mostly stable with small azimuthal variations of the centroid of the main spot by less than 10 • . We believe this is a good approximation for the upper limit of the noise source azimuthal variability for the high frequency range of 3-8 Hz. Indeed, vehicle traffic would probably have an even smaller range of spatial variations on a 10-d basis since the road and railway network is slightly modified over such timescales. Considering 10 • of noise azimuthal variations and the theoretical predictions of traveltime errors of ballistic arrival reconstructed from correlations of non-isotropically distributed noise sources from Weaver et al. (2009) and the further applications of Froment et al. (2010) and Colombi et al. (2014), we assess that, in case of a two sensors noise correlation, the error on the traveltime shift measurements (Fig. 5) would lead to a velocity change uncertainty of less than 0.5 per cent. In our case we emphasize that our traveltime shift measurements are obtained from azimuthally averaged noise correlations from the dense 417 stations array. We are thus confident that the observed velocity changes of + .5 per cent for the refracted wave is mostly related to physical velocity changes. However, the velocity decrease of 0.25 per cent for the direct wave might thus be partly biased by changes in the noise source properties.

INTERPRETATION AND C O N C L U S I O N S
In order to interpret the apparent velocity change of -0.25 per cent of the direct P wave, we can consider a simple 1-D model of wave propagation in the sedimentary cover and the ray theory that predicts that the observed apparent velocity change on the direct P wave corresponds to a real P-wave velocity perturbation averaged over the first 200 m depth (maximum distance between virtual source and stations of 2200 m). The only notable event during that time period is an episode of rainfall starting on day 51 in the region (40 mm of cumulated water height in 3 d). The decrease of velocity for the direct wave could thus be related to process of pore pressure diffusion as described by Rivet et al. (2015) and Wang et al. (2017). This interpretation is in good agreement with results from the companion paper Mordret et al. (2019) that uses ballistic surface waves obtained from noise correlations to monitor the shallow subsurface sediments. They find a clear S-wave velocity decrease following the rainfall episode starting day 51 that migrates at depth at a rate of 10 m per day down to 200 m.
Following again the ray theory, the increase of apparent velocity of 1.5 per cent for the refracted wave can be directly attributed to a change of P-wave velocity of the carbonate bedrock at 700 m depth. We checked for INSAR and GPS observations. These data do not show any clear transient anomalies during the time span of our analysis. We preclude the potential effects of swings in gas production in the main reservoir at 3 km depth below our study area due to the large distance to the probed carbonate layer. These induced pressure variations are of the order of less than 0.1 MPa locally and are thus too small to potentially induce a 1.5 per cent velocity increase 2.3 km above in the carbonate layer. We also rule out the effects of local induced earthquakes due to the low level of seismicity during the studied time period.
Considering this hypothesis of a bedrock velocity increase, one simple interpretation could be the effect of loading from rainfall on the bedrock that would increase the confining pressure and close cracks in the bedrock. However, it is unlikely that 4 cm of additional water height on days 51-53, corresponding to an increase of loading of 0.4 kPa, leads to a velocity increase of 1.6 per cent at 700 m depth. Indeed, following Yamamura et al. (2003) and Silver et al. (2007), and considering a high velocity-stress sensitivity of 10 −6 Pa −1 , the expected velocity change for a loading of 0.4 kPa should be about 0.04 per cent. There are three possible reasons for explaining this discrepancy. First, considering the error bars on Fig. 5, it is possible that the velocity increase for the P-refracted wave is smaller than 1.5 per cent. Secondly, considering that the carbonate bedrock is fully saturated, it is possible that its velocity-stress sensitivity is actually larger than 10 −6 Pa −1 . A large stress sensitivity value could be explained by high fluid pressure in this layer that reduces the effective confining pressure. And finally, it is possible that the average load on the study area is more than an equivalent of 4 cm of rainfall, possibly due to higher average precipitation on the area than the one locally measured. Canal drainage water flows could also be responsible for this increased load.
In conclusions, even though our results are hampered by high uncertainties and are spanning a too short period to be interpreted properly, we believe that this new passive monitoring approach has the potential for revealing seismic wave velocity temporal variations at localized areas at depth thus acting as a stress-strain probe. Even though the main drawback of this approach is its sensitivity to noise source changes and its requirement for dense seismic networks, we believe that future studies of the noise source bias effect together with recent step changes in seismic instrumentation (low cost sensors, Distributed Acoustic Sensing) will lead to groundbreaking advances in our understanding of natural and induced earthquakes, volcanic eruptions and will prove useful for reservoir management.

A C K N O W L E D G E M E N T S
We acknowledge the anonymous reviewers for their useful comments. This project received funding from the Shell Game Changer project HiProbe. We also acknowledge support from the French ANR grant T-ERC 2018, (FaultProbe), the European Research Council under grants no. 817803, FAULTSCAN and no. 742335, F-IMAGE and the European Union's Horizon 2020 research and innovation program under grant agreement No 776622 (PACIFIC). AM acknowledges support from the National Science Foundation grants PLR-1643761. The data were provided by NAM (Nederlandse Aardolie Maatschappij). We acknowledge M. Campillo, N. Shapiro, G. Olivier, P. Roux, S. Garambois, R. Brossier and C. Voisin for useful discussions. The authors thank Nederlandse Aardolie Maatschappij and Shell for permission to publish.