Quality parameters for passive image interferometry tested at the Groningen network

Passive image interferometry is widely used to retrieve velocity variations as function of time. The cross-coherence (spectrally normalized cross-correlation) of ambient noise recorded at two receivers is an estimate of the Green’s function that due to velocity changes will be stretched or compressed in time. The relative velocity change ( (cid:2) v / v ) is determined by the time stretching ( − (cid:2) t / t ) that yields highest correlation between the reference cross-coherence and stretched lapse cross-coherence. The estimation of (cid:2) v / v could be used, for example, to warn for hazardous situations developing in a volcano setting or due to degradation of a civil engineering structure. Before a warning would be issued, however, one would like to have a handle on the quality of the medium-change estimate. In the Groningen area in the Netherlands a well-sampled network of seismometers exists. This allows direct assessment of the quality of the (cid:2) v / v estimation by many different receiver paths sampling the same medium. We use this quality assessment of (cid:2) v / v to test other possible quality parameters that could also be used in settings with only a sparsely sampled seismic network. The quality of the (cid:2) v / v determination appears to be well described by the correlation coefﬁcient which is used to determine the velocity change. Consequently, in order to measure medium changes with a high certainty, it is important to have a high consistency between the lapse and the reference cross-coherence. A quality estimate based on this correlation coefﬁcient can also be applied when there is only one receiver pair available. If the quality is insufﬁcient, it can be improved at the cost of temporal resolution. This study also investigates possible causes of the measured medium changes: variations in temperature, soil moisture, air pressure, water table, gas production rate and subsidence. We ﬁnd a weak anticorrelation with temperature, and a weak correlation with the gas production rate and subsidence. The observed medium change is likely a complicated combination of different processes taking place.


I N T RO D U C T I O N
The possibility to detect seismic velocity changes is key to monitor natural phenomena with possibly hazardous implications. Before a volcanic eruption, when temperatures inside a volcano increase, and before an earthquake, when stresses increase, the seismic velocity changes. Monitoring these velocities may provide crucial information for predictions or better understanding of volcanic eruptions and earthquakes (e.g. Brenguier et al. 2008;Wegler et al. 2009;Obermann et al. 2014), just as the understanding of other tectonic deformations inducing stress changes (e.g. Brenguier et al. 2014;Takano et al. 2017). Also industrial activities in the subsurface lead to velocity perturbations. For example, the pressure front due to wastewater injection may be tracked by measuring velocity changes or the deformation development at an enhanced geothermal site (Hillers et al. 2015). Monitoring these velocity changes benefits both the economics and safety of these engineering projects.
Seismic velocities can be monitored actively or passively. Active sources such as explosives or seismic vibrators, however, are hard to use continuously. Furthermore, they are expensive and may be dangerous to use. Passive techniques provide viable alternatives. Passive image interferometry (PII; Sens-Schönfelder & Wegler 2006) is such a passive technique. It uses ambient noise to measure relative changes in the seismic velocity ( v/v). This procedure employs two major steps. First, the cross-coherence of ambient noise recorded at two receivers is calculated to estimate the Green's function (Wapenaar et al. 2010). The Green's function is the response at one receiver, caused by an impulse at the other receiver, generating waves for all possible paths to the first receiver. This Green's function estimate is constructed from one or a few days of data. We refer to it as the lapse Green's function. Second, the stretching or compressing of this lapse Green's function in time can be measured with respect to a reference Green's function and used to calculate v/v (Lobkis & Weaver 2003).
As outlined by Weaver et al. (2011) the estimated v/v contains a physical part, actual changes in the medium, but also there is a spurious contribution, for example, caused by changes in noise field. Both the estimate of v/v and its uncertainty vary with time. In order to use velocity monitoring as hazard or engineering tool, it is crucial to have a good estimate of the quality of the determination. This quality can directly be assessed when many station pairs sample the same medium with different azimuths. However, it is not obvious which quality parameter could be used when only one or a few station pairs are available.
To find a suitable generally applicable quality parameter, we apply PII to 1 yr of seismic noise recorded at the Groningen network . The Groningen network has a sufficient azimuthal distribution of station pairs to estimate the quality of the v/v determination directly, using the consistency of independent measurements sampling the same medium from different azimuths. This quality is compared with a suit of candidate quality parameters that could also be applied in networks with less optimal azimuthal coverage. The first candidate quality parameter is constructed from the correlation coefficients that are found after optimally stretching a lapse Green's function with the reference Green's function. The second is based on a quantification of spurious contributions and describes how well the lapse Green's function is estimated. The third is based on how isotropically the medium is illuminated from ambient seismic noise.
In the following we test the different quality parameters and we investigate possible causes of the medium perturbations as observed over Groningen. The next section outlines the setting of our field test. The section thereafter describes the implementation of PII and how the different quality factors are defined. This is followed by a description of the results of applying PII and the various quality parameters on actual data. Possible underlying causes of the observed medium changes are investigated in Section 5. We close the paper with a discussion and conclusions.

