Characterizing soil piping networks in loess‐derived soils using ground‐penetrating radar

Soil piping remains a relatively unexplored phenomenon despite its substantial impacts on water and sediment transfer at the watershed scale at numerous locations around the world. One of the main limitations regarding the study of this singular process is the characterization of the pipe networks (in terms of number, position, dimensions, and connectivity of pipes). In this context, noninvasive subsurface imaging using ground‐penetrating radar (GPR) has proven to be a promising technique. This study used two three‐dimensional (3D), high‐resolution GPR surveys performed in loess‐derived soils to characterize small pipe networks with little prior information about their location. The adopted methodology relies on high spatial resolution scanning, 3D subsurface imaging, and automated detection of reflection hyperbolas using a 200‐MHz center‐frequency antenna. Two small watersheds affected by piping were investigated at Sippenaeken and Kluisbergen (Belgium). Over the two scanned zones, results revealed significant subsurface continuous patterns. Even though the most obvious patterns corresponded to recent or past anthropic activities (e.g., artificial drainage pipes), validation tests confirmed that the chosen methodology can be used for pipe network characterization because some important continuity patterns were related to soil piping: for instance, a fairly small pipe (approximately 10–15 cm in diameter) was detected and validated for more than 100 m. Nevertheless, the high variability in size, depth, and orientation of the pipes imply that GPR may only be truly efficient when using very high spatial resolution scanning, which limits its application to specific conditions not always met in piped areas.


