Groundwater Monitoring with Passive Seismic Interferometry

Abstract

Passive seismic interferometry takes advantage of natural ambient seismic noise generated by the wind, the storms and the human activities (e.g. cars, trains and hereafter pumps) to recover the slight variations of the seismic wave velocity induced by changes in the groundwater level. Here we compare the seismic measurements with actual piezometric data acquired on the Crépieux-Charmy (Lyon, France) groundwater exploitation field. We show the excellent correspondance between variations in the groundwater level and seismic velocity variations. We present hereafter the time and space monitoring of an hydraulic dome formed to prevent biological and chemical pollutions to enter the exploitation field. The horizontal resolution is solely limited by the number of seismic stations used, and is about 30 m in the present study. The vertical resolution of seismic measurement is impaired by spurious artifacts linked to the intermittent sources of noise. In average, the sensitivity of the seismic velocity change corresponds to a 50 cm change of waterlevel. This study confirms the possibility of groundwater monitoring in an industrial context with ambient seismic noise.

Share and Cite:

Voisin, C. , Guzmán, M. , Réfloch, A. , Taruselli, M. and Garambois, S. (2017) Groundwater Monitoring with Passive Seismic Interferometry. Journal of Water Resource and Protection, 9, 1414-1427. doi: 10.4236/jwarp.2017.912091.

1. Introduction

1.1. Area of Study

The Crépieux-Charmy exploitation field (Figure 1) is located in the alluvial plain of the Rhône River (North-East of Lyon, France). With up to 375 ha it is

Figure 1. The Crépieux-Charmy groundwater exploitation field is limited by the Miribel channel to the North, and the Vieux Rhône River to the South. The area of water uptake is bounded by infiltration basins used to produce hydraulic domes in case of biological or chemical pollutions in either of the limiting channels. The image focuses on basin 5.2 with positions of the piezometers installed by the exploitant (red dots). The seismic stations are deployed around basin 5.2. The red lines correspond to the direct paths between pairs of stations of interest. The blue dots labeled from A to H represent the locations of the virtual piezometers, at the center of each red path.

one of the largest groundwater exploitation fields in France. Started in 1957, the Crépieux-Charmy field now covers the water needs of the city of Lyon (2 million inhabitants).

The Crépieux-Charmy field is localized in the well-known Eastern aquifer of Lyon, a free aquifer characterized by glacio-fluvial Würmien’s deposits. In the particular zone of study, the glacio-fluvial deposits are covered by recent alluvial deposits related to the Rhône River. They describe a sedimentary sequence with mostly sandy-gravel size grain with a heterogeneous distribution classically found in such sedimentary environment. The recent alluvial deposits have a variable thickness between 10 and 25 m. The depth of the basement, represented by Pliocene glacial claystones and molasses is around 158 m.a.s.l [1] . The sedimentary environment and the texture of this formation also contribute to the complexity of the hydrogeological regime in the area of study.

In order to prevent biological or chemical pollutions to enter the groundwater uptake area, different series of infiltration basins have been created along the path of the Miribel and Old Rhône channels (Figure 1). In case of upstream incoming pollution, these basins are filled up with water and form a series of hydraulic domes below them. Each dome represents a local increase of a few meters of the water table, supposed to invert the local water flow from the channel to the uptake area. In the present study, we will focus on the hydraulic dome formed under the infiltration basin 5.2, used for different scientific and experimental studies by the research platform of “Eau du Grand Lyon”.

1.2. The Hydraulic Dome of Basin 5.2

The infiltration basin 5.2 has dimensions of 370 by 100 meters. It is bounded by an earthen dam of 4 meters high. The basin is equipped with several piezometers (Figure 1) among which P96 that is located in the vicinity of the supposed maximum of the hydraulic dome formed when the basin is filled-up with water. Piezometer B10, located on the western bank, is used as a control of the elevation of the hydraulic dome. French regulations enforce a minimum infiltration distance for water of 1 meter. Because of the uncertainty on the position of the maximum of the dome that represents the minimum infiltration distance, the water supply to basin 5.2 is stopped when B10 indicates a groundwater level of 2 meters below the surface level (mbsl). It is started again when B10 indicates a groundwater level of 2.6 mbsl. As a consequence, the water level inside basin 5.2 oscillates around a mean elevation of 1.8 meters. Figure 2 presents a schematic

