Characterisation of geological thin layer by borehole radar

Abstract

To investigate a thin layer, we conducted borehole radar measurements at the GPR facility of Tohoku University, Japan. Both zero-offset profiling and tomography measurements were used to obtain water content distributions. Analyzing first arrival times with amplitudes, we found that a borehole radar can detect velocity changes caused by a high-water content, which acts as an electromagnetic waveguide. We conducted a 3D FDTD simulation and confirmed that the EM wave propagates as a guided wave in a thin layer related to antenna frequency. We determined a thin geological layer compared to an antenna frequency that cannot be detected by reflection measurement due to the poor radar resolution. By comparing the simulated data with measurement data, we could determine the geological boundary between host layers and a thin layer with high water content. In contrast, both the amplitude and the travel time of direct waves are affected by antenna positions and electrical parameters, including permittivity and conductivity. A vertical velocity profile was obtained by combining complementary zero-offset profiling and tomography data, which can be used to investigate water content distribution. The low-velocity thin layer was characterised by amplitude attenuation and late arrival time.

KEYWORDS:

• Borehole radar
• velocity analysis
• attenuation
• electromagnetic wave
• FDTD

Introduction

Ground-penetrating radar (GPR) can help to interpret the sedimentary environment by providing information on the layer geometry and electromagnetic properties. Borehole radar has the potential to detect water-containing structures (Sato and Thierbach 1991). Thin layers (layers that are not resolvable in terms of the applied wavelengths) are common in sedimentary deposits, especially in a hydrocarbon formation. In hydrogeology, groundwater flow and contaminant transport are considered thin-layer structures. Alluvial sediments consist of heterogeneous sedimentary deposits with a wide range of structures and textures (Kolterman and Gorelick 1996). A layer with high water content compared to a host layer and clay lenses are required for modelling groundwater flow and the solute transport process. Arcone (1984) first examined the electromagnetic properties of a thin ice sheet overlying water. As the ice–air and ice–water boundaries are strong reflectors, a thin ice sheet can be an electromagnetic waveguide as same as that we consider a horizontal layer between boreholes. van der Kruk, Streich, and Grenn (2006) numerically analyzed the high-permittivity sand layer as a waveguide structure. Layers, such as clay lenses, and preferential flow paths within aquifers can be of a limited thickness (Klotzsche et al. 2014). Thus, these layers can be a waveguide structure regarding electromagnetic wave propagation travelled along thin structures. In seismic method, guided seismic waves resulted from interference within the waveguide are used to predict the continuity and discontinuity of reservoir structures between wells (Perra et al. 2002). In this research, attempts are made to better understand the nature of radar responses to thin layers and analyze in quantitatively. Huo, J et al. (2021) described a guided wave that travels along the borehole in a high conductivity environment. Crosshole radar can be used to delineate subsurface geology in various environments. Klotzsche et al. (2012, 2013) introduced an amplitude analysis approach to detecting and identifying horizontal wave-guiding structures. Among the geophysical techniques, cross borehole GPR can provide an accurate estimation of 2D and 3D moisture content fields (Alumbaugh and Chang 2002). Crosshole radar techniques are based on the first arrival time of a direct wave between a transmitter and a receiver antenna placed in different boreholes. Direct, refracted, and reflected waves are recorded, depending on the electrical properties of the medium, borehole condition, and measurement configuration. Dense ray coverage achieved by changing the antenna position can define the layer structures. Amplitude information can also assist in evaluating the medium between the boreholes. Due to the amplitude loss of an EM wave, EM wave attenuation can determine the formation. Ray-tracing techniques are applied to find ray paths for a subsurface model (Matsuoka and Ezaka 1992). The initial purpose of crosshole radar is to obtain the most detailed in situ velocity profile of radar waves for specific investigations. The velocity of the radar signal and amplitude are strongly affected by the amount of pore water. Thus, the travel time of a radar signal and amplitude can assist to estimate sediment layers.

Conventionally, hydrogeological parameters are obtained by drilling techniques, including logging tools, core samples, and pumping tests. In recent decades, geophysical methods such as seismography, electrical resistivity, and GPR have been widely used and have characterised hydrogeological properties (Hubbard and Rubin 2000; Irving, Knight, and Holliger 2009; Mangel et al. 2012). GPR method provides values for conductivity and permittivity, which affect the velocity and attenuation of EM waves (Klotzsche et al. 2014). Due to the high contrast between permittivity in air $({\epsilon }_r=1)$ and in water $({\epsilon }_r=80)$ (Davis and Annan 1989), permittivity is strongly dependent on soil water content. Dielectric permittivity of soil is highly sensitive to the water content, and weakly sensitive to soil particles and types. This phenomenon permits that soil water can indicate dielectric permittivity of soil, regardless of soil types. Thus, GPR measurement intense pay an attention for water content, especially in soil layers.