S E T T I N G
Groningen is a province in the northeast of the Netherlands. It hosts the largest gas field in Western Europe. Production started in the early sixties after an accidental discovery in 1959 (de Jager & Visser 2017). In 1991 the first earthquake was recorded which could be ascribed to pressure depletion of the reservoir. For better monitoring the induced seismicity, a regional seismic network was constructed in 1995 ). Since 2015 a new network is operational dedicated to monitoring the Groningen field: the G-network. This network contains of 70 stations placed in the near surface. At each station there is an accelerometer at the Earth's surface and there are geophones at 50, 100, 150 and 200 m depth (Hofman et al. 2017).
All data from the G-network is openly available via KNMI (1993). Throughout this work we use the continuous recordings at depth of 50 m. These instruments well record surface waves from various noise sources and are less affected by local anthropogenic activity than the surface accelerometers. The sensors are 4.5 Hz threecomponent geophones which noise recordings well capture the Peterson New Noise Model down till about 0.3 Hz (Spica et al. 2018).
Below this frequency the instrumental noise starts to dominate over seismic noise.
From the 70 stations, we select a subset of stations in the middle of the Groningen field (Fig. 1a). This subset contains of 14 stations that sample an area of roughly 14 × 14 km 2 . An area has been chosen that may be considered 1-D for the upper 1700 m (Duin et al. 2006); It has unconsolidated Cenozoic sediments for the first 800 m, overlaying a thick layer of chalk (900 m). Below there is a thin blanket of Jurassic and Triassic sediments in the west (100 m) overlying a thick layer of Permian evaporites (1 km). In the east the evaporites are only half the thickness, but there are about 600 m of Jurassic and Triassic sediments. Below the evaporites, at about 3 km depth, the gas reservoir resides in a Permian sandstone formation. 91 unique station pairs sample the area with a rich distribution in azimuths (Fig. 1b).
Below 1 Hz, the noise field is dominated by the microseism (Longuet-Higgins 1950). The double-frequency microseism peak falls out of the sensitive range of the instruments. However, Groningen is close to the shallow the North Sea, which yields ocean-wave solid-earth interaction that can be recorded with frequencies beyond 1 Hz. The dominant noise direction is northwest, which corresponds to the location of the North Sea and North Atlantic with respect to the array (Kimman et al. 2012;Spica et al. 2018). From 1 Hz onwards anthropogenic sources start to dominate, like cars, trains and industrial activity. Wapenaar et al. (2010) derived that the cross-coherence H of the displacements at two receivers provides a fine estimate of the Green's function if the illumination is good. The cross-coherence is a crosscorrelation of the waveforms with additionally a spectral normalization; the cross-correlation is normalized with the noise amplitude spectra at the two receivers. This normalization works well to whiten the spectrally unbalanced noise cross-correlation. In this cross-coherence H(x A , x B , t) positive times represent the causal Green's function (i.e. the response at x B due to a pulse at x A ) and negative times represent the anticausal Green's function (i.e. the time inverted response at x A due to a pulse at x B ).