Figure 2. Schematic functioning of the hydraulic barrier in Crépieux-Charmy field [1] .

representation of the functioning of the infiltration basin. Some questions arise about the shape of the hydraulic dome, the location of its maximum and the speed of infiltration that is known to be variable. The scarcity of piezometers installed in and around the basin 5.2 prevents to answer those questions. In the following, we explore the potential of passive seismic monitoring techniques to retrieve pertinent and additional information about the hydraulic dome.

1.3. Passive Seismic Monitoring of Groundwater

Passive seismic monitoring leans on the recording of the mechanical waves that propagate into the subsurface (Figure 3(a)). These waves originate from natural sources such as the wind, the storms, the rivers, and from anthropogenic sources such as cars, trains and hereafter pumps. The mechanical seismic waves are classified as body waves (P and S waves) and surface waves (Rayleigh and Love waves). All of them are characterized by a given velocity that depends on different parameters of the medium such as the shear modulus, the stiffness and the density for elastic parameters; fluid composition, porosity and saturation for poroelastic parameters. The waves velocity is accurately determined by cross-correlating the recordings of the ambient seismic wave field at two seismic stations. The result is formed by a series of waves arriving at different times, depending on their velocity and on the path they followed. We differentiate between the ballistic waves that propagate directly from station A to station B and the coda waves that propagate indirectly from station A to station B. The passive seismic monitoring consists on the repetition of this measurement over consecutive time windows. The results form a correlogram (Figure 3(b)). Any change of the subsurface properties linked to water table variations for instance will result in a slight change of velocity: the arrival time of the waves will be changed accordingly. The first application was shown on the Merapi volcano [2] . The authors demonstrated that the late phases (for time lapse greater than 2 s) reconstructed by the noise correlations are dominated by body waves. They also proved that the temporal decay of these waves is reminiscent of the decay of the coda waves, opening the way for noise coda-wave interferometry. Using precipitation data at the regional scale and a hydrologic model, they linked the annual velocity variations to variations in the water table. More recently, Voisin et al. developed an original approach based on the reconstructed ballistic surface waves to monitor the water table variations within a landslide [3] . The local noise sources are dominated by a river, a railway and a highway, all three being stable in time and space. The daily correlations of ambient seismic noise revealed seismic velocity changes up to ±1.5% that follow a summer/winter cycle consistent with the pore-water pressures monitored locally. This study suggests that seismic noise correlations are sensitive to water table oscillations through saturation changes and could be used as a nondestructive hydrologic monitoring tool when the distribution of noise sources is stable in space and time [4] .

In the following, we apply the very same technique in an industrial context,

Figure 3. (a) One hour of seismic noise recording at the 9 operating seismic stations. Each trace is normalized by its maximum to facilitate the reading of the plot. Note the arrival of large wave packets between 1000 and 1500 s and between 2500 and 3000 s. They correspond to the noise of pumping and injection of water in the basin 5.2. The largest arrival at 3100 s corresponds to a train passing to the north of basin 5.2. The rest of the traces is formed by numerous small amplitude wave propagating inside the medium and carrying information on its state; (b) Normalized correlation of the records of stations 5 and 6 over the duration of the experiment (vertical axis from 0 to 210 hours). The horizontal axis gives the time lapse of the different waves propagating between the two stations. We distinguish the ballistic waves arriving between −1 and +1 s and the coda waves arriving at later times. The positive side corresponds to waves travelling from station 5 to station 6; the negative side corresponds to waves travelling from station 6 to station 5. Note the strong assymetry of the correlogram dominated by the strong arrivals between 0.4 and 1 s. This indicates that most of the energy carried by the waves is travelling from station 5 towards station 6. Note the slight changes of the arrival time of the ballistic waves over the duration of the experiment. These changes are related to changes of the velocity of the seismic waves inside the medium induced by changes in the water level below the basin 5.2. The black rectangle shows very strong waves arriving near zero time lapse, i.e. almost simultaneously to stations 5 and 6. They occur between hours 120 and 140 and correspond to the infilling of basin 5.2. The origin of these waves is found in the water pumping and injection into basin 5.2.