Crosshole radar is well suited to derive high-resolution images and recharge processes of aquifer systems (al Hagrey and Muller 2000; Alumbaugh and Chang 2002). Near-surface sediment is often layered. This layering can cause radar waves to propagate as guided waves (Ellefsen 1999) for crosshole radar measurement. The high-contrast thin layer caused by an increase in water content acts as a low-velocity waveguide layer when the thickness of the layer is smaller than the in situ dominant wavelengths of the GPR signal (Arcone, Peapples, and Liu 2003). The EM waves that travel within a thin layer may propagate over a large distance. The presence of waveguide crosshole measurements was studied for seismic data (Franssens, Lagasse, and Mason 1985). The phenomenon of the EM wave that occurs within a thin layer as a guided wave is complicated and not well studied.

This paper aims to demonstrate the interpretation of borehole radar data, along with a finite-difference time-domain (FDTD) simulation. To understand how a thin layer affects electromagnetic waves, a ray-tracing technique was applied to represent a thin layer with relatively high-water content hosted by lower water content layers. An electromagnetic guided wave only propagates the satisfied condition, which depends on applied EM wavelengths, layer thickness, and borehole conditions. First, we acquired crosshole data, both zero-offset profiling (ZOP) and tomography mode and converted it to the velocity profiles between the boreholes, assuming straight ray paths. Afterward, synthetic studies were examined based on the velocity profile. To characterise a layer structure observed in borehole data, we focused on working with a thin layer structure, EM wave propagation within the layer, and its boundary between layers. Even in high-attenuation (high water content) soil layers, we could observe direct waves in the crosshole measuremenrt. Therefore, we attempted to understand EM behaviour with a thin layer as a guided wave and to determine structural information.

Electromagnetic wave propagation between boreholes

Due to the applied EM frequency and the electrical properties of the material (conductivity $\sigma$, permittivity $\epsilon$, and magnetic permeability $\mu )$, EM waves cannot be observed over a large distance compared to the applied frequency in a highly conductive medium. By considering the soil layer as a homogeneous medium, the complex wavenumber $k=\beta -j\alpha$ (1) represents the attenuation coefficient $\alpha$ and the phase coefficient $\beta$ respectively, given by (1): (1) $\begin{array}{r}k=\beta -j\alpha =\omega \sqrt{\mu \epsilon }{\left[1-j\frac{\sigma }{\omega \epsilon }\right]}^{\frac{1}{2}},\end{array}$(1) In explicit expression, attenuation coefficient $\alpha$ and the phase coefficient $\beta$ can be written, respectively, as following: (2) $\begin{array}{rl}\alpha & ={\left[\frac{{\omega }^2\mu \epsilon }{2}\left\{{\left(1+\frac{{\sigma }^2}{{\omega }^2{\epsilon }^2}\right)}^{\frac{1}{2}}-1\right\}\right]}^{\frac{1}{2}},\end{array}$(2) (3) $\begin{array}{rl}\beta & ={\left[\frac{{\omega }^2\mu \epsilon }{2}\left\{{\left(1+\frac{{\sigma }^2}{{\omega }^2{\epsilon }^2}\right)}^{\frac{1}{2}}+1\right\}\right]}^{\frac{1}{2}},\end{array}$(3) where $\omega$ is the angular frequency.

In most GPR applications, an attenuation depends on frequency, and phase distortion can be neglected due to the electrical properties. Several wave phenomenons can be acquired between a transmitter (Tx) and a receiver (Rx) antenna, including direct waves, reflected waves, and refracted waves. By assuming a straight ray path, we calculated the velocity between boreholes at each depth. Therefore, the first arrival time indicates the fastest ray path among all the ray paths. We used straight propagation models to describe a velocity model between the boreholes. The derived velocity reflects the average velocity of all the different ray paths travelled with different velocities, as shown in Figure 1. The velocity of travel time in different layers is determined as follows (4): (4) $\begin{array}{r}{v}_{avg}=\frac{c}{\sqrt{{\epsilon }_{avg}}}.\end{array}$(4)

Figure 1. Schematic illustration of ray paths of crosshole measurement.

The shortest ray path between the boreholes is determined by Snell’s law, as described below (5): (5) $\begin{array}{r}\frac{\sin {\theta }_1}{\sin {\theta }_2}=\frac{v_1}{v_2},\end{array}$(5) where ${\theta }_1$ and ${\theta }_2$ are the incident and refracted angles, respectively, as shown in Figure 1.