Passive image interferometry
Typically a reference Green's function (H ref ) is estimated by stacking cross-coherences over 1 yr of data. Also for time-periods of a single day, cross-coherences are stacked to obtain an estimate of the state of the medium on that exact day: the lapse Green's function H lapse . Relative changes in propagation velocity would stretch or compress the lapse Green's function with respect to the reference Green's function. The time lapse seismic imaging technique (Lobkis & Weaver 2003;Weemstra et al. 2016) determines a correlation coefficient CC between H lapse stretched in time with factor (1 − ), and H ref : The stretching factor for which the correlation coefficient is highest serves as an estimate of the relative velocity change v/v = . t 1 and t 2 control the time window that is taken from the crosscoherence to estimate CC and . v/v is determined more accurately when this time window is chosen such that in the cross-coherence the part of direct waves is disregarded (Poupinet et al. 1984;Sens-Schönfelder & Wegler 2006;Weemstra et al. 2016), since these direct waves are most sensitive to the (daily varying) illumination of the medium. The time window including the direct wave is chosen as t 1 = −t max and t 2 = t max , where t max is the maximum lag time for which the cross-coherence is calculated. For the time window excluding the direct wave integrals where t direct is the traveltime of a wave travelling a direct path between the receivers, depending on the phase velocity and the distance between the receivers. Traveltime t direct is calculated with a marge to ensure that all energy of the direct waves is disregarded.
The average v/v over multiple receiver pairs is calculated as an azimuthally weighted average, together with an azimuthally weighted standard deviation. The weights are proportional to the range of azimuths that a receiver pair represents (i.e. the smallest angle between the two receiver pairs divided by two).
For the application of PII one could use the freely available Python package created by Lecocq et al. (2014) which is suitable to monitor velocity changes from ambient noise.

