Combining crosshole and reflection borehole ground-penetrating radar (GPR) for imaging controlled freezing in shallow aquifers

During test operation of a geological latent heat storage system as a potential option in the context of heat supply for heating and cooling demand, part of a shallow Quaternary glacial aquifer was frozen at the TestUM test site. In order to evaluate the current thermal state in the subsurface, the dimension of the frozen volume has to be known. As the target is too deep for high-resolution imaging from the surface, the use of borehole ground-penetrating radar (GPR) is being investigated. For imaging and monitoring of a vertical freeze–thaw boundary, crosshole zero-offset and reflection borehole GPR measurements are applied. The freezing can be imaged in the zero-offset profiles (ZOPs), but the determination of ice body size is ambiguous because of the lack of velocity information in the frozen sediment. Reflection borehole GPR measurements are able to accurately image the position of the freezing boundary through repeated measurements of ±0.1 m, relying on the velocity information from ZOPs. We have found that the complementary use of ZOPs and reflection measurements provides a fast and simple method to image freezing in geological latent heat storage systems. The presence of superimposed reflections from other observation wells and the low signal-to-noise ratio are problematic. The use in multiple observation wells allows an estimation of ice body size. A velocity model derived from multiple ZOPs enabled us to extrapolate geological information from direct-push-based logging and sediment cores to a refined subsurface model.

1 Introduction

A high proportion of space heating is still based on burning fossil fuels (Steinbach et al., 2020). For the transition to renewable efficient heating, alternative concepts have to be developed and improved. Additionally, the predicted climate warming will result in an increased need for cooling applications (Lozàn et al., 2019). Using heat pumps, latent heat storage (LHS) allows the extraction of high amounts of energy for heating with little temperature change and low storage volume (Agyenim et al., 2010). For phase change materials with a low freezing point, such as water, the residual frozen volume generated during the heating period can act as a cold storage and later be used for cooling purposes. Therefore, LHS can meet the demand for both energy-efficient heating and cooling. In conventional LHS, tanks containing the phase-change material are burrowed, which makes installation expensive and limits building opportunities, especially in densely occupied urban areas, where there is a high heating and cooling demand (Steinbach et al., 2020).

Accessing geological layers, specifically aquifers, as storage volume, instead of burying tanks, can increase application possibilities because there is less above-ground space required and reduced installation costs. With segmented heat exchanger probes that allow depth-oriented controlled freezing and thawing (Fig. 1), the geotechnical use of aquifers as LHS is possible (Dahmke and Schwarzfeld, 2022).

Figure 1 (a) Schematic sketch of geological latent heat storage at the experimental site in Wittstock with the concept of GPR reflection and GPR crosshole measurements. (b) Overview of the experimental site with the positions of borehole heat exchangers and observation wells.

To validate and improve thermo-hydraulic modeling of geological LHS, the current thermal state at operating depth has to be assessed. Data from borehole in situ temperature measurements and cooling fluid feed and return temperatures only give punctual information. Since the boundary between frozen and unfrozen ground indicates a major change in thermal energy, knowing the lateral extent of frozen ground around the heat probes can act as a proxy for the thermal state, so imaging the lateral propagation of freezing with geophysical methods is desirable.

Extensive experience can be referred to in the use of geophysical methods for the observation of thawing and freezing processes in permafrost soils and glaciers (e.g., Campbell et al., 2021; Schwamborn et al., 2002; Terry et al., 2020; Vonder Mühll et al., 2002; Weigand et al., 2020). Electrical resistivity tomography (ERT) has regularly been used for the monitoring of geological storage (Hermans et al., 2014, 2015; Lamert et al., 2012), although, being sensitive to the strong change in electrical conductivity between frozen and unfrozen sediment, the application is not suitable for a precise estimation of a vertical freeze–thaw boundary. The lateral resolution of surface measurements like ERT decreases with greater depths, and the high conductive layer on top of the monitored groundwater aquifer further decreases resolution in target depth. Even with borehole ERT and sharp boundary conditions applied in the inversion, the measurements at an LHS site would be affected by the installed heat exchangers and connecting pipes and cables.