where local noise sources are dominated by the pumping and injection of water, highly variable both in time and space. We address the following issues: is it possible to retrieve the water table variations in this configuration? If so, what would be the lateral resolution of the technique? And finally, what are the limitations of the technique imposed by the anthropogenic noise sources?

2. Comparison between Piezometers and Seismic Velocity Changes

2.1. Seismic and Piezometric Data

During 9 days from April 26, 2017 to May 06, 2017, nine seismic stations (Rauex by Sercel) equipped with vertical geophones with a nominative frequency of 5 Hzwere installed around basin 5.2. They continuously recorded the ground motion generated by the natural and anthropogenic noise sources at a sampling frequency of 250 Hz. The precise dating of the samples was insured by a continuous GPS synchronization of the seismic station. The waveforms recorded by each station are chopped into one hour windows (Figure 3(a)). During these 9 days of experiment, the basin 5.2 was emptied and filled alternatively to create a hydraulic dome under the basin 5.2. Evolution of the water height in basin 5.2 is monitored continuously (Figure 4). The height ranges from 169.4 m.a.s.l. (the bottom of the basin) to a bit more than 171 m.a.s.l. at its maximum. The water table level is monitored at different locations (Figure 1) every 10 minutes. Data are presented in Figure 4 for the piezometers B10, P96 and M25. Note the formation of the hydraulic dome that starts with the infilling stage of basin 5.2. The highest elevation is obtained for the piezometer P96, comforting its localization near the maximum of the dome. During the dry stage, the local groundwater flow is oriented from M25 (east) towards P96 and B10 (west). During the filled stage, the formation of the dome is associated with the reversal of the groundwater

Figure 4. Piezometric data (B10, P96, S20 and M25) and water level in basin 5.2. The red arrow indicates the period of seismic noise recording. Note the formation of the hydraulic dome when the basin is filled with water.

flow from P96 towards M25 and B10.

2.2. Methodology

2.2.1. Seismic Velocity Changes

For each seismic station and each one-hour trace, the signal processing operations are the following: the data is filtered between 3 and 20 Hz and is whitened in the spectral domain to correct for the differences of spectrum of the different noise sources. On a field experiment with standalone stations, it may happen that from time to time the seismic stations stop acquiring data; the missing data are replaced by zeros in order to complete the time series. Unfortunately in our experiment, station 1 stopped acquiring data early in the experiment. For this reason this station is not included in the analysis. For the pairs of stations of interest (linked by red paths in Figure 1) the one-hour data records are correlated to produce a correlogram similar to the one presented in Figure 3(b). The slight velocity changes are observable as slight changes of arrival time of the ballistic waves. Windowed crosscorrelation (or moving-window cross spectral method in the frequency domain) and trace stretching are two techniques commonly used to estimate local time shifts in seismic signals [5] . Within this study, we use the stretching technique [2] to measure seismic velocity variations between pairs of seismic stations. This technique requires the definition of a reference correlation that is computed by averaging all of the 1 h correlations over the duration of the experiment. In a second step, a time window is selected within the lapse time range of interest (where the strongest seismic phases arrived) to compute the velocity variations. For each hour, the time vector t of the correlation is multiplied by a stretching factor ε, i.e., t ← t × (1 ± ε) to align the arrival time of the seismic phases within the selected time window to their arrival time of reference. Many different values of the stretching factor ε are tested. The one retained is the one that maximizes the zero-lag crosscorrelation coefficient between the reference and the velocity variation of each hour. The optimal stretching factor ε is either positive (i.e., the seismic phases arrive later than their reference time) or negative (i.e., the seismic phases arrive before their reference time). Finally, the seismic velocity change curve is defined as the opposite of the vector of the hourly stretching factor.

2.2.2. Piezometric Changes

The velocity changes retrieved following the processing described in the previous section are relative to the reference, chosen here as the average of all 1 h correlations. In order to compare the results with the piezometric data, we simply need to remove the average of each piezometric curve provided for the instruments B10, P96, S20 and M25.

3. Results

3.1. Broadband Velocity Changes Versus Piezometric Changes

Figure 5 summarizes the results obtained for 6 different pairs of stations, orga-