Methodology

Crosshole measurement

Crosshole technique allows the velocity of EM wave propagation, to be obtained along direct path between antennas. Therefore, the existence zone of velocity anomalies associated with water content can be extracted. By analyzing each trace between antennas, we can be able to analyze subsurface layered structures. In crosshole measurement, Tx and Rx antennas are placed in different boreholes, and the positions of both antennas are changed to measure EM transmitted waves. Our experiment was performed at the GPR facility of Tohoku University, Japan. The facility has two testing boreholes 4.5 m apart, as shown in Table 1.

Table 1. Borehole dimension.

Parameter
Depth	20.0 m
Diameter (Inner)	100.0 mm
Casing type	Plastic case (PVC)
Water level	8.5 m
Distance between holes	4.5 m
Direction	Vertical (90°)

The data were acquired in April 2019, using a RAMAC 100 MHz slim borehole antenna. Tx antenna was set in Borehole 1 (BH1) at depths of 4–18 m, and the Rx antenna was set in Borehole 2 (BH2) in multi-offset measurement. To conduct zero-offset profiling (ZOP) survey, Tx and Rx antennas were first placed at the bottoms of the boreholes. Both antennas were then pulled up simultaneously in 0.1 m steps. ZOP data were acquired from the interval between 2 and 19 m. By the arrival time of recorded data and the known distance (4.5 m) between the boreholes, the velocity profile of the radar wave was determined by assuming a straight ray path.

Assuming a homogeneous medium and straight ray path, we can obtain the velocity profile along a borehole, as shown in Figure 2(a). By calculating a velocity with constant distance, we determined permittivity and converted it into volumetric water content (VWC) using Topp’s equation (Topp, Davis, and Annan 1980). However, the arrival time cannot evaluate subsurface structure due to the inhomogeneous subsurface and non-straight ray path.

Figure 2. (a) GPR profile of the ZOP. The black dashed line indicates the calculated velocity and the water content. (b) Modelled structure with assumed parameters. (c) Simulated GPR profile of ZOP mode.

We consider that the velocity analysis can be enhanced using an amplitude. We combined arrival times and amplitude information in the crosshole measurement to understand the medium between the boreholes. The velocity data using ZOP measurement roughly distinguished the horizontal layers, as shown in the GPR profile (Figure 2(a)). To define the exact boundary between the layers, tomography (multi-offset) measurements were used to calculate the arrival time of the direct and refracted waves. We applied the ray-tracing technique to identify soil layers distinguished by velocities. We also applied a numerical simulation using the FDTD method to confirm the measurement data. The amplitude changes confirmed that the layer boundaries depended on water content. Finally, we could evaluate subsurface structures in detail and obtain physical parameters such as velocity and permittivity.

Detection of a thin layer with high water content

We used simulated travel times of the direct and refracted waves to evaluate the boundary between layers. We assumed that the ray paths were straight and parallel to estimate the velocity. A thin layer with high water content hosted by a high-velocity layer (11.2 m to 16.2 m) appears in the ZOP data. We performed a 3D FDTD numerical simulation for the modelled structure using CST Microwave Studio software, based on the ZOP data shown in Figure 2(c). Gaussian pulse was used to excite a signal. Mesh density was controlled by the wavelength (six lines per wavelength). The layers were distinguished by their electrical parameters, and the two boreholes were modelled for the crosshole numerical simulation. We designed a half-wavelength dipole antenna with 1.5 m length (100 MHz) and used it for the FDTD simulation. First, the FDTD simulation was conducted as a ZOP measurement. The rapid changing of arrival times shows the boundaries between the layers. Due to the inhomogeneous subsurface layers and borehole conditions, the simulated ZOP data does not fit well with the measured data. However, the simulated arrival time and amplitude variation can indicate the interfaces of layers. A thin layer (high water content) is observed with a late arrival time and low normalised amplitude in Figure 2(c). In the measured ZOP profile, the low amplitude and late arrival time indicate a thin layer with high water content.

Consequently, we obtained tomography measurements with high-density ray coverage, as shown in Figure 3. In tomography measurement, the radar signals are generally unobservable in high-angle ray paths and layers with high water content due to the increased attenuation. However, we could observe the arrival times by increasing the measurement stacking in both ZOP and tomography measurements.

Figure 3. Tomography ray-coverage: (a) Tx: 4 m; (b) Tx: 6 m; (c) Tx: 8 m; (d) Tx: 10 m; (e) Tx: 12 m; (f) Tx: 14 m, and; (g) Tx: 16 m. Colour intensity indicates the apparent velocity.