The high contrast in electrical permittivity of ice (ε≈3–4) and water (ε≈81) makes ground-penetrating radar (GPR) very sensitive to the phase change, and it is therefore widely used for permafrost monitoring (e.g., Cao et al., 2017; Du et al., 2020; Hinkel et al., 2001; Sokolov et al., 2020; Steelman et al., 2010; Stephani et al., 2014; Stevens et al., 2008; Wang and Shen, 2019). Using borehole measurements gives the advantage of getting the source and receiver close and perpendicular to the imaging target. Successful application was shown in Kim et al. (2007) by mapping ice rings in an underground liquid gas storage. The downside of this setup is the non-directionality of the antennas used, making determination of the reflector position ambiguous.

So far, the use of GPR has been examined in monitoring freezing and thawing processes within natural geosystems. However, its potential application as an indispensable measurement method for the operational monitoring of geological LHS systems remains unexplored. We hypothesize that borehole GPR will bring additional value in terms of spatial information of the ice body and that it can thus be used for monitoring purposes.

Therefore, the aim of this study is to use borehole reflection and crosshole GPR to image an ice body in a shallow aquifer. Firstly, crosshole GPR is used to get an enhanced site characterization before installation of a geological LHS system, and then the concept of borehole reflection and crosshole GPR measurements is tested to image the ice development during plant operation.

2 Methods

2.1 Experimental site

The experimental site, TestUM, is located on a former airfield near Wittstock/Dosse in the northeast of Germany. The near-subsurface consists of unconsolidated Quaternary glacial sediments, which is a representative scenario for many areas in the North European Plain. The TestUM site was already intensively used for several injection experiments. These prior experiments investigated the effects of injecting heat (Heldt et al., 2021; Keller et al., 2021; Lüders et al., 2021), carbon dioxide (Peter et al., 2012; Lamert et al., 2012), and hydrogen (Hu et al., 2023; Keller et al., 2024; Löffler et al., 2022) into the subsurface. The subsurface investigations for these projects in the area show a high spatial geological variability typical of glacial sedimentation. There are layers of fine-, medium-, and coarse-grained sands; glacial till (boulder clay) characterized by a mixture of sediments, including high clay content; and gravel alternating in thickness and lateral extent. Therefore, prior to establishing the experimental LHS site, a comprehensive site investigation was conducted using an adaptive approach that integrated surface geophysical measurements, direct-push-based investigations, and traditional coring techniques. This multifaceted methodology was employed to identify the most suitable area for the experiment at the site. In particular, surface electrical resistivity and electromagnetic induction measurements provided a non-invasive option to analyze subsurface characteristics in terms of near-surface anthropogenic remains (old foundations, remains of the former airfield). Complementing this, direct-push-based investigations with the hydraulic profiling tool (HPT; Geoprobe, USA) enabled a detailed examination of the geological setting with the help of vertical high-resolution data of hydraulic and electrical conductivity (McCall and Christy, 2020; McCall et al., 2014; Vienken et al., 2012; Dietrich et al., 2008). Additionally, traditional coring was employed to extract core samples, facilitating a more in-depth analysis of the geological strata. An area was selected for the experiment where coring indicated a sandy aquifer at a depth of 10–16 m, covered above and below by clayey till.

To test geological latent heat storage, a subsurface volume is frozen and thawed in a controlled manner (Fig. 1a). A total of 16 heat exchanger probes, each 16 m deep, are installed in a 1 m×1 m grid (Fig. 1b). The probes are designed with two separate fluid circuits, with the upper spanning 3 to 9 m and the lower spanning 10 to 16 m depth. Operating the lower fluid circuit with a feeding temperature <0 °C and the upper circuit with a temperature >0 °C allows depth-oriented controlled freezing in the sandy aquifer while keeping the overlying geological units unfrozen. A total of 18 wells (outer diameter: 60 mm) and 9 multilevel wells are installed to have direct access to the subsurface for GPR measurements, in situ temperature measurements, and probing for hydrochemical and microbiological analysis. The water table varies between z_water=2–3 m b.g.s. throughout the year; thus the experiment takes place in a fully saturated environment.

The 2 in. wells are 18 m deep with mesh filtering at 10 to 16 m depth and are therefore filled with water. Well naming is adapted to groundwater flow direction from northeast to southwest, which flows at a velocity of 0.05–0.09 m d⁻¹ (Heldt et al., 2021). Upstream wells begin with the letter U, central wells being with C, and downstream wells begin with D. Wells for borehole GPR are positioned at a distance of 1.5 m from the outer heat exchangers to be close enough to register reflected signal despite the high attenuation caused by the high electrical conductivity in water-saturated sediments. Closer wells would have a higher probability of being affected by freezing, which would make them inaccessible. Also, if the reflector is too close to the antenna, the direct wave superimposes the reflected signal. Wells are only on three sides of the experimental field because an opening had to be left for logistical reasons, e.g., heat exchanger installation with heavy machinery.