Quality assessment parameters
Under the assumption of velocity changes being isotropic, the different receiver pairs sampling the same medium from different azimuths should detect the same v/v. Detecting differences in the velocity changes would then indicate that the velocity changes were not determined accurately. Thus, the standard deviation of the velocity change over the different receiver pairs sampling the medium from different azimuths is used as a base to derive a consistency quality parameter Q PII (t): where the daily varying standard deviation σ lapse v/v (t) quantifies the differences between the measurements, and the standard deviation of the yearly variations σ year v/v describes the yearly perturbation of the medium and acts as a normalization parameter. These standard deviations are defined as is the velocity change averaged over all n receiver combinations and v/v is v/v(t) averaged over a whole year. The quality factor approaches 1 when the measurements of all receivers pairs give very similar year v/v . A quality factor much smaller than 1 will be reached when there is a large inconsistency within the network of the estimated daily v/v i (t). Therefore, consistency quality factor Q PII (t) is a good measure for the reliability of the daily estimations of velocity change v/v.
The second quality factor Q CCF (t) is constructed from the highest correlation coefficient calculated using eq. (1) and averaging over all receiver combinations. The coefficient quality is given by Here CC i (t) represents the largest correlation coefficient CC( ) (eq. 1) used to calculate v/v(t) between receiver combination i, and n represents the amount of receiver combinations. This quality factor approaches 1 when the optimally stretched lapse cross-coherence is a perfect match of the reference cross-coherence. The quality factor will be lower than one when the lapse cross-coherence can not well be written as a stretched copy of the reference cross-coherence. Strong correlation between consistency quality Q PII and coefficient quality Q CCF would indicate that a high similarity between the lapse and reference Green's function is needed to estimate v/v well.
The third quality factor uses the power distribution of noise sources around the Groningen setting. We use the method of Ruigrok et al. (2017) to calculate the power of waves coming from different directions. Their method of cross-correlation beamforming requires as input the cross-coherence that we also used in the previous subsection. The beamforming results, however, are aliased due to the relatively large distances between the receivers. We limit the effect of the aliasing by only focussing on the beampower in the slowness range in which surface waves exist. To determine sensible slowness bounds we first determine the phase velocities of these surface waves.
For this purpose we apply multichannel analysis of surface waves (MASW) as introduced by Park et al. (1998). MASW uses the moveout of the surface waves over an array of receivers to obtain average phase velocities as function of frequency. As a first step we transform the estimated reference Green's functions to the wavenumber frequency domain. The amplitudes squared represent the power as function of frequency and wavenumber. The wavenumber axis is converted to a phase velocity axis using v ph (ω) = ω/k(ω), where ω is the angular frequency and k is the wavenumber. The connections between maxima of the power as function of phase velocity and frequency represent the dispersion curves of the surface waves. We smooth the dispersion curves and plot them on top of the power plot.
We use the slowness bounds of the dispersion curve with the highest power to improve the beamforming results by stacking the beamforming results of tight frequency and slowness ranges. With this stacked power P(θ ) as function of direction and its maximum P max , we define the illumination quality Q ILL (t) as the ratio between the average power of all directions and the maximum power: This quality factor will be 1 when all directions have the same power and approaches zero when all power is concentrated at one direction. Strong correlation between consistency quality Q PII and illumination quality Q ILL would indicate that an azimuthally balanced illumination is required to determine v/v accurately.
In a fourth quality assessment we compare, within the crosscoherences H(t), the power before the first arrival with the power after the first arrival. First, we use the highest velocity calculated with MASW, together with the distance between two receivers, to calculate the arrival time of the first surface wave, and we use the lowest velocity to calculate the arrival time of the last surface wave arriving through a direct path between the receivers. To allow the complete waveform to come in we take for τ 1 the first arrival time, minus a small time of 2 s, and for τ 2 the last arrival time plus an additional time of 10 s. Then we compare the average power before (causal and acausal) the first arrival P 0 = 1 2τ 1 τ 1 −τ 1 H (t) 2 dt with the average power of the surface waves P 1 = τ 1 H (t) 2 dt for every receiver combination for which τ 1 > 0 and calculate the average ratio P 0 /P 1 . Now we define the cross-coherence quality as The quality factor approaches 1 when the power of the response before the first arrival is much smaller than the power of the response of the actual waves, and will be negative when the power before the first arrival is larger than the power of the ballistic surface waves. Strong correlation between consistency quality Q PII and cross-coherence quality Q H would indicate that v/v can only be determined accurately if the Green's function is estimated well.

Green's function estimation and dispersion analysis
The estimation of the velocity changes and the quality there-off is done using lapse and reference Green's functions. Estimations of these Green's functions are shown in this section. 1 yr of continuous data is used in the analysis, starting on 2016 March 30. The reference Green's functions are further used to estimate geometric dispersion relations for surface waves in Groningen. These are needed for two of the quality parameters that are computed in the next section. Fig. 2 shows an example of lapse cross-coherences, together with the reference cross-coherence. For every day's calculation we stacked the cross-coherences calculated from 50 percent overlapping time windows of 20 min duration, where the first time window ranges from 0:00 to 0:20 UTC, the second from 0:10 to 0:30 UTC, etc. Fig. 3 shows the cross-coherences of all receiver combinations shown in Fig. 1, stacked over 1 yr. These reference cross-coherences are plotted as function of distance between the receivers for the radial (a), the transverse (b) and the vertical (c) component. The arrival times increase with distance, as expected. The transverse component shows arrivals with a larger moveout than on the radial component. This would indicate that Love waves are slower over the unconsolidated Neogene sediments than Rayleigh waves. Unlike the radial component, the vertical component shows clearly two arrivals, which implies two different modes of the Rayleigh wave.
Using the reference cross-coherences (Fig. 3), the phase velocities are calculated for the shallow subsurface of Groningen. Fig. 4 shows the power spectrum of the retrieved surface waves as function of frequency and phase velocity for Rayleigh waves (a) and Love waves (b). The power of the Rayleigh wave is constructed by summing the powers of the radial and vertical components. Several modes are found both for the Rayleigh and the Love wave, but different modes dominate for different frequency ranges. Only modes that exist consistently over the patch of receivers (Fig. 1) can well be picked. The velocities we find for the Rayleigh waves are comparable to the velocities found by Spica et al. (2018), but the Love-wave velocities we find are much lower.
From Fig. 4 we find that most energy propagates in the first higher mode with phase velocities between 850 and 1000 m/s, corresponding to frequencies of 0.65 and 0.45 Hz. Assuming that the maximum sensitivity of a wave is at one-third of the wavelength, the depths of maximum sensitivity range from 440 to 740 m. Therefore, most energy propagates within the 800 m deep package of Cenozoic sediments.

