High‐Resolution Geological Information from Crosshole Ground Penetrating Radar in Clayey Tills

Heterogeneous glacial deposits dominate large parts of the Northern Hemisphere. In these landscapes, high‐resolution characterization of the geology is crucial for understanding contaminant transport. Geological information is mostly obtained from multiple boreholes drilled during a site investigation, but such point‐based data alone do not always provide the required resolution to map small‐scale heterogeneity between boreholes. Crosshole ground penetrating radar (GPR) is suggested as a tool for adding credible geological information between boreholes at contaminated site investigations in industrial sites where infrastructure, such as electrical installations, can pose a challenge to other geophysical methods. GPR data are sensitive to the dielectric permittivity and the bulk electrical conductivity, which can be related to the distribution of water content and sand/clay occurrences. Here we present a detailed crosshole GPR dataset collected at an industrial contaminated site in a clay till setting. The data are processed using a novel inversion approach where information on changes in the velocity and attenuation of the radar signal are obtained independently. The GPR results are compared to borehole logs, grain size analyses, and relative permeability data from the site. The GPR data analysis provided valuable information on the understanding of the lateral geological variability. A silt layer with a thickness of a few decimeters, likely important for flow characterization, was confirmed and resolved by GPR data. Our findings suggest that crosshole GPR has the potential for contributing with high‐resolution geological information by filling the data gap between boreholes, thereby becoming a relevant tool in contaminated site investigations.