2.2 Measurement methods

2.2.1 Inclinometer

Observation wells can deviate from their designated vertical orientation. The extent of deviation is measured using the DevProbe1 inclinometer (Geotomographie, Germany). Readings of tilt and heading are taken for every borehole, allowing calculation of the well path in the subsurface. The mean of the occurring horizontal deviation at 17 m depth is 0.28 m, with a maximum of 0.99 m at well C05 due to a drilling obstacle. The corrected true well position in the subsurface is used for GPR interpretation.

2.2.2 Ground-penetrating radar

The use of GPR reflections can give specific information about the distance from a well to a reflector, assuming a well-known propagation velocity. Therefore, the use of GPR is favored in this setting. With our target being a vertical-layer boundary at a depth greater than 10 m, surface GPR is not applicable because the overlaying aquitard with high electrical conductivity limits penetration. With the help of boreholes, source and receiver antenna are brought close enough to the imaging target to record signals despite the high attenuation. Two omnidirectional Tubewave-100 (Radarteam, Sweden) borehole antennas, with a vertical offset between the center of the transmitting and receiving array of 0.3 m, in combination with a GSSI SIR-4000 (GSSI Geophysical Survey Systems, Inc., USA), are used for the survey. The peak frequency, measured inside the observation wells at the experimental site, is 190 MHz. We performed three measurements at three different stages of the test site operation: (1) before installation of the heat exchangers for enhanced site characterization, (2) a baseline in the final equipped site, and (3) a measurement after 54 d of running the LHS system when the existence of an ice body is verified by temperature measurements. For stage (1), we measured zero-offset profiles (ZOPs). To image the lateral extent of a body of frozen saturated sediment in stages (2) and (3), both reflection measurements and ZOPs were used.

For reflection measurements, one single antenna acts as transmitter and receiver. It is lowered down in a borehole, and traces are collected every 0.25 m over an interval from a 0.8 to 17.3 m depth of the midpoint between transmitter and receiver. Measurements are done in the wells D02, D04, D05, D06, D08, C04, C05, U03, U04, U05, and U06. The remaining wells outside the freezing area are occupied with other instrumentation, and the wells between the heat exchangers are inaccessible during freezing. A standard processing is applied, consisting of subtraction of DC shift, zero-time correction, bandpass frequency filtering using a third-order Butterworth filter with cut-off frequencies 50 and 600 Hz, and a gain correcting for spherical divergence. The maximum phase of the reflection arrivals is picked. Maxima are picked because, due to the low signal-to-noise ratio, first-break picking is too erroneous. To account for the time shift between the first break and the first maximum, the travel time of the first maximum of the direct wave is used for zero-time correction. Because the emitted signal is assumed to be reflected at the freeze–thaw boundary and travel back to the receiver (Fig. 1a), with known antenna position and propagation velocity, the distance to the appearing reflectors can be calculated. For distance calculation, borehole-specific velocity–depth profiles (see Table 1) are taken from the crosshole measurements described in the following paragraph.

Table 1 Velocity profiles from ZOPs used for time–distance conversion of reflections.

For crosshole zero-offset profiles, two antennas are placed in two boreholes, with one used as a transmitter and one used as a receiver. The antennas are lowered, simultaneously measuring at the same depth every 0.25 m over a depth interval from 0.8 to 17.3 m. The resulting ZOPs are processed with the same processing flow as the reflection profiles. Zero time is determined by picking the maximum of first arrivals in air measurements outside the wells. Source–receiver distance is increased from 1 to 3 m in 0.2 m steps in the zero-time measurements. Subtracting the calculated travel time for v=c in air leaves the value for zero-time correction. From the first arrival maxima in the ZOPs, we can calculate the bulk velocity of the soil between the antennas. When freezing happens, we assume higher velocity and lower attenuation that will result in lower travel times and higher amplitudes. To obtain a 3D distribution of wave velocity in the experimental field, different well combinations are realized (Fig. 2a). Well combinations covering all three sides of the test field ensure velocity information in close vicinity to all reflection measurements.

Figure 2 Crosshole ZOP for (a) site characterization and (b) freezing monitoring.