INTRODUCTION
In general terms, soil piping refers to the formation of subsurface pipes due to the erosive action of water flowing (Jones, 2010). With very few exceptions, wherever the piping process has been observed there is evidence that it contributes significantly to watershed hydrology (Carey & Woo, 2000;Chappell, 2010;Jones, 1994Jones, , 1997Jones, , 2010Smart et al., 2013;Uchida, Kosugi, & Mizuyama, 2001;Vannoppen, Verachtert, & Poesen, 2017;Wilson et al., 2017). However, despite these studies, the origin and hydrological functioning of soil pipes is not yet fully understood, in particular in the specific context of loess-derived soils (Vannoppen et al., 2017). This deficiency is partly due to the complexity of characterizing pipe networks (Bryan & Jones, 1997). Indeed, their underground and heterogeneous nature makes it very difficult to estimate the number of pipes, their dimensions, positions, and connectivity. Yet, without exact knowledge of the pipe network configuration, it remains very complex to assess precisely their role and their link with the different components of the hydrological cycle.
Until today, most research has relied on geomorphological mapping (observations of piping-related features such as pipe collapses, inlets, and outlets) to investigate pipe networks (Bernatek-Jakiel & Kondracka, 2016;Faulkner, Alexander, & Zukowskyj, 2008;Romero Díaz, Marín Sanleandro, Sánchez Soriano, Belmonte Serrato, & Faulkner, 2007;Verachtert, Van Den Eeckhaut, Poesen, & Deckers, 2010;Wilson, Rigby, & Dabney, 2015;Zhang & Wilson, 2013;Zhu, 2012). Unfortunately, these point-measurement techniques do not allow accurate subsurface characterization, and they probably underestimate pipe lengths and network densities (Bernatek-Jakiel & Kondracka, 2016;Holden, Burt, & Vilas, 2002). In their recent review, Bernatek-Jakiel and Poesen (2018) noted the strong need to develop new techniques to survey pipe networks. They emphasized, among others, the potential of noninvasive geophysical techniques such as ground-penetrating radar (GPR). Subsurface imaging using GPR indeed seems promising given its ability to detect heterogeneities in soil profiles when there is a sharp variation in electromagnetic properties (as should be the case with pipes). Holden (2004Holden ( , 2005 and Bernatek-Jakiel and Kondracka (2016),  already demonstrated the potential of GPR for soil pipe detection.  used GPR for the identification of subsurface piping in blanket peat in the northern Pennine Hills (UK). In their study, low-frequency GPR systems (100-and 200-MHz antennas) were used above 30 known pipe locations and, when pipes were detectable, transects were performed to assess their upslope extension. Their work confirmed that GPR can be used to identify these natural pipes (the smallest verified pipe identified by their GPR system was 9 cm in diameter) and to estimate their depths (with an accuracy of 20-30 cm). Results also suggested that pipe densities were much greater than could be assessed from the surface. This was later confirmed by Holden (2005) by applying a similar methodology in 160 blanket peat catchments, finding soil pipes in all Core Ideas • 3D GPR surveys were conducted in loess-derived soils to characterize soil pipe networks. • Small soil pipes extending for several tens of meters were successfully detected. • High spatial resolution scanning is essential. • Automated detection of reflection hyperbolas greatly facilitates pipe network mapping.
catchments with a mean frequency of 69 pipes per kilometer of surveyed transect. However, their results also showed that GPR on its own did not allow the pipe network connectivity to be assessed . To do so, Holden (2004) had to place GPR above known pipe locations and inject a tracer solution detectable by the radar (NaCl) in upslope pipe cavities. Although very effective, the methodology requires the pipes to contain pipe flow (Holden, 2004) but also requires prior information about pipe locations and inlets for tracer injection.
Bernatek-Jakiel and Kondracka (2016) used GPR for the assessment of piping activity in Cambisols in the Bieszczady Mountains (Poland). In their study, a higher frequency GPR system (500 MHz) was used above tens of potential pipe locations in combination with electrical resistivity tomography and geomorphological mapping. Their study also demonstrated that the subsurface density of pipes was significantly greater than that mapped on the basis of surface observations alone (for instance, the underestimation of pipe length reached a maximum of 49% for piping systems in an early stage of development) and that GPR may lead to a better estimation of the scale of piping. However, their work also underlined that surface indicators of piping activity are necessary to confirm often ambiguous geophysical data.
Soil piping has been reported in various parts of Belgium in loess-derived soils, predominantly on grassland and cropland. In the Flemish Ardennes, Verachtert et al. (2010) reported pipes (observed in collapses) with diameters ranging between 0.1 and 0.4 m. These pipes are mostly located between 0.7-and 1.4-m depth. Vanderputten (2006) first reported soil piping in the Limburg province. In this location, repeated field observations have shown the presence of pipes whose diameters can vary abruptly across short distances from a few to tens of centimeters or more at depths up to 1.4 m. The direction of these pipes can change rapidly both in the horizontal and vertical planes, sometimes by up to 90 • . Moreover, because piping constitutes a hazard for cattle as well as for agricultural traffic, farmers quickly fill up collapsed pipes with heterogeneous material including, e.g., bricks. All this makes for highly heterogeneous-and therefore challenging-conditions. Combined with the high silt contents resulting in important attenuations of the GPR signal when the soil is wet, these conditions strongly limit pipe detection based on single GPR profiles, as experienced by Verachtert (2011) in an explorative study in the Flemish Ardennes.
The aim of this research was therefore to evaluate the ability of high-resolution GPR scanning, combined with automated hyperbola detection, to characterize networks of highly variable pipes with little prior information about their location and heterogeneous soil conditions.

Sippenaeken
The first experimental site chosen for this study is located in Sippenaeken (Gueule Valley,eastern Belgium,50.746 • N,5.934 • E). Piping at this location was first described by Vanderputten (2006). This hilly region is characterized by a maritime temperate humid climate with mild winters and a mean annual precipitation of 930 mm, relatively well distributed throughout the year. A maize (Zea mays L.) crop was grown in the downstream part of the site, and pasture was found in the upstream part ( Figure 1). Its mean elevation is approximately 170 m asl and it presents an 8-18% northeast-facing slope ( Figure 1). Within this area, the soil mainly consists of Luvisols derived from Quaternary loess deposits (textural class: silt loam), but a hydromorphic soil strip also crosses the field (Service Public de Wallonie, 2013). According to the Belgian geological map, the site presents two underlying materials: Aachen Formation in the upper part and Houiller in the lower part of the watershed (Service Public de Wallonie, 2014a).
Before the measurement campaign, only two pipe collapses and one flowing pipe outlet had been observed in the lower and steeper part of the field grown with maize ( Figure 1). Even though these collapses were large (30 cm wide, 60 cm high, and 1 m long for the upper one), they were related to quite small pipes (diameter range 8-15 cm). A tracer test performed between the collapses and the outlet demonstrated that they were connected. Additional collapses are known to have occurred on this field in the past, but they are no longer visible because the farmer filled them and plowed his field yearly. Unfortunately, high spatial resolution scanning was impossible within the maize field because of the presence of the maize in spring and summer and excessive roughness in autumn and winter (plowed field). Therefore, a 170-by 110-m zone was delimited upslope in grassland for scanning (Figure 1), where no prior information on pipe location was available. In the northeastern part of the scanned zone, an artificial polyvinyl chloride (PVC) drainage pipe was recently installed.

Kluisbergen
The second study site is located in Kluisbergen (Scheldt Valley, western Belgium, 50.769 • N, 3.554 • E). Piping at this location was first described by Verachtert et al. (2010). This hilly region is also characterized by a maritime temperate humid climate with mild winters and a mean annual precipitation of approximately 800 mm, relatively well distributed throughout the year. The experimental site is covered by pasture. Its mean elevation is approximately 75 m asl and it presents a 4-29% north-facing slope. The soils mainly consists of Cambisols (textural class: silt loam) derived from Quaternary loess deposits, but poor draining Retisols are also found in the upper part of the field (Agentschap Geografische Informatie Vlaanderen [AGIV], 2001b). According to the Belgian geological map, the Tertiary lithology consists of homogenous blue massive clays containing >50% clay (Aalbeke Member) (AGIV, 2001a).
As shown on Figure 2, this second experimental site is characterized by the presence of one outlet and numerous pipe collapses. Some of these collapses were unfilled sinkholes, but the vast majority were filled. In this case, the filled sinkholes were discernible from the surface because the poor quality backfill material was preferentially colonized by specific types of vegetation contrasting with the rest of the field (e.g., stinging nettle [Urtica dioica L.]). Given the apparently random distribution of the collapses across the entire field, it was not possible to estimate the potential density of the network F I G U R E 1 Map and topography of the experimental field, location of known pipe collapses and pipe outlet, and positions of the 170-by 110-m scanned zone at Sippenaeken, Belgium. (Service Public de Wallonie, 2014b) or even the pipe orientations without additional information. Because two active pipe collapses (unfilled sinkholes, with water flowing inside) were observed near the thalweg and at a reasonable distance from the outlet, a 100-by 50-m area was delimited downslope for GPR scanning ( Figure 2). As at Sippenaeken, even though the collapses presented large dimensions (one of them was 2.5 m wide, 1 m high, and 4.5 m long), field observations demonstrated that they were related to relatively small pipes (approximately 15 cm in diameter). Although artificial drainage systems have been shown to enhance piping in this region (Verachtert et al., 2013), there was no record of such features at this study site.

Field equipment and measurements
At Sippenaeken, a time-domain GPR system (Model SIR-20, Geophysical Survey Systems Inc.) was used as an on-ground GPR equipped with transmitting and receiving 200-MHz center-frequency bowtie antennas. A 200-MHz center-frequency antenna was selected because it offers a good trade-off between horizontal resolution (able to detect potential pipes of about 5-10 cm in diameter or larger), depth penetration (around 4 m under field conditions but with an unfavorable signal/noise ratio after 2 m), and time consumption (higher spatial resolution is required when using higher center-frequency antennas). At Kluisbergen, a multifrequency time-domain GPR system (Model RIS MF Hi-Mod, Ingegneria dei Sistemi) was used as an on-ground GPR equipped with two transmitting and two receiving 200-and 600-MHz centerfrequency bowtie antennas, but technical problems during the field acquisitions prevented using the 600-MHz frequency.
At Sippenaeken, the 170-by 110-m zone was scanned during the summertime with a high spatial resolution perpendicularly to the slope in the x-axis direction as shown in Figure 1. Each transect was 140 to 170 m long and the separation between two consecutive transects was 0.5 m. Measurements were performed every 25 mm along the 220 transects. Three days were necessary to complete the survey.
At Kluisbergen, the 50-by 100-m zone was scanned during spring with a high spatial resolution perpendicularly to the slope in the x-axis direction as shown in Figure 2. Each transect was 50 m long and the separation between two consecutive transects was 0.25 m. Measurements were performed every 10 mm along the 400 transects. Four days were necessary to complete the measurements. A differential global positioning system (dGPS) was used to geographically position the scanned zones in the fields. It was also used to check the position of the radar profiles. The positioning of the radar profiles within the zones was done using measuring tapes laid on the ground surface. The position of the radar measurements along the profiles was recorded with the odometer of the survey wheel of the radar system. We used a local coordinate system for building the radar images, where z is the depth from the soil surface and x and y are the positions along the profiles and of the profiles following the soil surface topography, respectively. To set the scanning direction and ensure that the survey lines did not cross, a mason rope was stretched between the ends of the profiles.

Preprocessing
To enhance the characterization of pipe networks, standard data processing procedures were first applied to the individual radar images (different preprocessing configurations were tested on the whole dataset). First, a spatial frequency filter (f-k filter) was applied to remove relatively constant reflections from the signal (typically removing reflections with a spatial frequency lower than 0.005-0.01 m −1 ). These reflections mainly correspond to the soil surface and the multiple internal antenna reflections for the incident field, the latter being independent of the medium. To limit the effects of noise and clutter, a band-pass frequency filter was applied to the dataset, with varying bandwidths (typically around 25-500 MHz). Finally, a power gain function was applied to compensate for spherical divergence and attenuation losses using where S is the initial radar signal, t is the travel time, S g is the amplified radar signal, and a is the applied gain (typically 2.5-3 in this study). This preprocessing allows highlighting of the local heterogeneities, but it is also an important step for the subsequent automated hyperbola detection algorithm. Indeed, it decreases significantly the number of pixels to be analyzed after the application of the Canny edge detector (for instance by removing important constant reflections) and therefore significantly reduces the computational time (Mertens, Persico, Matera, & Lambot, 2016

Three-dimensional subsurface imaging
The numerous parallel transects recorded over the scanned zones were combined to construct two 3D images of the subsurface. To visualize and interpret these 3D images, they were represented as successive two-dimensional (2D) depth slices. Depth slices represent the amplitude of the signal at given depths as a function of the horizontal coordinates. Thus, they permit detection of horizontal structures in the 3D radar images, e.g., such as a pipe. For this study, because pipes are not horizontal and their depth may change in relatively short distances, we calculated depth slices by averaging the signal amplitude across a specific depth range, i.e., 0.01-0.50 m. It is noteworthy that depth slices must be analyzed as a whole and not only individually to assess the spatial continuity of a soil pipe. Indeed, field observations in the lower part of the field in Sippenaeken have revealed that the depth and diameter of pipes may vary significantly, even across short distances, and the same pipe may thus appear on several consecutive depth slices. Depth was determined by considering a constant relative permittivity of 8, estimated through field tests (outside of the scanned zones, over known pipe locations). As the relative permittivity is expected to vary with depth, space, and the presence of local heterogeneities, this approximation may lead to errors in depth estimation. It should not, however, affect 3D image interpretation and pipe network characterization.

Automated detection of reflection hyperbolas
As mentioned above, soil pipes act as discrete objects if crossed transversally, and hence they appear on radargrams as reflection hyperbolas. Detecting and mapping their occurrence in consecutive transects should therefore allow characterization of the spatial continuity of soil pipes (or other heterogeneities).
Recently, Mertens et al. (2016) developed a novel methodology for automated detection of reflection hyperbolas in complex media for which hyperbolas are not well defined (distorted, incomplete, etc.). This methodology is based on a detection algorithm calibrated based on the way humans recognize patterns. The methodology first relies on a Canny edge detector (Canny, 1986), which binarizes the image and highlights the reflections of interest. Once completed, the algorithm determines the apex of the hyperbolas by trying to fit a set of potential hyperbolas (physically plausible considering various combinations of relative permittivity and radius) to the profile edge dots defined by the Canny filter. The existence of a hyperbola is further determined using a set of carefully chosen criteria calibrated so as to fit the ambiguity zone for the pattern recognition process. The inherent distortion of field hyperbolas is further considered by defining a buffer zone around the theoretical hyperbola. This methodology, summarized quite simply here, is described in detail by Mertens et al. (2016).
The main advantage of this algorithm is its good performance in complex media where reflection hyperbolas are usually ill shaped (and therefore complex to detect automatically). The most sensitive step of the whole procedure is the application of the Canny edge detector because it determines whether a pixel will be further analyzed as part, or not, of a reflection hyperbola. It depends on two threshold values-if they are too restrictive, few hyperbolas are detected (only the strongest ones) but computing times and false alarm rate are lower and vice versa. Several threshold values were tested and used in this study.
To assess the spatial continuity of reflection hyperbolas, results of the automated detection algorithm were visualized on a single plan view. In practice, positions (in the x-y axes direction) of every detected hyperbola were recorded and combined on a single plan view of the scanned zone. Their detection depth was estimated considering a relative permittivity of 8 and then represented using a color index.
It is noteworthy that hyperbolas are identified in the 2D GPR profiles (there is no identification of hyperboloids in the 3D data). Therefore, heterogeneities located near the profile lines will generate hyperbolas with a shape that does not exactly correspond to their properties and position. However, this should not affect results interpretation and pipe network characterization, which focus on detecting significant continuous patterns in several consecutive transects.

Soil pipe network mapping
Results of each technique (3D subsurface imaging and automated detection of reflection hyperbolas) was first compared and then used to map subterranean structures and potential pipe networks. For each continuous pattern observed, positions and approximate depths were recorded and further used to construct the map. When continuity between two structures was difficult to assess with the two selected techniques, a visual analysis of the consecutive 2D transects of interest was also performed.

Validation
Once the mapping was completed, validation trenches were dug along the main observed continuity patterns, the key axes of a potential pipe network, to verify the presence of soil pipes.
In total, six to seven trenches (2-4 m long, 0.75 m wide, and 1-3 m deep) were dug at each study site.

Soil pipe network mapping
Three-dimensional subsurface imaging and automated detection of reflection hyperbolas both reveal interesting spatial patterns throughout the entire scanned zone. Regarding 3D subsurface imaging, an example of one depth slice is presented in Figure 3 (approximate depth of 0.8-1 m). Results of the automated detection are presented in Figure 4 for quite restrictive threshold values. As it appears on Figure 4, most of the detected hyperbolas are located between 0.5 and 2 m (mean detection depth of 1.12 m with a standard deviation of 0.47 m). This is in part due to the low signal/noise ratio below the 2-m depth (50 ns) such that only the strongest hyperbolas are detectable below this depth but also because automated detection was limited to the part of the signal located below the surface reflection to avoid an excessive false alarm rate in the 0-15-ns time range.
Conceptually, the scanned zone can be divided into three distinct geographical zones within which the radar signal presents similar features: the left zone (0-50 m along the x axis), the middle zone (50-110 m), and the right zone (110 m to the end of the transects). The most significant continuous patterns were observed with both methods in the middle part of the scanned zone (Figures 3 and 4). Here, several axes appeared continuous across several tens of meters. Moreover, even though it is difficult to assess with certainty, some of these axes seem somehow connected to each other and could therefore correspond to a potential pipe network.
In the left and right parts of the scanned zone (Figures 3  and 4), results are more difficult to interpret. Indeed, these two zones presented a large number of seemingly randomly distributed subsurface reflectors, often much stronger than the continuous ones observed in the middle part. Only one continuous pattern could be clearly observed with both methods. It is located to the bottom right of the scanned zone and corresponds to a known artificial drainage pipe location (x = 120 m, y = 15 m in Figures 3 and 4). Other potential axes are visible, generally more clearly with the automated detection method, but they are limited in length (maximum 5-7 m long) and appear isolated from each other (some examples are framed in Figures 3 and 4).
In addition to the previously described features, the 3D subsurface image also revealed the presence of large layers whose depth increased through consecutive depth slices. Such layers show up as strong reflections whose positions shift gradually with depth. An example of such a reflection can be found in Figure 3 [almost horizontal, from coordinates (60, 15) m to (150, 25) m], but other similar features were visible in other depth slices for different parts of the 3D image. These layers correspond to strong variations in the soil dielectric properties that may relate to the presence of clear distinct soil layers (for instance, a clay layer under the loam).
Results from both methods were combined to map the potential soil pipe network. As explained above, neither method revealed significant continuous patterns in the left and right parts of the scanned zone. Therefore, network mapping was limited to the middle part of the scanned area, where both methods revealed identical continuity patterns. The detected axes present similar levels of discernibility with both methods. This gives confidence when they are easily visible but also implies that there is no strong added value to using both methods when axes are less perceptible. To lower the uncertainty associated with the less obvious axes, a visual analysis of consecutive 2D radargrams was performed. However, the presence of either a large number of small reflectors or small breaks in the spatial continuity (absence of reflection for one or two transects) often lead to ambiguous interpretations (continuity or not). Nevertheless, a potential pipe network could be drawn ( Figure 5).

Validation
Several validation trenches were dug along the main continuous patterns suggesting the presence of a pipe ( Figure 5), except for one of them. Indeed, a few days before the validation, a nascent pipe collapse occurred within the scanned zone ( Figure 6a). Its position was recorded using a dGPS and compared with the network map. As it appears in Figure 5, the position of the collapse corresponds with one continuous pattern. To allow future monitoring, it was decided to safeguard the pipe integrity and no validation trench was dug. The collapse was thus considered as a confirmation of pipe presence. For the other continuous patterns, the validation results were the following (numbers have been attributed to the different trenches, as shown in Figure 5): • Trenches no. 1, 2, and 3: Three trenches were dug along this curved axis. It was thoroughly investigated as it behaved singularly on the radar data. Indeed, it presented patches of strong reflections clearly visible with both methods. These patches seem to be linked to each other by weaker reflections, but connectivity could be ascertained only by visualizing consecutive 2D radargrams. Two trenches were dug along these strong reflection patches and, after 40-50 cm, backfill material was found in both trenches. The backfill was composed of dark-colored soils mixed F I G U R E 3 Depth slice (approximate depth of 0.8-1 m) illustrating reflections from various types of features at Sippenaeken. Color scale represents the reflection intensity with pieces of anthropic materials such as bricks that made further digging impossible. These materials correspond to strong reflection hyperbolas, almost occurring in the surface reflection. This backfill material is very likely linked to a recently collapsed active pipe. Indeed, its morphology matches with the farmer's current way of dealing with pipe collapses: filling the sinkhole with poor quality materials (e.g., bricks) and covering the last 30-50 cm with good quality soil. As observed in other parts of the field, this filling method generally does not stop pipe flow because water usually finds its way through this coarse backfill material. The third trench partly confirmed this theory. Indeed, it was dug between two strong reflection patches and a flowing soil pipe was found at the 90-cm depth (7-8 cm in diameter, Figure 6b). • Trench no. 4: At the 80-cm depth, a singular 15-cm-wide by 15-cm-high stone line was found crossing the trench (Figure 6c). It was made of 5-10-cm size stones, the space between stones being largely filled with soil. It was oriented exactly as can be seen on the network map. Considering a relative permittivity of 8, its depth corresponds to the expectations (reflection hyperbolas located approximately 15 ns below the surface). Its origin remains uncertain but it seems linked with past anthropic activities. There was no sign of recent water flow.
• Trench no. 5: As for trench no. 4, after 60 cm, another singular 15-by 15-cm stone line was found crossing the trench. Again, it was oriented exactly as on the network map and its depth corresponds to the expectations (reflection hyperbolas located approximately 10-15 ns below the surface). • Trench no. 6: In this part of the scanned zone, no pipes or other obvious structures were found (the trench was 2 m in depth). This result could be partially expected because the trench was dug in the less discernible part of this axis. • Trench no. 7: As for trench no. 6, no pipes or other obvious structures were found (the trench was 2 m in depth). This result was more surprising because both methods revealed a clear continuity pattern.
Overall, five validation points out of seven were associated with identifiable reflectors. Two significant continuity patterns could be related to the presence of soil pipes. Are these two axes connected as suggested by the network map? Field observations made a few months after the validation campaign partly confirmed this idea. Indeed, during a winter storm event, two small ephemeral outlets were observed in this part of the field (Figure 6d). These ephemeral outlets consisted of small field mouse burrows connected to the underground pipe and discharging only during large storm Vadose Zone Journal F I G U R E 4 Automated hyperbola detection results for the entire scanned zone events (an example of such a connection observed in the lower part of the field is presented in Figure 6e). Their positions were recorded with a dGPS to compare them with the network map. As it appears in Figure 5, one ephemeral outlet is located at the junction between the two axes and may indicate that they belong to the same pipe system. The other one is located just above the lower axis, confirming our validation results. Given these validation results, a schematic representation of the confirmed subsurface structures is presented in Figure 7.

Soil pipe network mapping
Three-dimensional subsurface imaging and automated detection of reflection hyperbolas both reveal two distinct, spatially continuous patterns extending throughout the entire scanned zone at two different depths. Examples of two depth slices partly capturing these patterns are presented in Figure 8. Results of the automated detection are presented in Figure 9, again for quite restrictive threshold values. As at Sippenaeken, most of the detected hyperbolas are located between 0.5 and 2 m.
The first continuity pattern is clearly the most obvious with both methods and appears at an approximate depth of 0.7 m (Figures 8a and 9). Because it has a very straight shape, appears at the same depth all along the scanned zone, and is almost located in the thalweg of the study site, this feature was interpreted as an artificial drainage pipe even though no record of such structure had been made previously (Verachtert, 2011).
The second continuity pattern, well captured by the automated detection method, was more complex to reveal with the 3D subsurface imaging technique (Figure 8b): visualization of many successive depth slices was needed to capture its shape. It also extends throughout the entire scanned zone, but its approximate depth varies between 0.8 and 1.3 m and it presents rapid changes in terms of direction, behaving closely to what is expected from a natural soil pipe. In general, it is located relatively close to the previously described linear structure with which it appears to overlap across a distance of 20 to 30 m on the depth slices. Automated detection of refection hyperbolas nonetheless clearly show that the two detected subsurface structures, although they follow exactly the same direction on this part of the scanned zone, do not merge: the potential artificial drainage pipe is clearly located above the potential natural pipe. In some depth slices, this potential natural pipe also seems to be branching out with Vadose Zone Journal F I G U R E 5 Map of a potential pipe network at Sippenaeken combining three-dimensional subsurface imaging and automated detection of reflection hyperbolas. Validation trenches and piping-related features are also presented on the map shorter secondary axes, but they were less visible and connectivity was difficult to ascertain. All these secondary axes could not be discerned unambiguously with the automated detection method, even when using less restrictive threshold values. Finally, as was observed on the GPR data for Sippenaeken, the potential pipe is sometimes interrupted by large patches of much stronger reflections, especially in the lower part of the scanned zone, which may correspond to backfill material although no sign of past collapses could be inferred from the surface at these locations.
In addition to the previously described features, the 3D subsurface image revealed other important subsurface reflectors, often much stronger than the continuous ones but seemingly randomly distributed throughout the entire scanned zone. Because their signature on the radar data is very similar to the ones generated by filled collapses, they may correspond to past piping activities or may be connected to either smaller and/or deeper pipes impossible to detect with the GPR antenna used. As at Sippenaeken, other potential axes are visible but they are limited in size (maximum 5 m long), isolated from each other, and do not present any obvious signs of connectivity.
Given the present results, only the two clear continuity patterns discussed above were used for pipe network mapping. As previously done with the radar data for Sippenaeken, to lower the uncertainty associated with the less obvious axes, visual analyses of consecutive 2D radargrams were also performed. For the secondary axes branching from the main potential pipe, these analyses allowed ascertaining neither their spatial continuity nor their connection to the main axis because of the large amount of information contained in the signal. The results are presented in Figure 10.

Validation
Several validation trenches were dug along the main axes of the network map ( Figure 10). The validation results were the following: • Trenches no. 1 and 2: Two trenches were dug in the lower part of the potential pipe because this part behaved singularly on the radar data. Indeed, it presented patches of strong reflections linked to each other by weaker reflections, but these patches did not coincide with any surface indicator of pipe collapse. Two trenches were thus dug along these strong reflection patches and, after 40 cm, backfill material was found in both trenches. As at Sippenaeken, the backfill was composed of soil mixed with pieces of anthropic materials (bricks, etc.) that hampered further digging. This backfill material is very likely linked to active piping activities because its water content was much higher at the time of digging than the rest of the field, where soils were relatively dry at comparable depths. As observed at Sippenaeken, it was hypothesized that pipe flow could find its way through this coarse backfill material, ensuring the continuity up-and downstream of the backfill.
• Trench no. 3: This trench confirmed the presence of the artificial drainage pipe and natural soil pipe at different depths. As predicted by the radar data, a modern 8-cmdiameter PVC drainage pipe was found at an approximate depth of 0.6 m (Figure 11a). At the 1.2-m depth, a flowing soil pipe was found but it had developed around a much older artificial drainage pipe (in terracotta, not functional anymore) as presented in Figure 11b (the pipe was 15 cm in diameter). It appears that the soil pipe may have bypassed its natural path once it encountered this old artificial drainage pipe and followed it for 25 m before continuing on its natural course. This probably explains its very straight shape in this part of the subsurface image. • Trench no. 4: After 1 m, a pipe was found. The validation trench was submerged by 30-40 cm of water in a few tens of minutes, suggesting a relatively high discharge in the pipe (Figure 11c). This submersion made any further measurements of the pipe dimensions impossible, but it was nonetheless considered as a validation of the pipe presence. • Trench no. 5: The flowing pipe was found at a depth of approximately 1 m and it was about 8 cm in diameter ( Figure 11d). • Trenches no. 6 and 7: These two trenches were dug in the main secondary axis captured by the 3D subsurface image. Two small dry pipes were found in these trenches but their depth did not match with our expectations from the radar data (40 cm in the field against approximately 60-80 cm on the radar data). Moreover, they did not show any sign of past water flow and rather resembled small animals burrows. They were therefore not considered as a validation of pipe presence.
As mentioned above, numerous piping-related features were visible at the soil surface within the scanned zone. As shown in Figure 10, some of them match perfectly with the pipe positions all along its way (this is the case for the unfilled pipe collapses but also for some filled ones), therefore confirming our validation results. However, many of the filled collapses are not related to any important continuity pattern in our radar data. This could mean that these collapses are related to past piping activities and the pipes are now partially or totally infilled. It could also mean that we were not able to detect the pipes with the adopted radar configuration because they were too small, especially if there was no water flowing inside them during the radar campaign (lower electromagnetic contrast). Finally, similar to the filled collapses in trenches no. 1 and 2, some strong reflection patches on the radar data did not match with any piping-related features visible from the surface. This may suggest that there are more collapses than could be inferred from the surface.
Overall, validation results at Kluisbergen were promising. Indeed, even though the most discernible structure on the Vadose Zone Journal F I G U R E 7 Graphical representation of the subsurface structures detected with the ground-penetrating radar in Sippenaeken radargram was an artificial drainage pipe, it was possible to identify a small pipe (diameter of 8-17 cm depending on the validation trenches) and capture its path extending for >100 m at an approximate depth of 0.8-1.2 m. A tracer test (using milk) later demonstrated that this pipe is the one discharging water in the lower part of the field. Despite these positive results, we were not able to confirm any secondary pipes, and some clear surface indicator of piping activities remained unexplained. Given these validation results, a schematic representation of the confirmed subsurface structures is presented in Figure 12.

DISCUSSION
Results of both field campaigns suggest that the current methodology is promising for noninvasive characterization of pipe networks. Indeed, for the two scanned zones, the high-resolution 3D GPR surveys, combined with the automated detection of reflection hyperbolas, was able to capture the paths of small pipes extending over tens of meters. Nevertheless, these results also highlighted the current limitations of GPR when applied to heterogeneous field conditions with little prior knowledge of pipe location, which allows us to formulate practical proposals to address them.
To characterize pipe networks using GPR, high spatial resolution scanning is crucial. This was already suggested by , who could not confirm pipe connectivity between consecutive transects when using a low spatial resolution. It is even more important under highly heterogenous conditions and with relatively small pipes such as encountered at the two study sites. Indeed, it is impossible F I G U R E 8 Depth slices [approximate depth of (a) 0.6-0.8 m and (b) 1-1.2 m] illustrating reflections from various types of features at Kluisbergen. Color scale represents the reflection intensity on the basis of a single transect to ascertain the presence of a soil pipe without any prior information about its location. Moreover, very often, the strongest reflectors along each individual transect did not correspond to a pipe but rather to local heterogeneities in the soil profile, which can lead to substantial misinterpretation. In a sense, our results suggest that high spatial resolution scanning greatly reduces the ambiguity about pipe presence that may arise from single GPR transects, as previously encountered by Verachtert (2011) and Bernatek-Jakiel and Kondracka (2016). Comparison between the results obtained at the two study sites also suggests that the higher the resolution, the easier the interpretability: at Kluisbergen, using a similar 200-MHz center-frequency antenna but with a higher scanning resolution (one transect every 25 cm) allowed better and easier definition of a pipe having, at least in some parts of the scanned zone, comparable dimensions with the ones detected at Sippenaeken (one transect every 50 cm). This need for high spatial resolution is partly explained by the sudden changes in orientation that these natural pipes may present: lower resolutions may drastically limit continuity assessment. However, it implies some important drawbacks regarding its applicability. It is very time consuming, and many sites affected by piping may not be suitable for high spatial resolution scanning because obstacles (such as trees and roots in piped forests) and degraded landscape (collapses, pipe-related gullies, etc.) can considerably limit its feasibility.
Another important outcome is that, although high spatial resolution scanning is key when dealing with rather small pipes under heterogenous conditions, 3D imaging and visualization through consecutive depth slices is tedious and sometimes misleading. In this context, the use of an automated detection of reflection hyperbolas method greatly facilitated the identification of pipe paths. The Mertens et al. (2016) algorithm used in this study demonstrated its ability to work under field conditions where hyperbolas are often ill shaped. However, to clearly capture pipe paths, several Canny edge threshold values had to be tested and used in parallel: in this way, benefit is derived from distinct levels of information going from the more restrictive approaches (few but clear and more certain continuity patterns) to the less restrictive ones (more continuity patterns, but more uncertain due to a higher false alarm rate). In the end, this may be very time consuming depending on soil heterogeneity and the size and number of transects (in total, several weeks of computing time were needed for both scanned zones). This limitation partly comes from the chosen center frequency because this algorithm has been shown to better work with higher frequencies (higher detection rate and lower false alarm rate) (Mertens et al., 2016).
Apart from the clear continuity patterns revealed at both study sites, the two scanned zones presented short continuity patterns that appeared either isolated from each other or drowned in a large number of reflections. Although it was not possible to perform a validation for these features, it is possible that these correspond to more heterogeneous pipes, in both dimension and orientation. For example, we observed in another part of the Sippenaeken study site a pipe of about 10 cm in diameter located at a depth of about 90 cm that had expanded into a cavity 30 cm wide, 60 cm high, and 1 m long and then contracted again to a 10-cm-diameter pipe with a siphon. In this context, the small continuous patterns could correspond to locally larger pipes and some of the strong isolated reflectors that appear on the images could correspond to local pipe cavities. Weaker continuous patterns may also result from sudden changes in pipe orientation with respect to the scanning direction. Scanning in both directions, perpendicularly and parallel to the slope, could greatly alleviate this orientation issue but it will also significantly increase the time needed to both collect and process the data.
Higher center-frequency GPR systems should also be investigated as they may allow the identification of smaller pipes. However, using higher frequencies implies the detection of more heterogeneities and, if pipe reflections are not strong enough (for instance if they are small and empty), they F I G U R E 1 2 Graphical representation of the subsurface structures detected with the ground-penetrating radar in Kluisbergen could remain embedded in the surrounding heterogeneities and thus impossible to capture. Because using higher frequencies also implies a lower penetration depth, it could furthermore compromise the detection of the deeper pipes. In the future, operating at different frequencies may allow complex pipe networks to be better reconstructed as it should allow the detection of both smaller pipes (if they are not too deep) and deeper ones (if they are large enough). It also means that some pipe network configurations may exceed the GPR detection capabilities: to a certain extent, small and deep pipes will never be detected. Therefore, development of other approaches and novel technologies needs to continue to rise to the challenge of pipe network characterization. This missing link is of crucial importance to better understand the implication of piping in hydrological and geomorphological processes.

CONCLUSION
In this study, two 3D GPR surveys were performed in loessderived soils in Belgium (in Sippenaeken and Kluisbergen) to characterize small natural soil pipe networks with little prior information about their location. The results suggest that the adopted methodology, combining 3D subsurface imaging and automated detection of reflection hyperbolas, is appropriate for this purpose. They also demonstrate that high spatial resolution scanning is essential for its successful application, while relying on an automated object detection method is key to limit misinterpretations and facilitate pipe network mapping under heterogenous conditions. This study also points out some key recommendations for future applications and development. The use of multifrequency GPR systems may, for instance, allow the detection of smaller pipes without compromising the detection of the deeper ones. The refinement of the Mertens et al. (2016) algorithm for lower frequencies could also greatly reduce the time needed to process the data and ease its interpretation. Finally, enhancing the spatial resolution by scanning in two orthogonal directions could greatly improve the characterization of pipe networks, whose abruptly changing directions still restrict the performance of the proposed methodology.
To conclude, this study highlights the potential of GPR for pipe network investigations. However, its requirements are such that it will not be suitable for all piped areas, and other methods need to be developed for forested and highly degraded landscapes where high spatial resolution scanning is not feasible.

ACKNOWLEDGMENTS
We would like to thank the anonymous reviewers for their valuable comments, which greatly helped improve the paper.
We are grateful to the Agentschap voor Natuur en Bos (ANB, Flemish Government) and to Lambert Pinckers for granting access to the sites of Kluisbergen and Sippenaeken, respectively. Special thanks to Prof. Jean Poesen (KULeuven) for directing us to the Kluisbergen site and providing valuable background information, in particular regarding piping activity at this site. We are also grateful to Oscar Driessen and the late Wim van Geel for introducing us to piping at Sippenaeken and for their continuous support throughout the research. Finally, we thank Dr. Albéric Decoster, Bastien Notredame, Pierre André, and Sébastien François, who helped with field work at both study sites.