Introduction
Contaminated sites worldwide pose a risk to human health by threatening water resources, indoor environment and through direct exposure to contaminated soil.The European Union estimates approximately 2.5 million potentially contaminated sites need to be identified, investigated, and managed.In large parts of the Northern Hemisphere, the geology is characterized by heterogeneous glacial deposits (McCabe 1987;Shaw 1987;Houmark-Nielsen 2010), mainly clayey tills riddled with sand lenses of various sizes (Kessler et al. 2012;Klint et al. 2013).Heterogeneities, such as fractures and sand lenses, affect contaminant distribution and transport in these clayey tills, due to the contrast in hydraulic conductivity and biogeochemical properties (e.g., McKay et al. 1993;Jørgensen et al. 2002;Harrington et al. 2007;Damgaard et al. 2013).High content of sand in the clayey till matrix also increases the hydraulic conductivity (Aamand et al. 2022).Mapping the extent and geometry of these sandy features, as well as the presence of continuous thick clay layers, is critical for understanding contaminant flow paths in the clayey till.
Information at contaminated sites is often limited to the locations of the boreholes drilled as part of a site investigation.From such boreholes, geological logs are obtained as well as information from soil and water samples.Additionally, direct push probes provide logs on, for example, bulk electrical conductivity and relative permeability (McCall et al. 2014;McCall and Christy 2020;Steelman et al. 2020).Common for all these data types is, that they provide depthdiscrete information about the vertical variability of the subsurface, but little information about the lateral variability, which in practice is most often interpolated from these point data to build the conceptual geological model.
Geophysical methods can be used to obtain information on the subsurface variability between boreholes and have been applied to support conceptual site model development at contaminated sites.Høyer et al. (2019) used surface electromagnetic induction (EMI) data, direct current (DC) resistivity data, and transient electromagnetic (TEM) at a landfill situated in a clay till dominated setting.The interpretation of the geophysical data provided valuable information that was used in the development of a highly detailed 3D geological model for the area.In Maurya et al. (2018), results from an induced polarization survey in a major contaminant plume were used to support the geological modeling in the area.For these surface methods, the vertical resolution inevitably decreases with depth, and the penetration depth of, for example, surface ground penetrating radar (GPR) is limited in clayey till.Crosshole methods allow constant vertical resolution along the borehole depth, given appropriate spaced boreholes (Day-Lewis et al. 2005).Commonly used crosshole methods are seismic tomography (e.g., Wang and Rao 2006), electrical resistivity tomography (ERT) (e.g., Slater et al. 2000), and GPR tomography (e.g., Eppstein and Dougherty 1998).
Crosshole GPR investigations use electromagnetic (EM) pulses that are emitted from a transmitting antenna in one borehole and recorded by a receiving antenna in a nearby borehole.The recorded radar waveform carries a fingerprint of the medium through which the wave has traveled, depending on the spatial distribution of the subsurface dielectrical properties.The EM wave velocity can be related to the dielectric permittivity of the subsurface, which can be translated to subsurface water content (e.g., Huisman et al. 2003).The attenuation of the EM wave provides additional information on the bulk electrical conductivity, which is different for sand and clay (e.g., Reynolds 2011), and crosshole GPR hence has the potential to distinguish between sandy and clayey features between boreholes (Looms et al. 2018;Jensen et al. 2022a).GPR is minimally-invasive and can be used in already existing boreholes, and can therefore be integrated in any ongoing site investigation.Crosshole GPR signals are generally thought not to be adversely affected by electrical installations at the site, which allows for the use of GPR investigations at urban and industrial sites.Furthermore, the fast data acquisition potentially allow for practical implementation at operational scale.Finally, the sensitivity of the GPR method is highest in the middle between the boreholes (Day-Lewis et al. 2005), and thus contributes with information in areas of low data density, which is an advantage compared to standard ERT methods where the resolution is highest close to the boreholes.
Crosshole GPR has been used for a wide range of applications, such as monitoring infiltration using time-lapse tracer studies (e.g., Hubbard et al. 2001;Chang et al. 2006;Gueting et al. 2017), characterizing subsurface moisture content (e.g., Alumbaugh et al. 2002), and for characterizing porosities in aquifers (e.g., Klotzsche et al. 2013;Zhou et al. 2020).Until recently, GPR was not applied in clay tills due to the attenuation of the signal, but Looms et al. (2018) and Svendsen et al. (2023) showed that GPR can be used in these environments, given sufficiently short borehole separation up to 5 to 7 m.Furthermore, Jensen et al. (2022a), recently developed a linear inversion algorithm that allows for fast interpretation of both GPR traveltime and amplitude data, which makes the method attractive for use as part of contaminated site investigations.
The scale of GPR investigations is limited, and the method is hence most relevant in small-scale investigations of, for example, a source zone at sites with highly heterogeneous geology, such as flow till (Klint et al. 2013), where interpolation between boreholes meters apart can be challenging.In a subglacially weakly deformed flow till, small-scale sand lenses with lengths of 0.2 to 4 m and thicknesses of 2 to 50 cm have been identified (Kessler et al. 2012) and shown to be important for contaminant transport (Kessler 2012).Hydraulically important fractures most often have apertures smaller than what can be resolved using linear inversion of GPR data, hence the application of GPR to map hydraulically important heterogeneities is most relevant in layers where sand lenses act as preferential transport pathways as seen in, for example, Harrington et al. (2007) or Damgaard et al. (2013).
In this study, we present the first application of the linear inversion approach (Jensen et al. 2022a) to crosshole GPR traveltime and amplitude data at an industrial contaminated site situated in a clay till environment.The main aim is to assess whether the application and analysis of GPR data can provide supplementary high-resolution information on subsurface variations and delineate sandy structures, hence improve the understanding of the geology at the site at a small scale.Secondly, we identify and discuss the methods' limitations and advantages when used in an industrial contaminated site setting.

Field Site-Contamination and Geological Setting
Field investigations were conducted at an industrial site in Vassingerød, Denmark.The shallow geology is dominated by Quaternary deposits.The glacial sediments are composed of 30 to 35 m of sandy meltwater deposits, which is the main aquifer in the area, overlain by 10 to 15 m of clay till with embedded sand lenses.The clay till consists of two till units: a lower, consolidated unit (basal) and an upper, less consolidated till unit (flow or melt-out till) (NIRAS 2019).A transect from the geological model for the entire site is seen in Figure 1.
The field site has previously been thoroughly investigated in order to understand the spreading of an identified TCE contamination from the source zone in the clayey till to the underlying sand aquifer and to estimate the contaminant mass discharge (NIRAS 2019;Rosenberg et al. 2022Rosenberg et al. , 2023)).The GPR investigation in this study took place in the suspected source zone area in the clayey till.Further details of the site and contamination can be found in Rosenberg et al. (2022Rosenberg et al. ( , 2023)).