For the monitoring of freezing and thawing, only measurements running through the freezing area are of interest (Fig. 3b). Measurements between wells outside the freezing area serve as a reference to check for acquisition errors. If the measured velocity remains the same for the reference measurements, where there is no change in subsurface properties, the data acquisition is considered to be correct.

Figure 3 (a) Error of velocity estimation in crosshole measurements. Error depends on source–receiver distance and travel time. The black line indicates expected values for velocity of 0.061 m ns⁻¹. White circles highlight ranges for site characterization (A) and freezing monitoring (B). (b) Distance error for reflector position estimation via reflection measurements for different velocity errors.

2.2.3 Preliminary considerations – accuracy

The accuracy of velocity estimation in ZOPs is mainly influenced by positioning error and time error. Antenna positioning error consists of the GPS error and the inclinometer error. Time error is mainly composed of the picking error and the zero-time drift of the aperture. The time error was assessed by comparing travel times of repeated measurements on multiple days of the same well combination. Figure 3a shows the expected velocity error assuming a positioning error of Δs=0.05 m and a time error of Δt=0.5 ns. Greater distances between source and receiver minimize the velocity error. The downside of a greater source–receiver distance is that lateral variations are averaged. Smaller distances ensure that the measured velocities, which are then used for time–depth conversion of reflection measurements, represent the immediate vicinity of the observation wells. So, for site characterization, a trade-off between high velocity error and spatial averaging is selected with well distances of 2 to 4 m, resulting in an expected maximum velocity error of $\mathrm{\Delta }v<\pm \mathrm{0.0025}$ m ns⁻¹ (Fig. 3a). ZOPs for freezing monitoring, where the wave travels through the whole experimental field covering distances of around 6 m, reduce velocity error to $\mathrm{\Delta }v\approx \pm \mathrm{0.001}$ m ns⁻¹.

The low signal-to-noise ratio in reflection measurements results in an increased time error of $\mathrm{\Delta }t=\pm \mathrm{0.75}$ ns. Using the velocities from site characterization for the determination of reflector position, velocity error, and time error results in a distance error of $\mathrm{\Delta }d=\pm \mathrm{0.1}$ m (Fig. 3b) at a reflector distance of 1.5 m. The 1.5 m corresponds to the distance of the surrounding observation wells to the heat exchanger probes. As the freeze–thaw boundary moves outwards, and therefore the distance to the observation well decreases, the error decreases.

Vertical positioning was done by aligning depth markings along the antenna cable with the top of the well casing. We assume the error of this method to be smaller than 0.01 m; therefore it was neglected.

The estimation of the lateral extent of the frozen body from ZOP velocities is assessed. The distance s traveled by the wave consists of the distance traveled in the thawed medium and in the frozen medium.

Travel time t consists of the time traveled in the thawed medium and in the frozen medium:

$$\begin{array}{}\text{(1)}& s=s_{\text{ice}}+s_{\text{thawed}},t=t_{\text{ice}}+t_{\text{thawed}}.\end{array}$$

Inserting

$$\begin{array}{}\text{(2)}& t={\displaystyle \frac{s}{v}}={\displaystyle \frac{s_{\text{ice}}}{v_{\text{ice}}}}+{\displaystyle \frac{s-s_{\text{ice}}}{v_{\text{thawed}}}}\end{array}$$

and solving for s_ice yields

$$\begin{array}{}\text{(3)}& s_{\text{ice}}={\displaystyle \frac{s\cdot v_{\text{ice}}(v-v_{\text{thawed}})}{v(v_{\text{ice}}-v_{\text{thawed}})}}.\end{array}$$

If the permittivity, and therefore velocity, of the frozen saturated soil is unknown, the equation is underdetermined. Distance s is known from GPS positioning corrected for borehole deviation, v_thawed is taken from baseline measurements, and v is measured velocity. There is no access to direct information on v_ice in the test field because wells within the freezing zone are inaccessible during plant operation. Nonetheless, we expect to get qualitative information about the beginning of freezing and the shape of the ice body by looking at the changes in bulk velocity.

Figure 4 (a) Example profile from well C09 to C12. (b) Velocities calculated from the first arrivals of all ZOPs shown in Fig. 2a. (c) Amplitude of first arrivals. (d) Results of HPT and EC logging. (e) Sediment core at drilling MP055 (for the position, see Fig. 1b).

3 Results

3.1 Site characterization