Figure 5. Comparison of broadband seismic velocity changes (black curve, in %) and relative pieozmetric changes (blue curve, in m) for 6 different pairs of stations (7 - 8; 6 - 8; 5 - 6; 2 - 8; 2 - 5; 3 - 4). Computations were made in the 3 - 20 Hz frequency band. Note the excellent agrrement between both curves at the different locations. The letters (A, B, C, D, E, H) corresponds to the location of the virtual piezometers given in Figure 1.

nized from west to east. In each panel, the black curve represents the seismic velocity change retrieved during the experiment, and the blue curve represents the relative piezometric changes measured at the closest piezometer. The overall agreement between the curves is acceptable. Note the correspondence of scales in the present study: one meter of piezometric level change corresponds to 1% of seismic velocity change.

However, for some pairs of seismic stations, the details of the seismic velocity change do not match the evolution of the closest piezometer. The reasons for these discrepancies are partly found in the nonstationary noise source distribution. The activation and rest of the different pumps create spurious artifacts in the velocity change curves that are not related to changes in the medium.

3.2. Depth Dependence of the Velocity Changes

The seismic velocity changes presented in Figure 5 are measured on the ballistic waves that correspond to surface Rayleigh waves propagating directly from a station to the other. These Rayleigh waves have the interesting property to be dispersive, i.e. their velocity depends on the frequency [6] . This property is related to the sensitivity kernel of these waves that depend on the frequency: the lower the frequency, the deeper the kernel. By filtering in different frequency bands, it becomes possible to track the depth range of the measured changes. To address the problem of localization of the change in the medium, the crosscorrelation of the seismic records is computed for several frequency bandwidths (5 - 6, 6 - 8, 8 - 10 and 10 - 12 Hz). The broadband correlations are filtered with a Butterworth filter of order 4. The seismic velocity variations derived for limited frequency bandwidths are compared with those obtained from the broadest frequency bandwidth (3 - 20 Hz), which is used as a reference. The comparison is achieved by computing the autocorrelation of the 3 - 20 Hz seismic velocity curve and the crosscorrelation of the 3 - 20 Hz seismic velocity curve with the 5 - 6, 6 - 8, 8 - 10 and 10 - 12 Hz seismic variation curves, respectively. Figure 6 shows the results for the station pair 5 - 6. First, it is important to note that all limited spectral bandwidths exhibit velocity variations. The amplitude of velocity change in those frequency bands is of similar order (i.e. in the range of ±2%). However only the velocity variations computed within the 8 - 10 Hz frequency range match the reference computed over the 3 - 20 Hz band. The 6 - 8 Hz frequency band provides an acceptable fit to the broadband velocity changes, too. The other frequency bands (5 - 6 and 10 - 12 Hz) provide seismic velocity variations that do not match the reference seismic velocity variations and/or are not in phase with the broadband pattern. This implies that the velocity variations measured in the broadband frequency range are borne by a limited spectral content. In the present case of the pair 5 - 6, the frequency band that senses the environmental changes is 6 - 10 Hz.

A similar analysis is performed for the different pairs of seismic stations. Table 1 summarizes the results for the different pairs. The results of the column

Figure 6. Comparison of the autocorrelation of broadband velocity variations (reference in red) with the crosscorrelation of broadband velocity variation and band limited velocity variation (5 - 6 Hz in light gray; 6 - 8 Hz in black; 8 - 10 Hz in blue and 10 - 12 Hz in dark gray). The best match is obtained for the 8 - 10 Hz bandwidth that bears most of the environmental signal, and it is in-phase with the broadband velocity variations (i.e., with the piezometric measurements of P96). The other bandwidths have little in common with the broadband velocity variation and are apparently out of phase with the reference, except the 6 - 8 Hz bandwidth. The retained frequency band for this pair will be 6 - 10 Hz.

Table 1. Depth dependence of seismic velocity changes for the different virtual piezometers.

“best frequency range” are coherent at the exception of the last pair of stations (3 - 4) that shows a large frequency range. This is possibly due to the path between the two stations mostly formed by the eastern dam that limits the basin 5.2.