Figure 3 shows the ray coverage measured in each tomography mode, and colour intensity indicates the calculated velocities (Giroux et al., 2007). The tomography measurement (Tx: 6 m) clearly shows the velocity changes observed in ZOP data. Thus, the interface between the two different layers is recognised in both measurement methods. As we assumed layers with high water content, a thin layer compared to the antenna length (wavelength) appeared at around 14 m. We conducted several simulations to change the thickness of the layers with high water content compared to an antenna length. The important consideration is selecting a layer thickness when investigating waveguide structure. It is of interest that the received GPR signal depends strongly on these parameters, including antenna length and layer thickness. Due to the total internal reflection occurred in a thin layer, the antenna length associated with wavelength is considered.

We used the layer with high water content (${\epsilon }_r=16$) as a guided structure within hosted layers (${\epsilon }_r=8$) for the FDTD simulation. This simulation model could correspond to the saturated soil layer hosted by the layers with lower water content. To investigate the high amplitude and the late arrival spreading, we plot snapshots of the absolute electrical fields for different models, as shown in Figure 4. The thin layer can act as a waveguide structure when its thickness is close to the antenna length. We observed wave propagation behaviour through thin structure. In addition, the thin layer structure can act as a waveguide when internal reflection occurs within a layer. According to Snell’s law, total internal reflection occurs when the incidence angle is higher than the critical angle. The critical angle in the thin layer (high water content) is calculated by the following equation (6): (6) $\begin{array}{r}{\theta }_c=si{n}^{-1}\left(\sqrt{\frac{{\epsilon }_1}{{\epsilon }_2}}\right)\end{array}$(6) where, ${\epsilon }_2$ and ${\epsilon }_1$ are the dielectric permittivities of the thin and host layers, respectively.

Figure 4. Modelled layers between the boreholes. (a) Antenna is placed at the boundary of the layers (${\epsilon }_{r1}=8$, ${\epsilon }_{r2}=16$). (b) Antenna length is placed in a thin layer (${\epsilon }_{r1}=8,{\epsilon }_{r2}=4$, ${\epsilon }_{r3}=8$). (c) Antenna length is placed in a thin layer (${\epsilon }_{r1}=8,{\epsilon }_{r2}=16$, ${\epsilon }_{r3}=10$). The solid black line indicates a layer boundary.

The energy is trapped in the thin layer, as shown in Figure 4(c). Beyond the critical angle (4), multi-reflection occurred within the thin layer. As trapped energy, EM waves can propagate over a large distance. The delayed travel times were observed in the layer with high water content. The amplitude information in the measured data identified the low-velocity thin layer caused by high water content and its boundaries. We found that the contrast of the electrical parameters affects the waveguide. A thin layer can be a waveguide structure only with a satisfied condition that traps the energy. With Tx antenna located in the thin layer and the Rx antenna straddling the thin layer depth, trapped EM waves of high amplitude and late arrival times are recorded.

Figure 5 shows the simulated GPR profiles of ZOP with the different electrical parameters. The guided wave appeared in the thin layer as the high amplitude and the late arrival time in the simulated GPR profile in Figure 5(b). When internal reflections occur in the thin layer, the amplitude can be maximised, which causes high-amplitude waves. Klotzsche et al. (2013) showed Full-waveform inversion result to identify and characterise a low-velocity aquifer layer between the boreholes. In our research study, we only used the ray-tracing technique to clarify wave propagation between near-distanced boreholes. Due to the precisely evaluating boundary of the thin layer, a ray-tracing analysis gives a better understanding of wave propagation compared to the inversion technique.

Figure 5. GPR profiles of the simulated data. (a) Low-contrast layer. (b) High-contrast layer. The solid black line indicates the boundaries of the thin layer.

Figure 6 shows the amplitude profiles of both ZOP and tomography measurements. Due to the measurement distance and attenuation level, the amplitudes change over a wide range. We observed low-amplitude signals within the thin layer in the measured GPR profile, as shown in Figure 6(f). The reason is that the contrast of the electrical parameters could not satisfy the waveguide condition in the measurement data. We could observe the amplitude loss in the thin layer due to the high contrast conductivity. In the case of higher conductivity, no maximum amplitude for the thin layer is observed.

Figure 6. Amplitude profiles of measured data: (a) Tx: 4 m; (b) Tx: 6 m; (c) Tx: 8 m; (d) Tx: 10 m; (e) Tx: 12 m; (f) Tx: 14 m, and; (g) Tx: 16 m.