Logging with the direct-push-based hydraulic profiling tool (HPT) at drilling MP055 (Fig. 4d) shows a layer with high relative hydraulic conductivity, K_HPT≈10 $\mathrm{mL}{\min }^{-\mathrm{1}}{\mathrm{kPa}}^{-\mathrm{1}}$, and an electrical conductivity of σ≈10 mS m⁻¹ between 10 and 17 m depth. Above a layer with lower hydraulic conductivity, higher electrical conductivity is observed. The low conductivity spans up to 7 m depth, while electrical conductivity is higher at 10–9.5 m and at 7 m depth. Below 17 m, electrical conductivity increases to σ=15–20 mS m⁻¹ and K_HPT decreases to ∼0.5 $\mathrm{mL}{\min }^{-\mathrm{1}}{\mathrm{kPa}}^{-\mathrm{1}}$. Sediment coring at MP055 (Fig. 4e) matches the high hydraulic and low electrical conductivity with a sand layer and matches the low hydraulic and high electrical conductivity with a higher clay content. We assume a subsurface model with a sandy aquifer at 10 to 17 m depth covered at the top and bottom with an aquitard. The low-conductivity zone is very heterogeneous at this location (see supplement of Löffler at al., 2022). This can also be seen in Fig. 4b and c: over the small test field, we see most variations in amplitude and velocity precisely between 7 and 10 m, indicating the presence of a highly heterogeneous boulder clay. We correlate the low amplitude and low travel time with a highly heterogeneous layer, acting as a low hydraulic conductivity layer due to the presence of clayey material.

The ZOPs measured before the installation of the heat exchangers offer the possibility for an enhanced characterization of the experimental site itself. ZOPs were used to spatially extrapolate the information over the area of the test site. One example profile from wells C07 to C10 is shown in Fig. 4a. In the first 3 m the first arrivals are superimposed by a signal refracted at the surface. Up to a depth of 7 m, high amplitudes and longer travel times, converting to velocities of v=0.06–0.065 m ns⁻¹, can be seen. From ∼7 to ∼10 m and from ∼16.5 m on, there are layers with lower amplitudes and shorter travel times converting to velocities of v=0.07–0.075 m ns⁻¹. Between 10 and 16.5 m, very high amplitudes and longer travel times, converting to velocities of 0.06–0.065 m ns⁻¹, are measured. The layering corresponds to the identified layers in the core drilling and in HPT and EC logging. The aquifer appears as a layer with high amplitudes due to lower electrical conductivity and with longer travel times due to higher permittivity, and the aquitard appears as a layer with low amplitudes due to high electrical conductivity and with shorter travel times due to lower permittivity. Plotting all measured well combinations shows that the layering of an aquifer covered on the top and bottom with an aquitard is seen in the velocity and amplitude in all profiles but with the thickness and velocity of the upper aquitard being variable. Placing the velocities at the midpoint between transmitter and receiver gives with the 2D view an idea of the 3D distribution of the velocities over the experimental site (Fig. 5). Figure 5a is with a fixed north coordinate, and Fig. 5b is with a fixed east coordinate. A slight velocity increase is seen from southwest to northeast. The thickness of the upper aquitard also increases from southwest to northeast.

Figure 5 Wave velocities placed at the midpoint between source and receiver. Lateral 2D projection of the 3D distribution in an (a) east and (b) north direction.

3.2 Reflection profiles

Reflection measurements were performed at the final equipped site before the start of the freezing experiment and after 54 d of running the LHS system, with an ice body around the heat probes. In the baseline measurements, linear reflectors are visible in reflection profiles. An example is shown in the profile from well D02 (Fig. 6a). In the lower aquifer, reflectors are identifiable to a maximum travel time of 80 ns. In the aquitard, where attenuation is higher, only close reflectors are visible. From 0 to 4 m depth, the well reflections are superimposed by reflections at the surface and in the water table. Picking travel times of the reflectors and converting to distance, using the velocity estimation from ZOPs, this corresponds to the distance to other observation wells. So, in close vicinity, reflections of the surrounding observation wells occur in the data. Not all surrounding wells can be identified because of the low signal-to-noise ratio and superimposing signals. By identifying these reflectors as well reflections before the freezing is initiated, the baseline measurements enable us to distinguish between well reflections and the signal reflected at the freeze–thaw boundary.

Figure 6 (a) Reflection profile from baseline measurement at well D02 with visible reflectors. (b) Picked reflection travel times converted to distance using velocity from ZOPs. Reflector positions coincide with the distance to other observation wells.