Boreholes, Investigations, and Data Acquisitions Boreholes and Data
In this study, we use the available data within the local investigation area set in the source zone (NIRAS 2019;Rosenberg et al. 2022) to make a more detailed local geological model (Figure 1, top).Specific borehole logs and grain size analyses associated with boreholes B301, B302, and B303 were included in this study, while relative permeability data were obtained from MiHPT14, MiHPT15, and MiHPT16.The locations are shown in Figure 1.Crosshole GPR data were acquired between boreholes B302 and B303, since the largest lateral geological variability was expected between these two boreholes.The investigated tomographic transect is illustrated in Figure 1 (top-right).We include information from B301 and MiHPT14-16 in the following analysis to better understand the overall lateral variability at the site.Approximate locations of boreholes B301-B303 and relative permeability logs have been projected onto the geological profile in Figure 1.

Data Acquisition
The relative permeability data were collected using the Membrane Interface and Hydraulic Profiling Tool (MiHPT), but only the HPT data will be used in this study as a measure for the relative permeability of the subsurface (McCall et al. 2014;McCall and Christy 2020).Further details can be found in Data S1 and Rosenberg et al. (2022).
The grain size data is based on soil samples taken during the drilling of B301 to B303, every 0.5 m from 3 m below ground.As described by Lévy et al. (2022), the grain size distribution was analyzed using a sieve (0.063 to 2 mm) and a laser diffractometer (0.02 to 63 μm).The soil classification from the grain size samples follows the USDA particle size limit classification (e.g., Gee and Bauder 1986).
Crosshole GPR data were acquired using the PulseEKKO system from Sensors&Software with 100 MHz omnidirectional antennae and a sampling frequency of 0.2 ns.We used a 1000 V transmitter and the Ultra Receiver.The antenna housing dimensions are Ø29 mm and 1.44 m long.The GPR equipment is shown in Figure 2. The boreholes were drilled as a part of the contaminated site investigation in 2019 and were equipped with plastic (PEH) casing tubes with an inner diameter of d = 0.051 m.The distance between B302 and B303 was 3.81 m.Since the casing tubes were screened at 13.0 m below ground, a specially designed packer system was installed in the boreholes to remove water from the casing tube before the GPR investigations.This was necessary, to prevent water from affecting the propagating EM wave, as investigated in Jensen et al. (2022b).
Zero-offset-profile (ZOP) data were collected by employing the transmitter in borehole B302 and the receiver in B303, then moving the transmitting and receiving antenna simultaneously in vertical increments of 0.25 m from 1.5 to 12.0 m below-ground level.The data acquisition geometry is illustrated in Figure 2, where the lines indicate the straight wave paths for the propagating EM signal between a given transmitter and receiver position.The ZOP data are used to obtain a 1D profile of water content and amplitude variations along the borehole depths.The ZOP acquisition time was 15 min.
Multi-offset-gather (MOG) data were collected by fixing the transmitting antenna at one depth in B302 while lowering the receiving antenna in vertical increments of 0.25 m along the borehole depth in B303.The data acquisition geometry for MOG measurements is illustrated in Figure 2. The receiver positions were chosen so the collected data span ±45° from horizontal relative to the transmitter position, to avoid recording high-angle traces as they would be discarded in the post-processing (Peterson 2001).The transmitter was then moved 1 m down and the procedure was repeated.When the transmitter had covered the entire borehole depth, the receiver was moved in increments of 1 m, while the transmitter was moved in increments of 0.25 m.The MOG data are used in a geophysical inversion to estimate 2D subsurface models of water content and bulk electrical conductivity variations along the borehole depth and in the area between the boreholes.The MOG acquisition time was approximately 3 h.
Calibration gathers in air were acquired before and after recording the ZOP data as well as before, between and after each half MOG dataset, to calibrate the traveltime data for the zero-time of emitting the radar pulse, following the procedure described in Oberröhrmann et al. (2013).The field campaign took 1 d, which included installation and testing of the packer system, calibration gathers, ZOP and MOG data collection, and was carried out in June 2022.

Post-Processing of GPR Data
The traveltime and amplitude data were identified in the recorded waveform data and picked using a semiautomated approach.The first-arrival traveltime data were picked using cross-correlation between a reference wavelet and the acquired waveform data (Molyneux and Schmitt 1999;Hansen et al. 2013).The amplitude data were picked as the peak-to-peak waveform amplitude value in a defined window.No gain was applied to the data.The critically refracted EM wave traveling along the ground affected the shallower traces, and the first usable trace was obtained 2.75 m below ground.Prior to inversion, the measured MOG amplitude data were corrected for geometrical spreading of the energy, radiation patterns of the antenna and the antenna gain effect, A 0 .See, for example, Holliger et al. (2001) and Jensen et al. (2022a) for details on amplitude corrections.Traveltime data were time-zero corrected using the time-zero correction factor, t 0 , estimated from the acquired calibration lines.We obtained a MOG dataset of 464 traces in total.We also investigated using a subset of this full dataset, omitting all data from transmitter and receiver positions between 6.5 and 8.5 m below ground, resulting in a reduced dataset of 232 traces.The data geometry for the reduced dataset can be found in Data S1.

