Compaction of the Groningen gas reservoir investigated with train noise

S U M M A R Y Induced seismicity in the Groningen gas field in the Netherlands has been related to reservoir compaction caused by gas pressure depletion. In situ measurement of compaction is therefore relevant for seismic hazard assessment. In this study, we investigated the potential of passively recorded deep borehole noise data to detect temporal variations in the Groningen reservoir. Train signals recorded by an array of 10 geophones at reservoir depth were selected from the continuous noise data for two 5-month deployments in 2015. Interferometry by deconvolution was applied to the high-frequency train signals that acted as stable, repetitive noise sources. Direct intergeophone P and S wave traveltimes were then used to construct the Pand S-wave velocity structure along the geophone array. The resulting models agree with independently obtained velocity profiles and have very small errors. Most intergeophone P wave traveltimes showed decreasing traveltimes per deployment period, suggestive of compaction. However, the retrieved traveltime changes are very small, up to tens of microseconds per deployment period, with uncertainties that are of similar size, about 10 microseconds. An unambiguous interpretation in terms of compaction is therefore not warranted, although the 10 μs error per 5-month period is probably smaller than can be achieved from active time-lapse seismic surveys that are commonly used to measure reservoir compaction. The direct P-wave amplitudes of the train-signal deconvolutions were investigated for additional imprints of compaction. Whereas the P-wave amplitudes consistently increased during the second deployment, suggestive of compaction, no such trend was observed for the first deployment, rendering the interpretation of compaction inconclusive. Our results therefore present hints, but no obvious effects of compaction in the Groningen reservoir. Yet, this study demonstrates that the approach of deconvolution interferometry applied to deep borehole data allows monitoring of small temporal changes in the subsurface for stable repetitive noise sources such as trains.