Figure 7 displays profiles from well U04 before freezing and with a developed ice body during plant operation. A strong reflector appears at a travel time of ∼50 ns at a depth of 10 to 16 m. The hyperbolic increase in travel time at the upper and lower end of the reflection indicates that the reflecting object ends there, rather than that the signal is just not visible because of the high attenuation of the surrounding layers.

Figure 7 Section of a reflection profile from (a) baseline measurement and (b) during freezing. Appearance of a new reflector at 10–16 m depth.

In a second step, travel times of new occurring reflections are picked in all reflection profiles. After converting to distance, the extent of the ice body is estimated. Because the antennas transmit the signal omnidirectionally, the possible reflector origin is a sphere around the antenna position. Assuming the reflector origin is at the same depth as the antenna, ±0.2 m, the possible reflector positions of each measurement are plotted on a discretized 3D space. The positions closest to the test site are interpreted as the freeze–thaw boundary, and connecting the edges gives an estimation of the current lateral ice volume. Figure 8 shows the estimated ice boundary at a depth of 15 m with an extension of 4.2 ± 0.2 m in a southwest to northeast direction.

Figure 8 Estimated ice boundary based on possible reflector origin.

3.3 Crosshole ZOPs

For the monitoring of freezing and thawing, only the well combinations crossing the whole freezing area, D04-U04, D04-U06, and D05-U05, are of interest (Fig. 2b). Perpendicular measurements are not possible due to the lack of observation wells in the southeast. The profiles D04-D05 and U06-U05 outside the freezing area serve as a reference. They show unchanged velocities within the error range.

Figure 9 shows two profiles through the freezing area from well D04 to U04 and D05 to U05 before (Fig. 9a and c) and during (Fig. 9b and d) freezing. Due to the greater distance between wells and the high attenuation, no signal is registered in the aquitard. With the formation of ice, the first arrival times at 10 to 16 m depth drop significantly and reverberations occur. As an example, velocity at 15 m depth increases from 0.0605 ± 0.001 m ns⁻¹ to 0.099 ± 0.001 m ns⁻¹. The signal at a depth greater than 16 m arriving with an unchanged travel time of ∼100 ns indicates no freezing there. The changes shown also appear in the third profile crossing the freezing area, which is not shown here.

Figure 9 ZOP between D04-U04 and D05-U05 before (a, c) and after (b, d) 54 d of subsurface freezing.

4 Conclusion

We conducted borehole GPR measurements in a shallow Quaternary glacial aquifer before and during operation of an experimental geological latent heat storage. Prior to this study, it was difficult to predict the geological conditions in the immediate vicinity of the experimental site. We show that it is possible to extrapolate punctual geological information from one drilling point using a 3D velocity model derived from simple zero-offset crosshole measurements. This allows us to include detailed geometries of the geological layering in thermo-hydraulic modeling approaches.

Furthermore, while operating the LHS system, we aimed to investigate the feasibility of imaging a vertical freeze–thaw boundary using borehole GPR in the challenging (for GPR) context of water-saturated glacial sediments. We therefore carried out reflection and crosshole borehole measurements during an LHS storage experiment, in which a depth-horizontal subsurface volume is frozen. These experiments confirmed that both borehole crosshole and borehole reflection GPR enable us to image the frozen subsurface volume. However, only reflection measurements are able to quantify ice body size by determining the position of the freeze–thaw boundary with an error of ±0.1 m. Measuring at multiple campaigns has shown fast acquisition and good repeatability of the data. Resolution is mainly limited by the timing error of the wave arrivals caused by low signal-to-noise ratio because of high attenuation in a water-saturated environment. This prevents the use of higher-frequency sources for reflection imaging. At least two boreholes are required to obtain accurate velocity information and to map the extent of freezing in one direction. There is evidence for the vertical confinement of the ice body, but the clear determination of top and bottom is limited by Fresnel zone width.

For further projects, observation well positions that are in same distance range as the expected freeze–thaw boundary have to be avoided. In this study, the monitoring design was affected by having too many observation wells at a similar distance to the imaging target due to the requirements of providing access not only for geophysical monitoring.

Taken together, these results suggest that borehole GPR is a viable method for monitoring LHS systems. The combination of ZOPs and reflection measurements is a suitable setup for quick imaging of the lateral boundary of a freezing subsurface volume.