Inversion
The corrected MOG data were used in a tomographic geophysical inversion to obtain estimates of the 2D subsurface parameter distribution.The measured amplitude data were inverted to obtain tomograms for subsurface attenuation and the measured traveltime data were inverted to obtain the subsurface slowness (inverse velocity) distribution (as described in e.g., Giroux et al. 2007).We use a probabilistic linear inversion with a geostatistical Gaussian prior model included and we account for the forward modeling errors arising from choosing a linear scheme.The method is described in detail by Jensen et al. (2022a).The obtained tomograms presented here are the mean of all the equally likely subsurface representations that are consistent with the a priori information and the data.These mean results are overly smooth compared to each subsurface representation and the parameter value range is dampened (Tarantola 2005;Hansen et al. 2008;Looms et al. 2010).Further information on inversion parameters can be found in Data S1.
The obtained slowness distribution, s, from traveltime inversion can be related to the water content, v , through the dielectric permittivity, r and a petrophysical relation, for example, the empirial Topp's equation (Topp et al. 1980): where the permittivity can be estimated from r = s 2 c 2 , and c is the speed of light in vacuum (c = 0.2998 m ∕ ns).
The obtained attenuation distribution, , from amplitude inversion can be related to the electrical conductivity, , through (e.g., Annan 2005): where the dielectric permittivity, r , is estimated from the traveltime inversion as described above and Z 0 is the free space impedance (≈ 377 Ω) given by the free space magnetic permeability, 0 = 1.25 × 10 −6 H/m, and the permittivity of free space, 0 = 8.89 × 10 −12 F/m.

Existing Information
We construct a geological model for the area between the boreholes B302-B303 from an interpretation of the existing geological information in the area of investigation (Figure 3, right).Overall, we interpret two till units, consistent with the large-scale geological model for the area (Figure 1), separated by a silt layer at approximately 36 masl.
The silt layer is interpreted as coherent, however thinning toward B302.This interpretation is mainly based on borehole logs as the discrete grain size analyses from B302 has two silt peaks and it is difficult to verify the thinning from the sampling distance of 0.5 m intervals.The low required pressure observed between 37 and 35 masl in HPT15 and at 36 to 37 in HPT16 is indicative of a more permeable layer, such as silt, but this is not identified in HPT14.Furthermore, the thickness of the silt layer appears overestimated in the HPT15 log, compared to the borehole and grain size logs.This could be a result of sparse fracture networks in the transition between the upper and the lower till units, separated by the silt layer.
Below the silt layer, the variance between the data types (HPT logs, borehole, and grain size logs) from B301, B302, and B303, is low, and the till unit seems more uniform laterally.The lower till unit is perhaps divided into an upper and lower part, where the upper subsection above 34 masl has a higher sand content, as identified in the grain size logs, and the lower is more clayey.The subdivision of the lower till unit into a sandier upper and a more clayey lower part is not observable in the HPT logs or indicated in the borehole logs.A transition zone toward the sandy aquifer with an increasing amount of sand from approximately 32.5 masl, is apparent in all data types.A slight dip of the top of the transition zone toward B303, is interpreted from grainsize logs, and is also evident in the borehole logs.
The data from the upper till unit is inconclusive.The grain size log from B303 suggests a more sandy till in the entire depth interval 40 to 36 masl, which is contradictory to HPT data that indicate a geology with higher required pressure (lower relative permeability), that is, higher clay content, at least in the depth interval 37 to 38 masl.The borehole log for B303 suggests an interval of more sandy till between 37.2 to 38.8 masl.
The grain size log from B302 suggests an interval with higher clay and silt content at approximately 37.5 masl, which is supported by the high pressures recorded in the HPT logs in this depth interval, and a more sandy till above 38 masl.However, the borehole log does not reflect this variation in geology.In order to synthesize a geological model for the upper till, we mainly base our interpretation on the grain size data, since it can be challenging to evaluate small changes in, for example, sand content on site from a visual and mechanical inspection of a sample.
The interpolation between B302 and B303 is challenged by the differences in the data, which results in a highly uncertain conceptualized geological model in the depth interval 36.5 to 40 masl (Figure 3, right).We have included small heterogeneities in the interpretation of the upper till to honor this uncertainty as well as the grain size data.

