Geophysical methods help to assess potential groundwater extraction sites

Complex geology in glacial terrain makes groundwater exploration challenging for water supply companies. Abrupt lateral changes of geological conditions, for example, raise the risk of unsuccessful exploration drilling. At some locations, local water supply companies experience a large number of dry wells, raising the need for better exploration approaches. The use of geophysical techniques enables the detection and characterization of subsurface structures in terms of aquifers and aquitards. We present two examples in which seismic and resistivity methods are combined for subsurface characterization. A thrust structure is detected in the Varde study area in Denmark and a channel structure is mapped in the Hamburg-Sülldorf study area in Germany. The different imaging characteristics of seismic P‑ and S‑waves are demonstrated. The combination of seismic and resistivity methods enables delineating groundwater-bearing and groundwater-barrier layers down to a depth of 150 m.


Introduction
In areas with complex geology and limited knowledge of subsurface structures, the search for suitable well locations can be a challenging task. For sustainable management of groundwater resources, aquifers of reasonable size, protected by low permeability layers (e.g., clay or till) that prevent contaminant migration from the surface, are required. In glacially-influenced sedimentary areas, common in the North European Lowlands, the layered structure of aquifers and protective layers are disturbed by, e.g., glaciotectonics or deeply incised glacial valleys (e.g., Aber and Ber 2007;Jørgensen et al. 2009Jørgensen et al. , 2012Huuse and Lykke-Andersen 2000;Stackebrandt 2009). Interpolation of drilling results can therefore be difficult or impossible even over short distances, increasing the risk of unsuccessful exploratory drilling. For example, in the Hamburg-Baursberg/Sülldorf groundwater development area, experience shows that the probability of locating a suitable well site by exploration drilling is currently only 42% (Skowronek personal communication).
Geophysical methods for near-surface and hydrogeophysical applications have been improved significantly over the last few decades (e.g., Butler 2005;Rubin and Hubbard 2005;Pellerin et al. 2009;Binley et al. 2015). While resistivity methods are widely accepted for delineating aquifers and aquitards or groundwater salinity, seismic methods are only used sporadically in the exploration and development of groundwater resources. This is probably due to the high costs, ambiguity in the results and interpretative difficulties. Large-scale investigations can often only be carried out in connection with research projects (e.g., BURVAL Working Group 2009). Convincing and meaningful studies concerning application of S-and/or P-wave seismic reflections in northern Germany, Denmark and Canada were presented by, e.g., Polom et al. (2013), Krawczyk et al. (2012Krawczyk et al. ( , 2013, Pugin et al. (2009Pugin et al. ( , 2013; in combination with airborne electromagnetic surveys by Burschil et al. (2012a, b), Jørgensen et al. (2012), Høyer et al. (2013), Tamiru and Wiederhold (2020); and in combination with direct-push by Hoffmann et al. (2008). Recently Sun et al. (2020) demonstrated the imaging of shallow groundwater resources in glacial deposits by high frequency data using an explosive seismic source.
The primary hydrogeologic questions concern the detection of protective covering layers and the minimization of the exploration risk by understanding the subsurface structures and complexity. In this context the objective of our paper is to communicate the combination of seismic and resistivity methods and the advantage of geophysics to address hydrogeologic questions. We demonstrate in this paper the application of geophysical techniques (seismic and resistivity) to define the subsurface architecture in terms of permeable and low-permeable layers for two locations, one in Jutland (Denmark) and one at the northern boundary of Hamburg (Germany) (Fig. 1). P-waves are used at both locations complemented by a shear wave survey at the location Hamburg-Sülldorf. For 2D resistivity imaging we use electrical resistivity tomography and a newly developed mobile transient electromagnetic system (tTEM) (Auken et al. 2019). The added value of our work is the delineation of a buried valley structure in the Hamburg-Sülldorf area and of a thrust structure in the Varde area.

Study areas
The two selected study areas are located in glacially formed landscapes in western Denmark and in northern Germany (Fig. 1). Both areas have been overridden by glaciers during the Elsterian and Saalian glaciation but the Last Glacial Maximum (Weichselian glaciation) reached only a line eastward of our study areas ( Fig. 1a) (Ehlers et al. 2011;Houmark-Nielsen 2011). A good overview on a glaciotectonic thrust-fault complex in the North Sea is given by Winsemann et al. (2020) whereas González et al. (2021) show the internal structure of a push moraine in the Elbe-Weser region close to the Hamburg-Sülldorf location.
The Danish study site, Varde, belongs to the Varde hillislands where deeply rooted glaciotectonism is recognized. The hill-islands contrast the flat outwash plains as elevated, gently undulating landscapes often consisting of till. The elevated terrains are remnants of an old, pre-Weichselian landscape altered by fluvial, periglacial and postglacial processes (Høyer et al. 2013).
The German study site, Hamburg-Sülldorf, is characterized by an alternating, partly glacially compressed sequence of Quaternary sands, gravels, silts, clays and till. The Qua- ternary sequence, which is between 80 and 120 m thick, is underlain by more than 120 m of Miocene mica clay. Up to 150 and 350 m deep Quaternary channels with variable channel filling are known southwest-and southeastward of the area, respectively (Skowronek, personal communication).

Methods
In our study, we used resistivity methods, as they offer a layer characterization in terms of sandy and clayey material (if there is no saline groundwater), combined with seismic methods, as they provide information about the layer structure of the subsoil.
The method of geoelectrical resistivity imaging is widely known as electrical resistivity tomography (ERT). The field layout consists of a four-electrode arrangement, by which electrical current is injected between two electrodes, and the potential difference is measured between another pair of electrodes. Such measurements are repeated for different combinations of electrodes, thus 2D data coverage along a profile is obtained. ERT may be applied at a wide range of spatial scales to determine hydrogeological properties but also to determine states, such as volumetric water content or solute concentration (Ferré et al. 2005).
Electromagnetic measurements also provide information on the resistivity structure similar to ERT, but use a different technique (e.g., Christiansen et al. 2009). An electromagnetic field (AC or DC) is produced by an antenna and induces a secondary EM field in the subsurface. Measurements in the frequency domain (e.g., the Slingram approach) use the comparison of the AC primary and secondary electromagnetic field to derive the resistivity distribution of the subsurface, while measurements in the time domain use the rate of change of the secondary field after abruptly turning off the DC current in the transmitter coil (TEM, transient electromagnetics).
A new development in the field of time-domain EM is the tTEM-system, using a towed, ground-based, transient electromagnetic system, designed for detailed 3D geophysical and geological mapping of the shallow subsurface (0-80 m) in a fast and cost-efficient way (Auken et al. 2019;Maurya et al. 2020;Sandersen et al. 2021). An all-terrain vehicle carries the instrumentation and tows the transmitter frame and the receiver coil which are mounted on sleds for a smooth ride over rough terrain. Detailed 3D coverage is obtained by fast data acquisition (working speed is 15-20 km/h) and close line spacing of typically 10-25 m. This results in a dense lateral data spacing with an in-line distance between the soundings of 5-10 m. Such a high resolution cannot be achieved with electromagnetic airborne systems (e.g., SkyTEM), which are often favored for a 3D coverage of large areas (e.g., Siemon et al. 2009).

Seismic reflection method
Seismic reflection techniques have been well known and accepted for oil and gas exploration for about 100 years. Near-surface applications started with refraction techniques but since 1980 significant progress has been made in seismic reflection surveying (Steeples 2005). In the beginning, P-waves (compressional waves) were mainly created by explosives or weight drops, but nowadays seismic vibrators of different sizes are at hand and can be used for P-wave as well as for S-wave (shear wave) surveys (e.g., Krawczyk et al. 2012Krawczyk et al. , 2013Polom et al. 2013;Tamiru and Wiederhold 2020). These studies also show that seismic reflection methods work quite well in the unconsolidated sediments of the North European Lowlands despite only minor velocity contrasts between layers.
Seismic methods are sensitive to the velocity and density structure of the subsurface. The field layout is based on the transmission of seismic signals from an energy source and the measurement of the signals at receivers (geophones) at different distances. P-and S-waves differ in the particle mo-tion of the wave propagation when travelling through earth material. Due to their lower propagation velocity, S-waves have a shorter wavelength than P-waves leading to a better resolution of small subsurface structures given the same frequency content.
The subsoil material of our study areas mainly consists of sand, till and clay. Although these materials seem to be quite similar, in general surprisingly good seismic reflections at layer boundaries occur. The main condition for reflected seismic waves at a layer boundary is that the layers must have different impedance values (the product of density and seismic velocity). Density and seismic velocity depend on porosity, both decrease with increasing porosity. In general, sand, clay and till have different porosities and thus different densities and seismic velocities. The porosity of sand is 20-40% depending on the grain-size distribution. The porosity of clay can be up to 70% due to the high content of immobile water (low effective porosity). The porosity of till (mainly composed of sand and clay) depends on the clay content: if the volumetric clay content is low (up to ca. 30%), then the till porosity is lower than the porosity of sand (Marion et al. 1992), otherwise the porosity will approach the high porosity of clay. Usually, the clay content of till is low, so density and seismic velocity higher than the seismic velocity of sand can be expected, leading to good reflection results.
The field equipment for the P-and S-wave surveys in this study, and examples of seismic wave fields, are shown in Fig. 2. In the Varde survey, the P-wave Vibrator HVP30 (developed by the Leibniz Institute of Applied Geophysics (LIAG)) was used as the seismic source, while the seismic waves were recorded by an array of planted geophones. In the Hamburg-Sülldorf survey P-and S-vibrators ELVIS (originally developed by LIAG) were used with the geophones for P-wave and a SH-landstreamer (SH indicates horizontally polarized shear waves) for the S-wave data recording (Table 1).

Approach and results
The exploratory approach that we are pursuing to assess potential groundwater extraction sites is based on overview measurements with P-waves, possibly refinement with S-wave measurements (see case 2) and application of resistivity methods to determine the aquifer properties.

Case 1: Varde (glacial thrust structure)
The Varde study site is close to the town of Varde, located in the southwestern part of Jutland, Denmark (Fig. 1). The area is prone to a lack of clean groundwater. Nitrate and pesticides often contaminate the shallow aquifers and it has been a challenge in recent years to identify new uncontaminated aquifers. Recent investigations show that the area is glacially heavily deformed with thrusting and folding of subsurface units. This means that pathways for contaminant migration along these structures to the aquifers frequently exist in the area.
In the course of geophysical exploration, a seismic reflection survey was carried out in the Varde area by LIAG in 2016. The acquisition parameters are listed in Table 1. A tTEM survey, carried out by the HydroGeophysics Group of Aarhus University, followed in 2019. Further exploration included a SkyTEM survey and a pumping test. The first with a focus on deeper structures and the latter to check hydraulic properties.
The seismic section shows structural information down to a depth of about 1000 m. At depths below about -100 m (meters below m.s. l.) the seismic reflectors are more or less horizontal and continuous (Fig. 3b). With reference to Høyer et al. (2013), whose survey area is about 20 km further north, we interpret the deep reflectors as the top of the Arnum Formation (Miocene), the top of the Paleogene deposits and the top of the Cretaceous deposits (Fig. 3b). At depths above -100 m, the reflectors show stronger topo-graphic variation, which presumably indicate a glaciotectonic deformation.
The bulged structure in the northern part of seismic line 2 (starting at about profile meter 1000) might be interpreted as a thrust fault. From the appearance of the reflections in both seismic sections, the direction of stress seems to be from the NNE. It is not clear whether the buckling of reflector "A" in continuation to the south results from glaciotectonic or erosional forces. Unfortunately, the decollement or thrust plane is not clearly seen (reflector "B" in Fig. 3c).
The seismic velocities increase with depth, but there are only very slight lateral variations in the glacially affected part (Table 2). In the northern part the velocities are slightly higher. The interval velocities calculated from stacking velocities are relatively low in the first layer (above reflector "A") ranging from 1300 to 1500 m/s. Below reflector "A" the interval velocity values are significantly higher than 1500 m/s ( Table 2).
The resistivities gained from the tTEM survey show information down to nearly -120 m (below m.s. l.) (Fig. 3a). The maximum depth of reliable data-supported information is about 60 m. Two ranges of resistivities are clear: the high resistive top layer and the lower resistive bottom layer. The latter rises to the north. The layer boundary (marked with "A") correlates quite well with the prominent reflector "A" in the seismic section (Fig. 3a, c). The high-resistive top layer is interpreted as a buried valley with a sandy infill. The resistivity depth section reveals that no low-resistive (i.e. protective) covering layer is present.

Case 2: Hamburg-Sülldorf (glacial valley)
Similar to the Varde area, new well fields are required in the metropolitan area of Hamburg, to ensure a reliable water supply. The targeted aquifer is fine-grained and has low yields. In addition, due to the location in an area of pushed moraines, tilted to steeply standing compressed layers and possibly embedded mica clay lenses are to be expected. Accordingly, there is a considerable exploration risk and a significant risk of total failure in this area (i.e. of finding no productive groundwater zones). With the use of geophysical techniques, a substantial reduction of the exploration risk is expected.
In June 2019, two seismic profiles with a P-wave seismic source (ELVIS P) were carried out by LIAG in the Hamburg-Sülldorf area. A follow-up survey with S-waves (ELVIS S) was carried out in 2020 with the aim of achieving a higher resolution of the layered structure. For acquisition parameters, see Table 1. Several ERT lines are available from a service company.
The difference between the P-and S-surveys is clearly shown in the recorded seismic wave field (Fig. 2c, d, h, i). In the P-wave survey, the traces of the geophones close to the seismic source are dominated by noise from the source, ground roll (surface wave energy) and the air blast (Fig. 2c). No reflected signals are visible in this zone. This means that the minimum depth of layers which can be imaged by P-seismic is about 30 m. In the S-wave survey, early (=shallow) reflected signals are clearly seen close to the source position (Fig. 2h). The reason is that the noise travelling in the longitudinal direction and the air blast are not recorded by the SH-geophones which are more sensitive to horizontal oscillations. Working on paved ground leads to suppression of the surface waves. Even shallow layers can therefore be detected by S-wave surveys. P-and S-waves also differ in their signal frequency as shown in the power spectra. In our example, for P-waves most seismic energy is in the 50-150 Hz frequency range, for S-waves it is between 30 and 100 Hz.
The P-wave seismic reflection profiles show a two-part picture (Fig. 4c, d). The upper part, down to a depth of about 100 m below m.s. l., is characterized by chaotic layering, while the lower part is characterized by a uniform, parallel reflection pattern. This continuous layering is presumably attributable to Neogene or Miocene strata. Further up, horizontal reflections end at the presumed edge of an incised channel (e.g., line 1 at the north side and line 2 at the east side). The channel-like structures are marked in violet. Line 1 runs more or less perpendicular to the strike of the channel. The channel seems to have different phases, the younger one shows asymmetric filling. Strong horizontal reflections inside the valley (marked in blue) might indicate thick till layers as known, e.g., from borehole R99 (for location see Fig. 1). The shear wave section shows much less depth penetration but represents the till layer in more detail. The base reflection of the till is especially well displayed (blue line in Fig. 4b at about -50 m relative to m.s. l.), and the S-wave section shows additional information in the upper 20 m. Here a small additional valley structure in the distance range from 120 to 200 m is a good example of the near-surface resolution that is possible with shear waves (black and white color in Fig. 4b).
The seismic velocities derived from stacking velocities are not robust enough for a material interpretation. Here the electrical resistivity is a much more meaningful parameter; the data are displayed in Fig. 4a, d and Table 2. For better correlation, the valley as interpreted in the P-wave section and a prominent deeper reflection, are also shown in the resistivity section. The deep reflection at about -120 m correlates quite well with the low-resistivity layer (< 20 Ωm) that represents clay and thus is interpreted as the top of the Miocene mica clay. Reflections above and outside the valley might be interpreted as Pleistocene and probably Pliocene sands (yellow reflectors). Inside the channel, resistivities higher than 50 Ωm indicate freshwater-bearing sands, a possible aquifer. Lower resistive regions above, in the depth range down to elevations of -50 m, correlate well with the strong reflections that were interpreted as glacial till. The resistivities support this interpretation. Further above, a thin sand layer is well protected by an overlying clay layer. The very high resistivities on top represent the unsaturated zone and the layer boundary correlates well with the refractor depth indicating groundwater level.

Discussion
The investigations add valuable knowledge on the subsurface conditions that could not be obtained in this way from scattered boreholes. The area north of Varde was, based on the geophysical investigations (seismic, tTEM and additional SkyTEM), interpreted to be affected by rather deep glaciotectonics. This is in accordance with results of a test drilling and a pumping test, which revealed an irregular aquifer within thrusted layers resulting in surface pollution entering the aquifer. Due to these conclusions, this area was not found feasible for a future water supply. Instead, the target area was shifted towards a buried valley system northwest of Varde, where a more homogeneous sand aquifer within the valley was mapped by the SkyTEM survey and a new well field within this valley is now being developed.
In the Hamburg-Sülldorf area, the valley structure is recognized in line 1 by a sub-horizontal reflection pattern (Fig. 4c, d). Line 2 probably runs at an acute angle to the strike of the valley and the northern flank might be observed in the data. In the S-wave section the till layer and especially its base are well-displayed at about -50 m m.s. l. (Fig. 4b).
Deeper reflections are weak. The fact that the strong reflection at about -50 m m.s. l. is weak or interrupted along the distance ranging from 150 to 250 m might be caused by the small depression shown above that scattered the seismic energy and shielded the mapping of deeper horizons. S-waves are more sensitive to scattering due to their shorter wavelength. However, the resistivities also show a gap within the presumed till layer with resistivity values between 20 and 30 Ωm. In the P-wave section, corresponding reflections are disturbed. Both indicate an interruption of the glacial till with the possibility of reduced aquifer protection.
Seismic velocities are better resolved for S-waves than for P-waves due to the lower velocities and longer travel times of S-waves. Therefore, P-wave velocities show only little variation, while the S-wave velocity analysis indicates higher interval velocities for the glacial till. This is supported by logging results of the exploration borehole ABBAU-2 located about 1000 m southwest of observation well R99 (Fig. 1). An upper till layer and a lower silt/clay (mica clay) layer, both embedded in a sandy environment, are clearly indicated by higher gamma activities due to their clay content. The density log shows increased densities for the till layer. However, no significant density contrast between sand and interlayered silt/clay is observed (Fig. 5a). For comparison, we show in Fig. 5b seismic velocities for sand, clay and till layers from VSP (vertical seismic profiling) measurements. The measurements were carried out on the island of Föhr (Burschil et al. 2012b) and in the glacial valley Ellerbeker Rinne (Rumpel et al. 2005) about 25 km north of the Hamburg-Sülldorf area. As expected, P-wave velocities for till layers are clearly higher than velocities for sand layers. With one exception (circle and triangle at 40 m depth around 500 m/s in Fig. 5b) this is also valid for S-wave velocities. As already indicated in the density log, significant velocity differences between sand and clay are not observed in these VSP data (Fig. 5b).
The resistivity values are much more meaningful with respect to the aquifer properties than the seismic velocities. Due to the surface electrical conductivity of clay minerals, even a minor clay content of the material reduces the electrical resistivity significantly ( Table 2).
The various data were processed independently. The relatively good match of the depth of the Miocene mica clay in the P-wave seismic and resistivity data suggests trustworthy processing. Slight discrepancies of resistivity and

Conclusion
We showed the application of geophysical methods to locate potential groundwater extraction sites for two areas and presented a qualitative interpretation. The geological structure was revealed in both areas by the geophysical investigations. In the Varde area a sandy channel filling could be separated from a glaciotectonic thrust structure. Test drilling and a pumping test revealed an irregular aquifer with only a limited protective cover which would result in a considerable vulnerability of the aquifer due to contamination from the surface. In Hamburg-Sülldorf a buried valley structure was identified, the existence of which was not formerly known. Neogene and Quaternary facies were distinguished, and sedimentary or erosional features are now better understood, a prerequisite for a successful drilling campaign. The differentiation of the top of the Miocene clay, Neogene sands, and a Quaternary channel with internal structure was possible with the presented geophysical investigations.
In both areas the combination of seismic depth sections with resistivity data has greatly increased knowledge on the subsurface architecture. At the Varde site, for example, only the geophysical surveys could provide the information that was required for rejecting the site for water supply. For the Hamburg site, only the geophysical survey could provide the required information which was necessary to understand the subsurface architecture. Here, too, the findings are guiding the development of the new well field and the site of the next exploration drilling. The approach using a P-wave survey for the first overview, a refinement by a Swave survey and additionally the use of resistivity methods for addressing the material properties was successful and is recommended for future seismic surveys in comparable areas.
Acknowledgements The TOPSOIL project is co-funded by the Interreg North Sea Region Programme (J No 38 2 27 15). The manuscript was significantly improved by the comments and questions from an unknown reviewer and C. Juhlin.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.