Estimating medium changes and the quality there-off
In this section we show the results of the application of PII and the different quality parameters as introduced in Section 3. Fig. 5 shows v/v obtained using the radial, transverse and vertical component, when in the stretching technique (eq. 1) the full cross-coherences are used (a, c, e), or the late arrivals only (b, d, f). Each curve represents an average of v/v of receiver pairs binned with orientation bins of 22.5 • . In each computation, both the anticausal and causal parts of the cross-coherences are used. When the direct waves are omitted (panels on the right-hand side of Fig. 5) larger velocity variations are estimated. In this case the estimated velocity variations are more consistent for receiver pairs with different azimuths. The calculation parameters used to obtain Fig. 5 are given in Table 1.
The averages of v/v over all station-pair azimuths are shown in Fig. 6. These are azimuthally averaged means that are calculated with the late arrivals only. Also the azimuthally weighted standard deviations σ lapse v/v (t) are indicated. We now have arrived at a point where we can assess whether the quality of the v/v determination is correlated with the other quality factors. From the data underlying Fig. 6 we compute consistency quality Q PII (t) in accordance with eq. (2) and coefficient quality Q CCF (t) in accordance with eq. (3). The comparisons between these quality factors are shown in Fig. 7 for the three different components, using the procedure excluding the arrivals of the direct surface waves. A first inspection reveals a strong correlation between Q PII (t) and Q CCF (t) for the radial and transverse components. For the vertical component the sizes of the qualities do not match. The correlation coefficients between the quality factors, however, show a strong correlation for all components. The correlation coefficients are shown in Table 2 for both the procedures excluding (Q PII, ex ) and including (Q PII, in ) the direct wave in the determination of v/v. The correlations are lower for the Qualities determined including the direct wave.
In order to compute illumination quality Q ILL (t) (eq. 4) we first need to compute the azimuthally varying beampower P(θ ) and maximum beampower P max . The beampower is calculated for frequencies for which the main phase velocity is determined well (Fig. 4). For every frequency bandwidth of 0.02 Hz, slowness bounds are calculated by inverting the phase velocity at the minimum and maximum frequency. Using these frequency bandwidths and corresponding slowness bounds, we calculate the power as function of time and angle with the north (clockwise). Fig. 8 shows this power of the waves recorded at the radial and vertical components, stacked over all frequency bandwidths between 0.34 and 0.73 Hz ([0.34 0.36] Hz, [0.35 0.37] Hz, etc.). At most frequencies the dominant source is in the direction of angles between 280 • and 350 • . This can be due to microseisms on the North Sea, which is located in the northwest, as seen from Groningen. A large change in illumination can be seen between days 13 and 98, which is a period with a large number of not working stations. Beampowers from this period, therefore, have not been used in further calculations.
From the data underlying Fig. 8 we derive the quality of the power distribution in accordance with eq. (4): Q ILL . We compare this quality to the consistency quality of the determination of v/v: Q PII . Fig. 9(a) shows the time-series of these quality parameters. The correlation coefficients between these qualities are shown in Table 2. These coefficients are below 0.1. Hence it can be concluded that a correlation is absent. Over the course of the deployment, different stations were inactive and therewith the array response function changed. This resulted in a varying amount of aliasing, effecting the Q ILL parameter, without actual change in illumination. A large effect can especially be seen between days 13 and 98 in Fig. 8. For this reason, we did not include the values of the quality Q ILL of these days in the calculations of the correlation between Q PII and Q ILL .
Due to an increase in spatial aliasing, there is some inconsistency between the frequency bandwidth used for computing the beampower ([0.34 0.73] Hz) and in the determination of the velocity variations ([0.30 1.0] Hz). We have assumed that outside the [0.34 0.73] Hz range the energy field was comparable to the field inside the range. We do however not expect a large effect on the correlation between the two quality factors if there is some deviation from this distribution outside the range mentioned.
From the lapse cross-coherences we compute the quality parameter that is indicative of the quality of the Green's function estimate. We compute this cross-coherence quality Q H in accordance with eq. (5) and compare this quality to consistency quality Q PII . The correlation coefficients between these qualities are shown in Table 2. These coefficients are below 0.2. Hence it can be concluded that a correlation is absent.