1D Profiles (Zero-Offset-Profiling Data)
1D profiles obtained from zero-offset-profile acquisition of GPR data between boreholes B302 and B303 are presented in Figure 4.The recorded waveform data is equivalent to what can be seen in the field during data acquisition, and contains information on (1) traveltime (x axis), hence radar propagation velocity that can be translated to subsurface water content, as well as (2) signal strength identified as waveform amplitude (colorscale), hence attenuation of the radar wave, which provides additional information on subsurface bulk electrical conductivity.
The EM wave velocity, calculated by dividing the ray path length with the time-zero corrected traveltime data, and the picked waveform amplitude is shown in Figure 4b and  4c, respectively.
The transition from the unsaturated to the saturated zone is seen as a substantial decrease in EM velocity (i.e., an increase in water content) between 40.5 and 39.5 masl (Figure 4b).The exact location of the water table in the formation is difficult to determine from the GPR data due to the retention properties of the formation that create a capillary fringe zone.The silt layer at approximately 36 masl clearly impacts the velocity data, which is low in the layer, due to the high water-saturated porosity of the layer.The silt layer appears to act as a low-velocity waveguide, which traps the propagating EM wave.This is confirmed when inspecting the raw waveform data.Waveguide effects are seen as high-amplitude arrivals later in the signal in the waveguide depth interval (marked in Figure 4a), as well as destructive signal interference for receiver positions close to both the upper and lower boundaries of the waveguide, similar to observations in Klotzsche et al. (2012).This interference arises from different wavefronts arriving simultaneously at the receiver.These depths are indicated with asterisks in Figure 4c.
The observed transition zone to the sandy aquifer is identified from the increase in signal amplitude below 33 masl, that is, lower attenuation hence lower electrical conductivity, which is interpreted as a decrease in clay content.The velocity data is only affected to a lesser degree due to lack of variation in water content.The 1D GPR data support the original geological model in Figure 3 suggesting that the lower till unit can be subdivided into an upper more sandy interval and a lower, more clayey interval.This is interpreted from a section of higher amplitudes from just above 35 masl down to 34 masl, indicating lower electrical conductivity, and lower amplitudes around 34 to 33 masl, indicating higher electrical conductivity.The 1D GPR data does not inform on the uncertain interpretation of the upper till unit in the original geological model.We observe intermediate GPR amplitudes around 4 mV at 37 masl, indicating a lower clay content, but a predominant layer of sandier clay till is not interpreted at this depth in the original model.The high amplitudes above 39 masl could be caused by either an increase in sand content, a decrease in water filled porosity or both.
The reason that the GPR method does not show abrupt changes is due to the sensitivity kernel of the propagating EM wave, that is, the propagating radar wave is affected by a volume of at least 0.5 m to all sides along the travel path of the radar wave, causing a smoothing of the recorded signal.