I N T RO D U C T I O N
The Groningen gas field in the Netherlands is one of the largest onshore gas fields in the world. Production started in 1963 and the first recorded induced earthquake occurred in 1991 (van Eck et al. 2006). Seismicity remained rather low in the following years, with about five earthquakes (1.5 ≤ M L ≤ 2.7) per year, but increased with increasing production between 2003 and 2013, the year of the largest (M L =3.6) earthquake to date (Van Thienen-Visser & Breunese 2015). It has been suggested that reservoir compaction due to pore pressure reduction by gas extraction plays a dominant role inducing the seismicity (e.g. Bourne et al. 2014;Candela et al. 2018). The region of highest seismicity indeed roughly corresponds to the region of largest subsidence (Fig. 1). Monitoring the temporal variations of the reservoir is essential to calibrate mechanical models that relate compaction to the extraction of gas. Experimental studies have shown that the mechanical behaviour of the reservoir rock in response to pressure depletion is more complex than simple elastic shortening (Pijnenburg et al. 2018(Pijnenburg et al. , 2019. It is therefore important to have in situ measurements of reservoir compaction. Some measurements have been obtained from downhole radioactive markers (Kole 2015) and, more recently, from distributed strain sensing (DSS, Cannon & Kole 2018). It would be beneficial to investigate if compaction can be determined independently from seismic data.
A large number of studies have shown that it is possible to detect temporal changes in subsurface properties from changes in the seismic wavefield. Seismic velocity changes associated with earthquakes have been measured in various ways: actively in a cross-well experiment (Silver et al. 2007;Niu et al. 2008) as well as passively using different approaches (e.g. Yamada et al. 2010;Nakata & Snieder 2011;Pei et al. 2019). Time-lapse 3-D seismic surveys revealed traveltime changes in and around oil or gas reservoirs (e.g. Hatchell & Bourne 2005;MacBeth et al. 2018) and reservoirs for CO 2 storage (e.g. Boait et al. 2012). Minor changes in the medium can be inferred from the coda of ambient noise autoor cross-correlations, and studies of coda-wave interferometry have revealed seismic velocity changes caused by the solid-earth tides (Sens-Schönfelder & Eulenfeld 2019) and earthquakes (e.g. Wegler & Sens-Schönfelder 2007;Brenguier et al. 2008), in buildings (e.g. Nakata et al. 2013) or other structures (e.g. Salvermoser et al. 2015). However, coda-wave interferometry has limitations localizing the sources of the detected changes in the medium, although approaches have been suggested to solve this problem (Larose et al. 2010;Obermann et al. 2019). Compared to coda-wave studies, interferometric studies that use direct (ballistic) waves to identify temporal variations provide a better spatial resolution. Recently, Brenguier et al. (2020) and Mordret et al. (2020) inferred temporal velocity variations in the shallow Groningen subsurface from ballistic body and surface waves obtained from ambient noise data of a dense surface geophone network. Variations in the 3-km-deep subsurface can not be obtained from surface sensors, but Behm (2017) suggested that it should be possible using noise recordings from deep downhole arrays.
In 2013, two monitoring wells in the Groningen gas field were equipped with geophone arrays in the reservoir at 3 km depth. Zhou & Paulssen (2017) determined the P-and S-wave velocity structure along the geophone array of one of these boreholes with noise interferometry by cross-correlation. They further found that diurnal variations in anthropogenic noise had a strong impact on the interferograms, impeding the detection of temporal velocity variations in the reservoir. However, they suggested that high-frequency repetitive noise signals generated by passing trains might allow more stable and higher accuracy traveltime measurements. Trains are well-known sources of strong vibrations (Chen et al. 2004;Fuchs & Bokelmann 2018), and seismic interferometry has previously been applied to train signals (Quiros et al. 2016;Zhang et al. 2019). Moreover, Brenguier et al. (2019) showed the potential of train-signal interferometry for monitoring temporal velocity changes associated with active faults. To our knowledge, the study presented here is the first that employs train noise recorded by a deep borehole geophone array.

B O R E H O L E DATA A N D T R A I N S I G N A L S
The data that are used in this study are from a deep borehole geophone array at reservoir level at 3 km depth. The operator of the gas field, the Nederlandse Aardolie Maatschappij (NAM), deployed a string of 10 geophones in well SDM-1 during two 5-month periods in 2015 (January 23-June 29 and July 3-December 1). The well is located in the Loppersum area of the Groningen gas field, a region that experienced a high level of seismicity and subsidence (Fig. 1). Three-component 15-Hz geophones were positioned at depths from 2750 to 3000 m (Fig. 2a), with 30 m intergeophone spacings from the top geophone (GP01) to geophone 9 (GP09) and 15 m between GP09 and lowest geophone (GP10). The geophones are positioned at the same depths for the two deployments and there is only a small time gap of 4 d between the two deployment periods. The setting is presented in Fig. 2(b) showing the location of the borehole at ∼500 m distance from the railway track which links the towns Stedum and Loppersum. Note that the vertical distance between railway track and geophone array is six times larger than the horizontal distance, implying close-to-vertical propagation of train noise along the geophone array.
Two examples of train noise recordings are presented in Fig. 3   ( Fig. 3c) and the next train traveling in the opposite direction (Fig. 3d). The spectrograms of the bottom geophone (GP10) are shown in Figs 3(a) and (b). Despite the large depth of the geophone array (3 km), train signals can be clearly distinguished in the spectrograms (Figs 3a and b). In the time domain, the signals can not be recognized in the unfiltered data, but they can be identified after bandpass filtering between 30 and 120 Hz (Figs 3c and d). The spectrogram of Fig. 3(a) shows that the characteristic frequencies increase from ∼20 to ∼80 s (after 12:09:09 on July 3, 2015) with maximum amplitudes between 60 and 80 s when the train is close to the borehole. The timing fits the train schedule which says that trains in the direction to Loppersum should depart from Stedum station at 12:08. The increasing frequencies are likely produced by train acceleration. Note that preceding the dominant wavetrain, there is a smaller amplitude signal at ∼30 s in Figs 3(a) and (c). The spectrogram of the next train is shown in Fig. 3(b) (starting at 12:15:55). In this case the dominant signal begins with high frequencies which gradually decrease, probably due to the deceleration of the train approaching Stedum. The timing agrees with the train schedule: trains from Loppersum should arrive at Stedum at 12:17. The secondary arrival is now observed after the dominant train signal, at ∼110 s in Fig. 3(b). The timing of the secondary wave trains with respect to the dominant signal for trains in opposite directions suggests excitation approximately halfway between Stedum station and borehole SDM-1. A satellite image shows that there is a switch at this location from the single track to the double track at Stedum station, which likely excites the secondary arrivals. Two additional figures are presented in the supplementary material to demonstrate the similarity of the recorded train signals.

T R A I N S I G N A L D E T E C T I O N A N D D E C O N V O L U T I O N
Unlike previous studies which used trains signals buried in continuous noise (Quiros et al. 2016;Brenguier et al. 2019;Zhang et al. 2019), we use isolated train signals to exclude other types of noise and signals. This requires identification and extraction of the train signals from the continuous noise data.
The train signals have the character of a wave train without a clear onset, their amplitude is slightly higher than the noise level in the 30-120 Hz frequency band and their duration is up to 100 s (Fig. 3). A first quick-and-dirty detection is obtained from the 30 to 120 Hz continuous, vertical component data of geophone 2: the average of the (absolute) amplitude of a 30 s moving time window should be higher than a threshold which is just above noise level. Then, the spectrogram around each potential event is obtained. For each spectrum (as a function of time) the total power within the 30-90 Hz frequency band is calculated. This time-dependent signal is then smoothed to allow clear identification of the main and the secondary train signals (their time difference is normally 30-40 s, see Fig. 3). When the secondary signal arrives before the dominant signal, the event is identified as a train from Stedum to Loppersum, if the order is reversed it is identified as a train in the opposite direction. In this way, out of the roughly 9000 detected trains for each 5-month period, the travel direction of approximately 7000 trains could be determined. Finally, a 20 s time window is determined around the maximum of the main signal. The choice of 20 s is based on the observation that the dominant signal often has this duration (from 60 to 80 s in Fig. 3).
The P-and S-wave response can be obtained by interferometry using cross-correlation or deconvolution (Snieder et al. 2006(Snieder et al. , 2009. Rather than applying cross-correlation as in Zhou & Paulssen (2017), deconvolution was chosen in this study because of the ability to preserve the high frequencies while eliminating the transient train time signal. Akbar et al. (2018) showed that the deconvolution method performs well on the geophone data of SDM-1.
In the frequency domain, the deconvolution of the jth component of geophone R (R j ) by the ith component of geophone S (S i ) is given by where D ji RS may be interpreted as an estimate of the Green's function response of a virtual source at the location of S acting in the idirection recorded by a receiver at the location of R in the j-direction. In practice, to preserve stability of the deconvolution, a water level is applied and the deconvolution is approximated by: where S i * (ω) is the complex conjugate of S i (ω), and i (ω) is its autocorrelation with a water level: The water level is taken as 0.01 per cent of the maximum spectral power of the autocorrelation (c = 0.0001). If the autocorrelation for a certain frequency has a value below this water level it is replaced by it. As train signals are strongest in the 30-90 Hz frequency band (Figs 3a and b). The deconvolutions are calculated for 30-90 Hz bandpass filtered data to avoid influence from the higher amplitude noise at frequencies below 30 Hz.

P -A N D S -WAV E R E S P O N S E O F T H E R E S E RV O I R
We retrieved P-and S-wave responses from the intergeophone deconvolutions. The P-wave responses are obtained from the vertical component deconvolutions, as the P wave predominantly propagates in the vertical direction along the deep geophone array. In the following, we illustrate the method for the data of the second deployment (July 3-December 1, 2015). Fig. 4(b) presents the vertical component deconvolutions of a single train using the signal of the top geophone [D Z Z RS (t), with S = GP01 and R = GP01-GP10]. It shows that a stable deconvolution response is already obtained for a single train signal. The direct downgoing P wave from GP01 is retrieved robustly, despite the slight oscillations that are present after the main peak. By averaging over 60 trains (approximately the number of identified trains per day), we obtained stable deconvolutions that are nearly identical to the ∼9000-train average for the entire deployment period (Fig. 4c). By stacking, incoherent oscillations are significantly reduced and reflections from the bottom and the top of the reservoir can be distinguished. Although the P wave coda is determined with sufficient accuracy to allow a more detailed investigation, we focus here on the direct P wave because it has a higher amplitude and is more stable.
Other train noise interferometric studies retrieved P or Rayleigh wave responses (Quiros et al. 2016;Brenguier et al. 2019;Zhou & Paulssen 2019). Here we show that the S-wave response can also be obtained from train signals. To obtain the S-wave response with the highest signal-to-noise ratio, we searched for the direction that gives strongest S-wave amplitude. To find this direction, we performed horizontal component deconvolutions for all azimuths for the geophone pair GP02-GP10, covering the largest distance in the reservoir. The horizontal component pairs were rotated from west (-90 • ) to east (90 • ) to find the largest amplitude. Fig. 5 shows that the largest amplitude is obtained at approximately 0.11 s for an azimuth of -35 • or, alternatively, 145 • . The timing matches the expected S-wave traveltime for the velocity model of Zhou & Paulssen (2017). Considering the uncertainty in the observed azimuth, the uncertainty in the orientations of the horizontal components (±15 • ), as well as the fact that the trains are moving sources, it is likely that the azimuth of -35 • roughly agrees with the direction from SDM-1 to the dominant source direction. If the S waves were excited by the horizontal movement of the train on the railway track, the maximum S-wave polarization is expected to be roughly parallel to the railway (the frictional forces would excite horizontal shear displacements parallel to the track). Instead, the obtained direction is roughly perpendicular, corresponding to the source-receiver direction, that is within the vertical plane of dominant wave propagation. With the polarization in this plane it is likely that the recorded S waves are SV waves generated by conversion from P waves in the overburden  of the reservoir. Because the relative orientations of the horizontal components with respect to each other are fairly well constrained (Zhou 2020), the horizontal components of all geophones were rotated to 145 • and 55 • . We found that the train signal deconvolutions between the geophone pairs for the 145 • direction indeed produced clear S-wave responses. Fig. 6(b) shows that the downgoing direct S wave, with an average velocity of 2100 m s -1 , can be obtained from a single train deconvolution. Similar to the P-wave response, the S-wave stack of 60 trains gives nearly identical results to the one of ∼9000 trains (Fig. 6c), although the S-wave response is smaller in amplitude compared to the P-wave response.

P -A N D S -WAV E V E L O C I T Y I N T H E R E S E RV O I R
The previous section showed that intergeophone P-and S-wave responses can be retrieved from deconvolutions of individual train signals. Thus, the P and S wave traveltimes for all possible geophone pairs can be measured from the peaks of the downgoing direct waves. The velocity structure can then be calculated from the traveltime data. Similar to Zhou & Paulssen (2017), a kernel density estimation (Botev et al. 2010) is used to obtain the probability density function of the measured traveltimes for each geophone pair and its maximum likelihood value is taken as the intergeophone traveltime. The traveltimes from all geophone combinations are then used in a linear least squares inversion to obtain the P and S wave traveltimes between neighbouring geophones with uncertainty. Subsequently, the P-and S-wave velocity profiles are calculated. Fig. 7 presents the P-and S-velocity profiles obtained from the two data sets for 2015. We note that the uncertainty in the P-velocity profile is very small with a maximum value of 20 m s -1 (0.4 per cent). The uncertainty in the S-velocity profile is somewhat larger, varying between 30 and 80 m s -1 (1-3 per cent). The small uncertainties illustrate the high accuracy of the traveltime measurements and their internal consistency. Furthermore, we investigated the difference in the P velocity profiles calculated for trains traveling in the two opposite directions and found this difference to be negligible, that is, overlapping within the uncertainty (see also Zhou & Paulssen 2019).  The P-velocity structure obtained by train noise interferometry matches the sonic log data provided by NAM very well demonstrating the suitability of the approach and the accuracy of the results. However, we observe that our average P-wave velocity for the Ten Boer claystone (the interval between GP02 and GP04) is lower than the sonic log measured in 1963. This was also found in our previous study which used ambient noise cross-correlations (Zhou & Paulssen 2017), but the uncertainty in the current study is much smaller. The cause of the smaller P velocity in the claystone obtained from our data compared to the sonic log data is not clear. Rather than an underestimate of the true velocity, an overestimate of the (apparent) velocity may be expected from interferometry. Such an overestimate will occur if the dominant wavefield does not propagate along the array but is incident at an inclined angle. It is also unlikely that the velocity in the claystone has decreased over the past 50 yr (for instance by the opening of cracks associated to earthquakes) because we do not observe a traveltime increase for 2015 in the time lapse data, as will be presented in the next section. Potentially, there is an effect caused by the difference in scale length. The sonic log measurements typically use a 10 kHz signal corresponding to a wavelength of ∼0.4 m, whereas our measurements have frequencies from 30 to 90 Hz corresponding to wavelengths of 40-120 m. This implies that our data are sensitive to a wide area outside the borehole. The fault map provided by NAM shows a fault with an offset of 40 m at a distance of 50 m from the borehole. We speculate that heterogeneity in the vicinity of the borehole, potentially related to faulting (and their damage zones), may explain the discrepancy between our results and the sonic log data.
The S velocity in the reservoir is just above 2000 m s -1 , which is in good agreement with the average value of 2280 m s -1 given by NAM for the entire Groningen reservoir (Romijn 2017). The velocity profile also agrees with the previous study by Zhou & Paulssen (2017) using east-component ambient noise the cross-correlations. The largest difference is again obtained for the formation of the Ten Boer claystone (GP02-GP04), where the current study yields a lower velocity than obtained from ambient noise. There are two effects that can play a role. The ambient noise data were obtained for lower frequencies and might therefore be more strongly affected by the overlying high velocity anhydrite layer. Secondly, there can be an effect of anisotropy. The S-wave anisotropy over the entire (mostly sandstone) reservoir was found to be approximately 4 per cent (Zhou & Paulssen 2017) and it may be larger in the claystone. In the presence of anisotropy, the difference in inferred S velocity may be caused by the difference in horizontal direction that was used for the two studies.

T E M P O R A L C H A N G E S O F P WAV E T R AV E LT I M E S A N D A M P L I T U D E S
In the previous section, we used isolated train signals to calculate intergeophone P and S wave traveltimes, and from those we calculated the velocity structure. In this section, we show that the data can also be used to determine temporal variations.
To obtain time-lapse intergeophone traveltime measurements with a higher accuracy than single trains, we stacked the deconvolutions of 30 consecutive trains travelling in the same direction. Traveltimes of the direct P waves are then measured from the 30-train stacks. Since there are typically about 30 identified train signals per direction per day, this approximately corresponds to one measurement per day for each direction. For the two deployments of 2015, this allows analysis of temporal variations over the year. Fig. 8 gives an overview of all P wave traveltime measurements. Each row shows the traveltime diagrams obtained from the stacked Downloaded from https://academic.oup.com/gji/article/223/2/1327/5890340 by Utrecht University user on 07 November 2020 Figure 9. (a) P-wave traveltimes from GP01 to GP02 and (b) from GP02 to GP08. The measurements are obtained from stacked deconvolutions for trains from Stedum to Loppersum (red) and from Loppersum to Stedum (blue). The two deployments are indicated by their gray scale background. Linear fits are applied to the data and the accumulated traveltime change per deployment period (with standard deviation) is noted. vertical-component deconvolutions with, from top to bottom, the signal of geophone GP01-GP09 (i.e. the denominator in the frequency representation of the deconvolution). These geophones can be interpreted as the virtual sources. The columns represent the virtual receivers, the signal which is deconvolved. Within each diagram, the background indicates the deployment period: light grey for January 23-June 29 and darker grey for July 3-December 1. The traveltime measurements for trains from Stedum to Loppersum are in blue and those for the opposite direction are in red. There are no data for GP09 for the second deployment because the vertical component was out of order. Furthermore, the data of GP10 are affected by a change in the level of the gas-water contact related to drilling activities during the second deployment (Zhou 2020) and are therefore not shown.
Some diagrams show distinct steps in traveltime from the first deployment to the next, for instance for the geophone pairs that have GP04 as virtual source or receiver. These are likely caused by small changes in geophone position between the two deployments. To investigate the long-term traveltime changes that might be related to reservoir compaction, we applied linear fits to the traveltime data. We did this per deployment period for each geophone pair and obtained independent estimates for trains from Stedum to Loppersum and the opposite direction. The traveltime changes over the entire deployment period are then estimated from the linear fits. Fig. 9 shows examples for the geophone pair GP01 → GP02 across the anhydrite-claystone interface (Fig. 9a) and for geophone pair GP02 → GP08 along the reservoir (Fig. 9b). Compilations of all traveltime changes per 5-month deployment period are presented in Tables 1 and 2. To assess the reliability of the results, we not only calculated the traveltime changes for the deconvolutions with the virtual source above the virtual receiver (data above the diagonal in Tables 1 and 2), but also for the virtual source below the virtual receiver (data below the diagonal). It can be verified that most of the measurements are in agreement within the error. The tables show that most intergeophone traveltimes decreased a few or a few tens of microseconds per 5-month deployment period. The 1-standard deviation uncertainties are typically around 10 μs. Since the traveltime changes are generally smaller or just slightly bigger than the Downloaded from https://academic.oup.com/gji/article/223/2/1327/5890340 by Utrecht University user on 07 November 2020 Table 2. P wave traveltime changes for July 3-December 1 in microseconds. Values for trains from Stedum to Loppersum are shown in red at the top and those for Loppersum to Stedum in blue at the bottom. uncertainty, it is difficult to draw hard conclusions. However, the general pattern of decreasing traveltimes (negative values in Tables 1 and 2) is in agreement with an increase in velocity caused by compaction.
We note that our temporal traveltime changes of up to tens of microseconds per 5-month period are much smaller than the traveltime changes that are typically obtained for time-lapse 3-D seismic surveys. These vary between 0.1 and 1.5 ms yr -1 for various sandstone reservoirs, although reservoir slowdown has also been observed (MacBeth et al. 2019). We point out that the traveltimes changes of the two methods can not be directly compared. The active time-lapse seismic surveys include the effects of strain, that is reservoir thinning (or possibly extension), as well as of changes in medium velocity. The traveltimes measured by our interferometric approach, on the other hand, are obtained relative to fixed geophone positions within the borehole. They do not include changes in layer thickness and only reflect velocity variations in the medium.
Compaction may also affect the amplitudes of the deconvolutions. Deconvolution interferometry largely eliminates the source signal, so when the noise sources are repetitive and stable, amplitude variations of the deconvolutions can reflect changes in the medium. The (maximum) amplitude of the direct P wave was measured from the stacked train deconvolutions. An overview of all the data is presented in Fig. 10. Generally, there are only small amplitude variations during the first half of 2015, whereas all geophone combinations show a gradual amplitude increase during the second half of 2015. During this period, the amplitude typically increases from ∼0.04 to ∼0.045, corresponding to a 12 per cent increase. It might be tempting to explain the amplitude increase by compaction, but the absence of such a trend during the first half of 2015 suggests that this interpretation is not warranted.

C O N C L U S I O N S
This study shows that it is feasible to very accurately determine direct P and S wave traveltimes along a geophone array by deconvolution interferometry of train signals. The study was carried out for a deep borehole geophone array in the Groningen gas reservoir to investigate effects of compaction. Intergeophone traveltimes measured from isolated train signals were first used to determine the velocity structure along the geophone array. The P velocity profile, retrieved from vertical component deconvolutions, was found to be in good agreement with sonic log data except for the Ten Boer claystone layer in the reservoir. There, our data yield a roughly 10 per cent smaller velocity than the sonic log, clearly exceeding the 0.4 per cent P velocity uncertainty estimated from the train data. The discrepancy can be caused by the difference in scale length of more than a factor of 100 between the sonic log data and our train signal data, and may reflect non-homogeneous structure around the borehole. The S velocity structure was inferred from horizontal component train signal deconvolutions. The S waves have maximum amplitude in the direction perpendicular to railway, suggestive of an SV wavefield generated by P-to-S conversions at lithological interfaces between the surface and the geophone array. The S velocity model from train signals is in good agreement with that obtained from a previous study using ambient noise cross-correlation interferometry (Zhou & Paulssen 2017) with the largest difference again for the claystone layer of the reservoir. Yet, we find that the S-velocity models are overall in good agreement considering the uncertainties in the previous study as well as potential effects of anisotropy (∼4 per cent, see Zhou & Paulssen 2017).
The train signal deconvolutions were then used to investigate temporal variations of the P-wave response. Most of the (daily stacked) intergeophone P wave traveltimes showed a slightly negative trend for the two 5-month deployments in 2015. Whereas this is in line with an interpretation of reservoir compaction increasing the P-wave velocity, many of the intergeophone traveltime changes were too small to be significant considering the uncertainties. It is noteworthy that the changes, which are up to several tens of microseconds per 5-month period, are smaller than obtained by time-lapse 3-D seismic surveys for a large variety of reservoirs. Our interferometric approach is only sensitive to variations in medium velocity (assuming that the geophones within the steel casing of the borehole do not move along with the surrounding rock), whereas the traveltime changes from the time-lapse seismic surveys include the additional effect of changes in reservoir thickness. The P-wave amplitude data of the deconvolutions also did not show clear evidence for compaction. The amplitudes consistently increased during the second deployment, as expected for compaction, but this was not observed for the first deployment. Thus, no unambiguous effects of compaction in the Groningen gas reservoir were found for the two 5-month deployments analysed in this study. Nevertheless, our results demonstrate that very small temporal changes in the medium can be detected by train noise interferometry applied to borehole geophone array data.

A C K N O W L E D G E M E N T S
We thank NAM for providing us with the data. Chris Spiers, André Niemeijer, Jianye Chen, Ivan Vasconcelos and Elmer Ruigrok are thanked for useful discussions. We thank Augusto Casas and an anonymous reviewer for their constructive comments that helped to improve the paper. Base map in Fig. 2(b) is retreived from OpenStreetMap. This project has been funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie actions grant agreement No 642029 -ITN CREEP.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Train signals recorded on 3 July 2015 at 12:41 (a, c) and 12:46 (b, d). These train signals are half an hour later compared to the ones shown in Fig. 3 of the main text. Figure S2. Train signals recorded on 27 November 2015 at 17:10 (a, c) and 17:16 (b, d). These train signal are more than 4 months later compared to the ones shown in Fig. 3 of the main text.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.