I N T E R P R E TAT I O N
In this section we search for underlying causes of the observed medium changes (Fig. 6). In the frequency band we used ([0.30 1.0] Hz) for detecting the medium changes, the surface waves sample primarily the top 800 m of unconsolidated Cenozoic sediments. These sediments might be affected by gas production from, and compaction of, the underlying reservoir. From above, these sediments undergo influence from the atmospheric state. Here we investigate a possible correlation both with gas production and meteorological parameters.
Seismic velocities could be changed by meteorological conditions through changes in water saturation, pressure and temperature (e.g. Wang et al. 2017). Fig. 10 shows the comparison of v/v with meteorological parameters from a nearby weather station (Eelde; Fig. 1a, [lon,lat,alt]=[6.585 • , 53.125 • , 5.20 m]). A weak anticorrelation (R = −0.32) can be seen between temperature and v/v. When temperature goes up, velocity goes down. Air pressure does not show any substantial correlation with v/v. The third meteorological parameter is soil moisture. Soil moisture is difficult to measure in situ. Instead a proxy is used: the measured precipitation minus a modelled evaporation. The soil moisture goes down during the growing season of plants, when much evaporation takes place  Table 1. The colours from this figure match the col ors of the lines between receivers in Fig. 1. The used reference cross-coherences are shown in Fig. 3.  Fig. 10c). This yearly cycle of soil moisture depletion and replenishment cannot be traced in v/v. We also test a possible correlation of v/v with the ground water table. Daily water levels in the area are obtained from http://www.di noloket.nl. Ground water variations vary from site to site. However, the yearly trend at all sites follows the same pattern as the soil moisture proxy. The water table fluctuations are less than 1 m over a year. As with soil moisture, no correlation can be seen between v/v and the ground water table.
Gas from the Groningen gas field is produced from 20 production clusters (http://www.namplatf orm.nl). Three of these production sites largely stopped production in 2014 (Muntendam-Bos et al. 2017). From the 17 remaining clusters, we select the ones that are well within the area that is sampled with PII (Fig. 1). We leave out production from clusters at the southern edge of this area. Fig. 11(a) shows both the total monthly production over these sites and v/v. Some correlation can be noted. Especially the production peak between days 200 and 300 coincides with a peak in v/v during the same period.
Over the year, overall there is an increase in v/v (Fig. 6). This could potentially be attributed to compaction. The production of gas from the reservoir leads to compaction of the reservoir rock, which is largely accommodated by subsidence of the Earth's surface. This subsidence is routinely monitored with spirit levelling, InSAR and GPS (van Eijs & van der Wal 2017). The current total subsidence is about 30 cm in the middle of the subsidence bowl, which is expected to grow to about 47 cm in 2080 (van Thienen-Visser & Fokker 2017). Potentially, the compaction at reservoir level also results in compaction of overlying Cenozoic sediments, especially near the middle of the subsidence bowl. The surface waves that we used in the application of PII sampled a significant portion of this 800 m package of unconsolidated sediments. Compaction would lead to a smaller distance of the grains and thus to an increase in seismic velocity. Fig. 11(b) shows the subsidence as derived from InSAR. The upgoing trend corresponds with the v/v signal, but the oscillations over the year show little correlation.

D I S C U S S I O N
In Section 4 we extracted velocity variations from noise crosscoherences over radial, transverse and vertical components. For all components a consistent v/v was found. The v/v functions, however, were most pronounced on the horizontal components (Fig. 6). This was also reflected in overall higher consistency qualities Q PII for the horizontal components (Fig. 7). Figure 6. Azimuth-weighted average of the velocity variation plotted per day with its standard deviation (blue squares with error bars) and smoothed with 5 d (red line), estimated with late arrivals only (a) for the radial, (b) for the transverse and (c) for the vertical components. Detailed calculation parameters are given in Table 1.
From the three quality factors we compared with Q PII in Section 4 we found that only coefficient quality Q CCF is descriptive of the quality of the velocity determination. This notion corresponds with findings of Hadziioannou et al. (2009) who found in a laboratory experiment that retrieving a high quality Green's function is not important for estimating velocity changes. Instead, what is important is that the coda of the Green's function is retrieved in a stable  Table 2. Table 2. Correlation coefficients (at zero lag) between quality factor Q PII (t) and quality factors Q CCF (t), Q ILL (t) and Q H (t). Quality factors Q PII, ex and Q PII, in (eq. 2) give the consistency qualities of the velocity measurement when the first arrivals are excluded and included, respectively. Quality factors Q CCF (eq. 3), Q ILL (eq. 4) and Q H (eq. 5) are the correlationcoefficient-based, illumination-based and the cross-coherence-based qualities, respectively. fashion. Also Hobiger et al. (2012) show, for a field study, a direct link between the quality of the v/v estimate and the correlation coefficient. Fig. 2 shows for the Groningen case that the lapse cross-coherence does not have a high quality. It is not symmetric in time and has many spurious arrivals at acausal times, before the arrival of the surface wave. Hence, it is not a good estimate of the actual Green's function. This bad quality lapse cross-coherence can be explained by the biased illumination from the northwest as found in Fig. 8. Nevertheless, the coda that is retrieved in Fig. 2 is very stable and is consistent with the reference cross-coherence. This yields the overall high coefficient quality values (Q CCF ) when excluding the direct wave (Fig. 7). The match of Q CCF with Q PII confirms that it is the similarity between the lapse cross-coherence and the reference cross-coherence that determines the quality of the estimated medium change. This further supports the PII uncertainty parameter derived by Weaver et al. (2011) that is based on the CC (eq. 1).
The direct wave is often excluded in PII studies (e.g. Hobiger et al. 2012) as it is more sensitive to changes in illumination. One might expect therefore that the PII signal becomes less oscillatory when the direct wave is left out. In Fig. 5 it can be seen that the variance over the receiver pairs decreases when the direct wave is left out. This is consistent with the reduced sensitivity to illumination. However, the amplitude of the PII signal increases when the direct wave is left out. This unexpected observation might be a specific feature of applying PII in a soft-sediment setting where most energy of the ballistic Rayleigh wave is in the first overtone (Fig. 4). The coda waves might be dominated by fundamental-mode reverberations which are more sensitive to shallower structure where changes are more likely to occur.
Similarly remarkable in Fig. 5 is the higher amplitude on the horizontal components. This might indicate that the changes are especially accrued in the surface-wave coda and that the contribution of P-wave coda on the vertical results in a damping of the v/v signal.
We did not find a strong correlation between meteorological variables and v/v. Only a weak anticorrelation of v/v with tem-  (2) (red; direct wave excluded), illumination quality Q ILL in accordance with eq. (4) (a; black) and cross-coherence quality Q H in accordance with eq. (5) (b; black) for the radial component. The correlation coefficients between the quality factors are shown in Table 2. perature could be noted (Fig. 10). The absence of a strong correlation with the meteorological state might be attributed to the low frequencies used (<1.0 Hz) which have a limited sensitivity to the top soil that is changed most by meteorological conditions. Moreover, changes in the water table are quite small over the year in Groningen (less than 1 m). This is much less than the approximately 40 m variation that was found on Merapi volcano (Sens-Schönfelder & Wegler 2006), where these ground water fluctuations were found to be the major contributor to the observed v/v. We did find some correlation between v/v and the monthly gas production within the area sampled with PII. The large peak in the winter season could potentially be explained by a rise in production. Also we found an upgoing trend in v/v which might be linked to cumulative gas production through pressure reduction at reservoir level and compaction of both the reservoir and overlying sediments. For both observations, application of PII over a few years of data would be needed to exclude a coincidental match.

C O N C L U S I O N S
In the course of this research we determined dispersion relations of seismic waves for the shallow Groningen subsurface within fre-quency bandwidth [0.30 0.95] Hz. Both for Rayleigh and Love waves we found multiple modes.
Variations of the seismic velocity ( v/v) in the shallow subsurface of Groningen were retrieved using PII. We employed two methods: one that included and the other that excluded the direct arrivals. Having station pairs sampling the same medium with many different azimuths allowed us to compute the quality (expressed as azimuthal consistency) of the determined v/v. The consistency quality of the results using the method excluding the direct arrivals were much better than those including them. We found consistent v/v for all components. However, the consistency of the estimates from the horizontal components were better than those from the vertical component.
We tested a quality factor based on the correlation coefficient between the lapse cross-coherence stretched in time and the reference cross-coherence. This quality factor did show high correlation with the azimuthal consistency of v/v. Therefore, also in a setting with a single receiver combination, a quality factor can be assigned to the estimated velocity change. The similarity of the lapse crosscoherence and the reference cross-coherence increases when the temporal resolution is decreased; one can improve estimation of medium changes at the cost of temporal resolution. Soil moisture proxy (mm) Figure 10. Comparison of estimated velocity change with meteorological variables detected in Eelde, NL (Fig. 1a). The black lines in the three panels are the mean medium change as estimated over the radial component (Fig. 6a). It is compared with the daily average temperature (a), the daily average air pressure (b) and a proxy for the daily average soil moisture (c). In the upper right-hand corner of the panels the normalized zero-lag cross-correlation coefficients (R) are listed.
Variations within the power distribution of the noise were determined using cross-correlation beamforming between all shown receiver pairs. This also resulted in a quality factor. We investigated the correlation between the illumination quality and the consistency quality. It appeared that these qualities were not correlated.
We also evaluated the power within the cross-coherences from before and after the expected first arrival. This also resulted in a quality factor that did not correlate with the consistency quality.
As in previous research we found that a balanced illumination is not needed to retrieve velocity variations of good quality using PII. Nor is it necessary to have a good quality lapse Green's function estimate. It is important, however, to have a large resemblance between the lapse Green's function and the reference Green's function.
The observed medium change is likely a complicated combination of different processes taking place. We did not find one parameter that can largely explain the observations. We found a weak anticorrelation with temperature and we found a weak correlation with local gas production data and subsidence data. Subsidence since 01-01-2015 (mm) Figure 11. Comparison of estimated velocity change with gas production data (a) and the subsidence in Loppersum (b), which is at the northwestern edge of the used subarray (Fig. 1b). The black line is the mean medium change as estimated over the radial component (Fig. 6a). The red line in (a) depicts the monthly production in billion Normal cubic metres (10 9 Nm 3 ). The production is a sum of six production clusters that are within the area that is spanned by the used network (Fig. 1); the names of these production clusters are 't Zandt, Leermens, Bierum, Ten Post, Overschild and Amsweer. The red line in (b) shows the subsidence obtained with InSAR (bodemdalingskaart.nl).