2D Tomograms (Multi-Offset-Gather Data)
The results from the inversion of MOG traveltime and amplitude data collected at the field site between boreholes B302 and B303 are shown in Figure 5a and 5b.The water content tomogram in Figure 5a and 5c is estimated from the obtained slowness tomogram (not shown) and the electrical conductivity (EC) tomogram in Figure 5b and 5d is estimated from the obtained attenuation tomogram (not shown).
We observe diagonal streaking in the EC tomogram throughout the model area, most noticeable around 34 to 36 masl as indicated by the arrow in Figure 5b.These artifacts as well as the observed high contrasts in radar wave velocity caused by the silt layer, suggests that the linear inversion approach is challenged due to the presence of the low-velocity waveguide.
To avoid waveguide effects caused by the silt layer, we repeated the inversion with a reduced dataset.The water content and EC tomograms obtained from inversion of the reduced dataset are shown in Figure 5c and 5d, respectively.The area with no data coverage is cross-hatched and the variations in the tomograms at these depths are only a consequence of the chosen a priori information included in the inversion.The diagonal streaking is greatly reduced compared to the full dataset tomogram and the structures recovered with the reduced dataset now correspond to the structures in the water content tomogram.The mean value of the water content tomogram is slightly shifted due to the large amount of data that are removed, but the recovered structures are consistent with the full dataset tomogram in Figure 5a.
We base our geological interpretation of the GPR results on the tomograms shown in Figure 5a and 5d.Overall, the vertical variations captured in the water content tomogram, Figure 5a, correspond well with the variations observed in the 1D velocity profile in Figure 4b.The transition from unsaturated to saturated conditions is again observed from the water content model between 40.5 and 39.5 masl, and the silt layer is recovered and clearly delineated as a layer with high water content at 36 masl.The estimated water content ranges between approximately 0.25 and 0.35 in the clay till, and 0.40 in the silt layer.In the saturated zone, the water content translates directly to porosity.The added information of the 2D water content tomograms, is primarily the thinning of the silt layer toward B302 that supports the observations in the borehole logs and the suggestion of a layer with higher water content, possibly connected, at 38.0 to 38.5 masl in B302 and at 38.5 to 39.0 masl in B303.
From the EC tomogram, Figure 5d, the subdivision of the lower till unit is confirmed, with a coherent layer of higher clay content (higher conductivity) at 34.3 to 33.3 masl and higher sand content (lower conductivity) above 34.3masl.The transition zone to the sandy aquifer is not as evident as in the ZOP results (Figure 4).Above the silt layer, the high EC values above 20 mS/m at 38.0 to 38.5 masl in B302 and at 38.5 to 39.0 masl in B303, suggest a possibly connected layer.Apart from this, we do not see clear structure in the top till unit.

Updated Local Geological Model
An updated geological model between boreholes B302 and B303 was constructed based on the interpretation of the 1D GPR profiles as well as the 2D tomograms for water content and electrical conductivity.The updated model is seen alongside the original interpretation in Figure 6.The silt layer coherency between the boreholes B302 and B303 and thinning toward B302, is confirmed from inversion of GPR traveltime data.The subdivision of the lower till unit is supported by the GPR amplitude data with a sandier clay till above 34.3masl and a coherent layer of clayey till with higher clay content between 34.3 and 33.3 masl.The original interpretation of the transition zone below 32.5 masl is confirmed mainly by the GPR ZOP results.
In the upper till, the GPR results and the grain size data from the boreholes do not concur, and we therefore evaluate that these heterogeneities are not large enough to be resolved by GPR data.The size of these lenses has, therefore, been reduced in the updated geological model.
The only exception to the above is the layer observed at approximately 38 to 39 masl that is visible in both the water content and the EC tomogram.Whether the increase in EC is caused by either an increase in porosity, clay content, or both, is difficult to determine from the GPR data.However, the depth of this layer coincides with observations of higher clay/silt content in the grain size analyses from both boreholes.We therefore interpret a possibly connected clayey layer at this depth.