The relationship that links the frequency bandwidth with the “equivalent depth” is depending on the local shear wave velocity profile with depth. In the context of this study, different active seismic experiments have been performed around the basin 5.2 to retrieve the velocity profile in the first tens of meters. The outcome of these profiles is that the shear wave velocity in the subsurface is quite low but highly variable. For the sake of simplicity the shear wave velocity is considered as constant over the first 5 meters and equal to 300 m・s−1. The depth z is computed according to the relation z = 0.15*c/f, where c is the local shear wave velocity and f is the frequency [7] . This depth corresponds to the maximum of energy of the fundamental mode of the surface Rayleigh waves with the given frequency. With a constant shear wave velocity profile with depth; this maximum of energy forms a good proxy for the sensitivity kernel of the surface waves that indicates the maximum sensitivity of the waves to a change of shear wave velocity. The 6 - 10 Hz bandwidth has a maximum of energy at a depth of 4.5 to 7.5 meters below the surface. This range overestimates the piezometric values of P96 (Figure 4) that range from 2 to 5 meters below the bottom of basin 5.2. We note a similar behavior for the virtual piezometers A, B, C, D. They form a first group with a minimum depth of 4.5 meters obtained for points B and C, i.e. around the piezometer P96. The virtual piezometers E, F, G indicate deeper levels for the water table. Finally the virtual piezometer H forms a group on its own, indicating a water table oscillating between 2.25 and5.6 meters, to compare to the values of piezometer M25 that indicates a small oscillation of the water table between 4 and 5 meters below the bottom of basin 5.2. The local path between the stations 3 and 4 is radically different from the path of the other stations, dominated by the structure of the eastern dam of the basin. This might explain the difficulty to reduce the frequency bandwidth to a few hertz, due to a strong contrast of velocity between the dam and the underground. In a general manner, the depths of water table inferred from seismic measurements overestimate the true values measured by the piezometers. This discrepancy has to be related to the approximation of a unique medium with constant velocity. A more complete description of the velocity profile below basin 5.2 would be required to properly compute the sensitivity kernels and the equivalent depth of the velocity changes. Of course, taking into consideration a more realistic velocity profile might enhance the importance of the higher modes of the surface waves and mislead the interpretation in terms of depth. Similarly, the filling of basin 5.2 imposes a supernumerary layer that modifies the free surface conditions. For these reasons, the simple approach of the maximum energy and the results of equivalent depth presented in Table 1 must be considered as gross approximations.

4. Additional Constraints on the Hydraulic Dome

4.1. 3D Representation of the Hydraulic Dome

For practical reasons the seismic stations were installed around the basin 5.2. The seismic velocity change curves obtained in Figure 5 integrate the information on the changes of the medium below basin 5.2 all along the path between the seismic stations. The virtual piezometers lettered from A to H are conveniently placed at the midpoint of each path. A distance is computed from the location of piezometer B10 situated on the western bank of the basin 5.2. Figure 7 presents a 3D representation of the results obtained in this study. It includes on the same scale the relative piezometric changes (in meters) and the seismic velocity changes (in percentage) as a function of the distance to B10 and as a function of the time of the experiment. As for Figure 5, the overall agreement between the two independent measurements is acceptable. There seems to exist a difference in the behaviors of the virtual piezometers. Points A, B, C, D and E follow more or less the same variations, similar to the one measured by piezometers B10, P96 and S20. The draining stage of the basin is clearly visible in all virtual piezometers although apparently delayed for point A (see section 4.2 for an explanation). The filling stage is also extremely marked for all these virtual piezometers. The fastest filling rate is observed for the virtual piezometer B, located very close to piezometer P96.

On the contrary, points F, G and H exhibit a different behavior than the previous virtual piezometers. For these three virtual piezometers, the draining stage is not pronounced or is not even present in the data. This has to be related to the initial water level conditions of the basin 5.2 at the beginning of the seismic recordings (Figure 4). It is very likely that the hydraulic dome is locally much less expressed (as seen in piezometer M25 during the filling stage) and vanishes quite rapidly, hence the small relative changes of seismic velocity mostly governed by the regional groundwater flow. The filling stage of basin 5.2 turns into a small

Figure 7. 3D representation of relative piezometric changes (thick ribbons) and of relative velocity changes (thin ribbons) as a function of distance (in meters from B10) and of time (in hours). Note the good agreement between the independent measurements. The virtual piezometers help to constrain the shape of the hydraulic dome in areas deprived from piezometers. Points A, B, C, D, and E have a similar behaviour to the piezometers B10 P96 and S20. Points G and H have smaller amplitudes responses. They do not seem affected by the draining of the basin in the first stage and respond only at the end of the experiment when the hydraulic dome is well developed.