Figure 7. GPR profiles of the fixed transmitter. The solid black line indicates the simulated travel time of the direct waves. The black dashed line indicates the boundary between the layers. (a) Tx = 4 m. (b) Tx = 12 m. (c) Tx = 14 m.

Figure 8. (a) Measured GPR profile. (b) Simulated GPR profile. (c) Measured signal. (d) Simulated signal. The solid black line indicates the selected ray trace.

In contrast to the lower conductivity, the maximum amplitude is detected within the thin layer as a guided wave. We can see the layers distinguished by their differing amplitude levels in Figure 6. The amplitude intensity distinguishes the horizontal layers. Due to the high angular data in the tomography, we see abrupt changes in the amplitude intensity that used to define the boundary between horizontal layers. The thin layer that appeared in the measured data cannot show a high amplitude signal because of the high conductivity caused by the high water-content. In Figure 6(a), the maximum amplitudes overlap each other. The amplitudes of tomography measurements became lower than those of ZOP measurements due to the increase in antenna separation. As shown in Figure 6(b), the amplitude changes indicated the boundary between layers where Tx antenna was placed at such a boundary. When detecting a thin layer, ZOP measurement could not define the boundaries. Only the minimum amplitude indicated a thin layer with high water content. The amplitude of tomography measurement shows the boundaries of layers, as shown in Figure 6(f).

We tried to fit the GPR profiles of tomography measurement to the simulated travel time which used constant velocity. Using a known antenna distance, we assume a travel time in a homogenous layer. We roughly define a layer boundary if there is a difference between assumed travel time and the real observed arrival time. At 3, 5.1, and 6.2 m, we can see the changes in arrival time in Figure 7(a), which shows the interface between the different layers. The layer velocity can be defined as the difference between arrival time and the simulated arrival time. Thus, the velocities of the layers that started at 3 and 5 m must be lower than the velocity used in the simulated travel time.

On the other hand, the velocity of the layer which starts at 12 m must be higher than 0.063 m/ns, as shown in Figure 7(b). By fitting the simulated data to the measured data, the layers have been characterised by defining velocity and water content. However, the high-contrast thin layer cannot be characterised by the simple analysis that we used. We found that EM wave propagates as a guided wave through the thin layer distinguished by electrical properties, especially dielectric property. To satisfy the waveguide structure, we need to apply an EM wave that can travel by internal reflection. By numerical analysis, it has been confirmed that EM energy is concentrated when internal reflections occur in a layer, as shown in Figure 5(b).

Wave propagation through a boundary

A physically long-length dipole antenna is used for borehole measurement to generate a low frequency for a high penetration. Most current borehole GPR systems use frequencies below 100 MHz, achieving a penetration range of 20–50 m in crystalline rock (Slob, Sato, and Olhoeft 2010). Dipole antenna cannot generate EM wave as a point source. We observed waveform discontinuities in the GPR profile, as shown in Figure 8(a). If the feeding point of the half-wavelength dipole antenna is placed at the boundary of different velocity layers, as shown in Figure 4(a), the waveform discontinuity appears in the GPR profile both for measured and simulated data, as shown in Figure 8. Due to the EM wave travelling through the different layers, the waveform appears as separated waves. Figure 8(c, d) shows the selected trace along the solid black lines. The time delay between the separated waveforms depends on the velocity contrast.

Conclusion

This study investigated EM wave phenomenon in a thin layer caused by high water content. ZOP and tomography measurements were combined to investigate horizontal layers between boreholes. For synthetic data, the thin layer was considered by the wavelength and electrical parameters. In practice, a dipole antenna is used for borehole measurement due to its well-fitting to the borehole. The guided wave can interpret a layered structure by analyzing an arrival time and amplitude information, especially a thin layer. We can detect a thin layer and its boundaries by comparing simulation data with measured data. In ZOP and tomography measurements, the antenna position strongly affects the amplitude and travel times. When antennas are placed within a thin layer, the amplitude directly indicates a waveguide.

In some cases, the antenna length can be higher than the layer thickness. Thus, some of the energy directly radiates into the host layers. As a result, energy concentration is lower than the case the dipole antenna is in a layer. Moreover, the layer thickness can strongly affect the measured signal. We tried to evaluate horizontal layered structures using the amplitude and arrival time. Comparing the applied methods provided the most accurate estimates of the layer boundaries. Crosshole radar could provide a quantitative estimate of the medium between the boreholes. The main objective of the research work was to evaluate horizontal soil layers distinguished by electrical properties. The crosshole radar survey discussed in this paper assisted in determining a thin layer with high water content. These surveys confirm that the borehole radar survey is effective in identifying structural information.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Japan Society for the Promotion of Science [grant number JP20K20990].