Discussion
The Use of Crosshole GPR to Refine the Geological Model We used the results from the GPR data analysis to update the original interpretation of the geology and connect the borehole information as the GPR data provided valuable information on the subsurface geology between wells and thus allowed for a robust interpretation of the lateral connectivity.This study illustrated the added value of obtaining both a water content and an electrical conductivity tomogram.The two-dimensional water content model provided essential information about the subsurface porosity, and the obtained porosity values for the clay till are in line with reported values for clay till in Denmark (Aamand et al. 2022).The water content tomogram aided the delineation of the silt layer due to its high porosity, which arises from the well-sorted material.The water content tomogram did however not inform substantially on the geological variations of the remaining subsurface due to the fairly constant porosity in these areas.In contrast, the conductivity tomogram showed larger variability in these areas, which we interpreted as a result of variation in the clay content.Using this information, we were able to identify a continuous section with high clay content in the lower till unit, as well as sections less likely to be continuous in the upper till unit.
We did, however, observe challenges with the simple inversion routine, as the high contrast in water content (i.e., in EM velocity) was not included in the prior description and therefore caused artifacts in the conductivity tomogram (Hansen et al. 2014).To solve this, we omitted the data most affected by waveguide effects in the inversion.The subsurface area with no data coverage was therefore not resolved in the conductivity tomogram, instead, the water content tomogram informed on the geology in this depth interval.The consistency between the full and reduced water content tomograms implies that the linear inversion of traveltime data is robust, even in the presence of pronounced velocity contrasts in the subsurface.Waveform data with late-arriving and high-amplitude elongated wave trains occur when a lowvelocity layer is embedded in a geology with higher EM velocity (Klotzsche et al. 2012) and the waveguide effects therefore contain important information about the geology.Klotzsche et al. (2013) presented an amplitude analysis approach where the presence of waveguides as well as their upper and lower boundaries can be identified from direct analysis of the recorded waveform amplitudes, i.e. without the need for inversion of the data.This approach could be employed in this and future studies when lowvelocity waveguides are suspected.
The investigated area between boreholes B302 and B303 is quite small compared to the area of the field site.This study is intended as a pilot study, and for future applications, the possibility of larger field scale GPR investigations should be explored.For this particular site, the lateral extension of the silt layer could be investigated across the site since it is easily detected in the ZOP data, given that boreholes with appropriate distances exist (up to 4 to 8 m depending on the geology).The depth of investigation is limited by borehole depth and the GPR cable length, which is 30 m for the equipment used in this study.
Another cross-borehole geophysical method, frequently used for environmental application, is electrical resistivity tomography (ERT), often used in combination with induced polarization (IP).Similarly to the crosshole GPR method, crosshole resistivity methods have been applied for various purposes such as monitoring of in-situ remediation (e.g., Mao et al. 2016;Bording et al. 2021;Lévy et al. 2022) and mapping of possible contaminant flow paths in clay till (Bording et al. 2019).As mentioned previously, the resolution of resistivity tomography is, in contrast to GPR tomography, highest close to the boreholes (Day-Lewis et al. 2005).Crosshole resistivity and radar methods can therefore be used in combination since they provide complementary information for inversion and hydrogeological interpretation (e.g., Linde et al. 2006;Winship et al. 2006;Looms et al. 2008;Haarder et al. 2012;Cassiani et al. 2014).The electrodes used for ERT methods must be installed in specially made boreholes, and high contact resistance in the boreholes can negatively affect the data quality (Bording et al. 2019).While the initial installation of electrodes may be expensive and time-consuming, ERT data can be obtained almost autonomously after electrode installation.This allow for long-term monitoring campaigns as, for example, Lévy et al. (2022), where 10 data acquisition rounds were carried out over 2 years to monitor spreading of a remediation agent.Contrarily, crosshole GPR can be used in existing investigation wells, which make GPR attractive as a low-cost method to obtain additional information in ongoing site investigations.

Antennae Separation
One limitation to consider when applying GPR at a contaminated site is the distance between the boreholes, since the signal is attenuated as it travels through the formation.In this study we measured a sufficient signal-tonoise level at a borehole distance of 3.81 m.Even for high angle traces in the MOG gather where the traveled distance is at least 5.3 m, the signal was sufficient.We estimate that in this type of geology, the boreholes could have been up to 6 m apart for ZOP investigations and 4.5 m for MOG investigations.

Water in Wells
The presence of water in the screened wells enforced the need for installing a packer system before conducting the GPR measurements.Dry boreholes are important for obtaining high-resolution tomograms in the saturated zone in clay till settings as the water downshifts the frequency content of the traveling EM wave and therefore reduces the possible resolution (Jensen et al. 2022b).The boreholes were pumped dry using a packer system that was made especially for this field campaign using off-the-shelf products and deployment was a bit tedious and time consuming.We consider the development of a packer system that is easy to install in existing, screened wells as an important step toward field scale use of GPR at contaminated sites.Alternatively, an additional sealed inner tube could be inserted in the well given that the inner diameter of the well casing allows it.The antenna housing diameter is Ø29 mm and the minimum recommended borehole diameter is 33 mm.Despite this limitation, the fact that GPR can be deployed in existing investigation wells (with non-metallic casing), is a clear advantage.

On-Site Evaluation of Data
An advantage of the method when applied in the context of contaminated site investigations is the fast data acquisition and the possibility to obtain information about the subsurface between boreholes directly in the field, by visual inspection during data collection of one-dimensional ZOP data.Subsequent preliminary data processing of the 1D ZOP data can be done in 10 to 15 min.