relative change of the velocity, more pronounced for point F than for points G and H. Point F appears as the transition between the western side of basin 5.2, characterized by large and rapid variations of the water table and of seismic velocity, and the eastern side of the basin, characterized by limited changes both in water table level and seismic velocity.

4.2. Limitations of the Approach

If the overall agreement between piezometric measurements and seismological measurements is acceptable, some discrepancies appear when it comes to the fine details of the time evolution of the water table. As illustrated in Figure 3(b), the monitoring of the water table based on seismological measurements requires a stable distribution of the sources of noise. This requirement is hardly achieved in the industrial context of the exploitation, where water pumping and injection happens at different times and locations in the 11 hectares of the field. The nonstationary distribution of the sources of noise creates spurious artifacts in the velocity changes. This might become an issue when trying to identify the depth dependence of the velocity change, because most of the industrial noise sources have a frequency content of a few Hertz. Activation of such sources creates strong energetic seismic body waves that are refracted inside the medium and may arrive simultaneously with the surface waves of interest for the monitoring of the water table. An example of this is provided in Figure 3(b).

5. Conclusion

We have shown here the potential of the seismic noise for the monitoring of the water table in an industrial context. The seismic velocity changes presented in this study integrates the information along the path between the pairs of stations, i.e. horizontally. They also integrate information vertically because of the finite frequency band of analysis that scans the medium at different depths that include fully and partially saturated horizons. However, the velocity change is largely governed by the position of the saturation front (i.e. the water table): an analysis in different frequency bands helps to constrain the depth of this front through time. The nonstationary source of noise distribution impacts the retrieval of the details of the water table evolution. The interpretation of the seismic velocity changes might be impaired by spurious artifacts. Nonetheless, the seismic noise monitoring allows retrieving the major features of the water table evolution. The lateral resolution of the technique is limited by the distance between the seismic stations, with a minimum in this study of 30 meters. The spurious artifacts limit the vertical resolution of the retrieved changes.

Acknowledgements

The authors acknowledge the financial support of Labex OSUG@2020. The authors thank the research platform of “Eau du Grand Lyon” for authorizing entrance on the exploitation field of Crépieux-Charmy and BENEFICIARIO COLFUTURO 2017.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1] Loizeau, S. (2013) Amélioration de la compréhension des fonctionnements hydrodynamiques du champ captant de Crépieux-Charmy. Ph.D. Thesis, Université Grenoble-Alpes, Grenoble.
[2] Sens-Schonfelder, C. and Wegler, U. (2006) Passive Image Interferometry and Seasonal Variations of Seismic Velocities at Merapi Volcano, Indonesia. Geophysical Research Letters, 33, Article ID: L21302.
https://doi.org/10.1029/2006GL027797
[3] Voisin, C., Garambois, S., Massey, C. and Brossier, R. (2016) Seismic Noise Monitoring of the Water Table in a Deep-Seated, Slow-Moving Landslide. Interpretation, 4, SJ67-SJ76.
https://doi.org/10.1190/INT-2016-0010.1
[4] Zhan, Z., Tsai, V.C. and Clayton, R.W. (2013) Spurious Velocity Changes Caused by Temporal Variations in Ambient Noise Frequency Content. Geophysical Journal International, 194, 1574-1581.
https://doi.org/10.1093/gji/ggt170
[5] Mikesell, T.D., Malcolm, A.E., Yang, D. and Haney, M.M. (2015) A Comparison of Methods to Estimate Seismic Phase Delays: Numerical Examples for Coda Wave Interferometry. Geophysical Journal International, 202, 347-360.
https://doi.org/10.1093/gji/ggv138
[6] Beaty, K.S. and Schmitt, D.R. (2003) Repeatability of Multimode Rayleigh-Wave Dispersion Studies. Geophysics, 68, 782.
https://doi.org/10.1190/1.1581031
[7] Graff, K.F. (1975) Wave Motion in Elastic Solids. Ohio State University Press, Columbus.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.