Equipment
The GPR equipment is compact and easily transported.It fits into the back of a standard SUV, and can be set up in 10 to 15 min (see https://youtu.be/VJDFd JBLvJ0).The equipment can even be used indoors if required, for example, within a storage facility.

Reproducibility
Lastly, the GPR data are consistent when measurements are repeated.We evaluated the reproducibility of GPR data by plotting overlapping data acquired from (1) the first and (2) the second half of the MOG data acquisition against (3) the acquired ZOP data.The overlapping data were recorded approximately 1.5 h apart.The data fell along a straight line, as expected, and a linear model was fitted to the traveltime and amplitude data, respectively.For the traveltime data fit R 2 = 0.994 and for the amplitude data fit R 2 = 0.993.This confirms that the acquired GPR data are reproducible and are not arbitrarily affected by industrial installations or by small deviations in antenna position and orientation.

Conclusions
Crosshole GPR measurements were carried out at a contaminated site in a clay till setting to evaluate whether the GPR data analysis can provide valuable information about the geology between two selected boreholes located in a source zone.The GPR data analysis provided a strong guide for interpretation and validation of the existing data types (HPT data, grain size logs, borehole logs), and the analysis confirmed presumptions about the continuity of structures between boreholes in the lower till unit.The GPR data analysis resulted in a better understanding of the geological heterogeneity in the upper till.
The original interpretation of the geology between the boreholes was updated on the basis of the GPR results and the uncertainty of the conceptual geological model in our local investigation area was reduced, however not eliminated.The GPR investigation conducted at this site provides information that could be used to design a more targeted uncertainty analysis of the hydraulic properties of the two till units and their impact on transport pathways.
With this pilot study, we show that the use of crosshole GPR at contaminated sites with complex geology can provide high-resolution information about the geology between the boreholes.The fast data acquisition and the novel efficient data processing make GPR investigations attractive to use as a low-cost additional source of information in ongoing contaminated sites investigations.Furthermore, GPR can be used at industrial and urban sites where infrastructure, such as pipes and electrical installations, can pose a challenge to other geophysical methods.
A pilot study like this could be extended by upscaling of the GPR method and the results could, as a logical next step, be related to contaminant transport at clayey till sites.

Figure 1 .
Figure 1.Upper: Overview of the field site in Vassingerød, Northern Zealand, Denmark (approximate location on map).The local area of investigation is shown in zoomed view (top) with the colored points that are considered in this study.The tomographic transect investigated by GPR is illustrated as the green dashed line.Lower: Geological profile for transect A-A′.The transect illustrates the general geology at the site, with sand lenses embedded in the clay till, underlain by a sandy aquifer.Transect adapted from NIRAS (2019).The approximate locations for MiHPT logs (red) and investigation wells (blue) are projected onto the transect.

Figure 2 .
Figure 2. Left: Photo of the field site with the GPR equipment.Antennas are laid out on ground during warm up of the equipment.Middle: Schematic illustration of GPR transmitter and receiver boxes connected to the 100 MHz antenna https://www.sensoft.ca/products/pulseekkopro/borehole-gpr/ (with permission from Sensors&Software).Right: Data acquisition geometry of ZOP data and MOG data.Asterisks indicate transmitter positions and circles indicate receiver positions.The lines indicate the straight wave paths for the propagating EM signal between a given transmitter and receiver position.

Figure 3 .
Figure 3. Existing information: Borehole logs for B301 to B303, grain size analyses for B301 to B303 and HPT logs from MiHPT14-16, see Figure 1 for locations.Interpreted geological model between boreholes B302 and B303 from existing data (right).

Figure 4 .
Figure 4. Zero-offset-profiling (ZOP) data collected between B302 and B303: (a) Recorded waveform data (radargram).Green asterisk indicates high amplitude arrivals later in the signal.(b) EM velocity.(c) Amplitude data, obtained from the raw waveform data.The orange asterisks marks destructive signal interference at both the upper and lower boundary of the silt layer.The ZOP data can be considered a 1D model, representing the average subsurface between the transmitter and receiver locations.

Figure 6 .
Figure 6.Original geological model (left) and updated conceptual geological model constructed on the basis of knowledge obtained from GPR data (right).

Figure 5 .
Figure 5. Two-dimensional tomograms obtained from inversion of Multi-Offset-Gather (MOG) data.Estimated water content tomogram and estimated electrical conductivity (EC) tomogram.Tomograms (a) and (b) are results from full dataset inversion.An area with diagonal streaking artifacts is indicated with the red arrow.Tomograms (b) and (c) are results from reduced dataset inversion where data between 34.9 and 36.9 masl are omitted to avoid influence from the silt layer.