Evaluating the joint use of GPR and ERT on mapping shallow subsurface features of karst critical zone in southwest China

The soils and underlying weathered carbonate rock in karstic regions play an important role in the infiltration, storage, and retention of water and nutrients. Because of significant heterogeneity of the karst, the use of individual geophysical techniques is often not sufficient for unambiguous assessment of the irregular distributions of soils and underlying fractures. In this study, ground penetrating radar (GPR) and electrical resistivity tomography (ERT) are jointly used with additional observations to delineate the shallow subsurface structure in two exposed profiles. The results show that ERT is effective for detecting the soil–rock interfaces, even for irregular terrain and fracture structures, such as a funnel‐shaped doline, as the soils and rocks show a large resistivity contrast. Although ERT may be able to sense the presence of extensive fracturing, it cannot detect individual small aperture fractures. Joint use of different frequencies of the GPR antenna (e.g., 100 and 500 MHz in this study) allowed the detection of most fractures at different depths in the study sites. However, forward modeling of typical weathered rock features illustrates that the GPR data cannot resolve any reflection signals of the vertical fractures, so the features of vertically enlarged fractures filled by soils cannot be seen from the GPR images. Moreover, large uncertainties of resistivity at the interface between soils or fractures and bedrock limit the identification of an irregularly distributed subsurface structure. Despite the limitations of individual techniques, the combination of ERT and GPR enhances the delineation of the soil–bedrock interface and identification of the fracture network, which can allow an enhanced geological interpretation of shallow subsurface features in the karst areas.


TAO ET AL.
Vadose Zone Journal development of karst features, such as solution-widening fractures, grooves, cavities, and conduits, as well as the surface caving and depressions. Soils are generally 30-50 cm thick and are sporadically developed on carbonate rocks. The composition and structure of shallow subsurface soils and limestone fractures control infiltration and percolation, erosional rates, and patterns in the landscape. Therefore, capturing the shallow subsurface structures can improve our knowledge of the complex hydrodynamic functioning of both unsaturated and saturated zones (Bakalowicz, 1995;Ford & Williams, 2007;Goldscheider & Drew, 2007;Mangin, 1975;White, 2007).
Geophysical methods, such as electrical resistivity tomography (ERT), ground penetrating radar (GPR) and refraction seismics, have been widely used to survey subsurface structures in karst areas. Compared with classical hydrogeological methods (boreholes and pumping tests), geophysical methods can be applied to survey karst terrain and geological features over a large area. Seismic surveys can provide relevant information in karst media (Šumanovac & Weisser, 2001), such as detecting epikarst, and mapping karst near-surface heterogeneities (Bosch & Müller, 2001Guérin & Benderitter, 1995;Ogilvy et al., 1991;Turberg & Barker, 1996), but this technique provides a limited resolution compared to ERT and GPR.
ERT provides two-(2D) and three-dimensional (3D) images of the variations in electrical resistivity (inverse of electrical conductivity) using electrodes typically placed on the ground surface. Electrical resistivity tomography can be effective for detecting cavities, sinkholes, and shallow conduits or enlarged fractures filled by material that provides a large resistivity contrast with respect to the host material (Ellis & Oldenburg, 1994;Guérin et al., 2009;Šumanovac & Weisser, 2001;Valois et al., 2010;van Schoor, 2002;Zhou et al., 2002;Zhu et al., 2011). However, ERT suffers from resolution limitations; consequently, it has mostly been used to locate the upper part of a sedimentary buried karst system and characterize the overlying sedimentary covering (Carrière et al., 2013), as the soil resistivity is typically much lower than the carbonate rock resistivity. Increase of electrode spacing enables ERT to sense the resistivity deeper. However, large electrode spacings result in decreasing spatial resolution and thus use of the ERT image may fail to identify fractures positions.
The GPR technique uses electromagnetic (EM) waves to detect contrasts in electrical properties in the subsurface. When an EM wave encounters a fracture in a solid rock, reflection and refraction occur. Changes in the direction, phases, and amplitudes of wave propagation can be used to quantify properties of individual fractures (e.g., fracture aperture and fracture filling; Deparis & Garambois, 2009;Tsoflias & Hoch, 2006). The polarization properties of EM wavefields can be used to identify steeply dipping fractures due to the phase

Core Ideas
• Using synthetic and field data, we evaluate ERT and GPR for detecting subsurface structures. • We evaluate different GPR antenna frequencies for characterizing rock fracture zones. • We compare ERT and GPR methods for delineating the soil-rock interface and rock fracture zones.
difference between orthogonal pairs of polarization data sets (Tsoflias et al., 2004). Ground penetrating radar typically has a higher resolution than ERT; however, the resolution of the acquired GPR data and the penetration depth of GPR waves are dependent on the frequency of the antenna and the electrical conductivity of the subsurface, respectively. There is a trade-off between resolution and penetration depth (e.g., higher-resolution data but shallower depths of investigation profiles, or lower-resolution data but deeper depths of investigation). Furthermore, the ability of GPR to detect fractures can be limited due to unfavorable fracture orientation, the presence of fracture areas that are smaller than the size of the first Fresnel zone, and limited penetration depth (Dorn, 2013). When GPR is used to identify and locate deep subsurface karst features (e.g., cavities, conduits, and fractures), it often fails to give information about the material forming the structure (Orlando, 2013). The joint use of GPR and ERT can be effective for an enhanced characterization of geological features in karst (Elawadi et al., 2006;El-Qady et al., 2005). The combination of GPR and ERT surveys has been used to identify bedrock (Diallo et al., 2019), and gypsum deposits in urban areas (Gołębiowski & Jarosińska, 2019), and geological structure of karst unsaturated zone (Carrière et al., 2013), and to assess the risk of subsidence of a sinkhole collapse (Gómez-Ortiz & Martín-Crespo, 2012). Most of these studies show that the GPR method has advantages in the imaging of vertical and inclined fractures near the surface, and that the ERT method has advantages in delineating horizontal structures. Because of the strong heterogeneity of the subsurface, the choice of adequate methods for characterizing heterogeneities in the karst environment is very challenging and remains mainly site related (Chalikakis et al., 2011). In the karst areas of southwest China, there are different combinations of soils and rocks and fracture networks that results in contrasts between fill soils or fractures and the surrounding rocks. Until now, there has been lack of quantitative assessments regarding the effectiveness of GPR and ERT for detecting different subsurface structures, and limited studies on whether the joint use of GPR and ERT can improve detection of the strongly heterogeneous subsurface structures in the karst areas of southwest China. F I G U R E 1 Location and geology of the study area (modified from Cheng et al., 2019) The objective of our study is to assess effectiveness of ERT and GPR (and their joint use) for identification of the subsurface structure, including the soil-bedrock interfaces and fracture networks, in the karst region of southwest China. To reduce interpretation ambiguity, the relative dielectric permittivity of the materials is derived from measurements on four typical karst profiles. The capacity of GPR with high and low frequency antenna, and ERT to resolve typical rock features is assessed by using synthetic and field data. The efficiency and complementarity of the joint use of GPR and ERT surveys for the karst interpretation are evaluated by comparison with the detailed visual surveys in two typical exposure profiles.

THE SELECTED SUBSURFACE PROFILES
The test site is located in the Puding County of Guizhou Province, in the center of the southwest China karst terrain. The area has a humid, subtropical monsoon climate, with an annual mean temperature of 15.2˚C. The mean annual rainfall is 1,315 mm, with 85% falling during the wet season (May-October). The lithology is ∼90% Triassic argillaceous limestone and dolomite ( Figure 1). Soil properties and thickness are closely related to rock composition and topography. More than 90% of soil is less than 0.4 m thick in the study area. Exposed rock usually occurs in limestone areas and steep hillslopes, and thick soils are distributed in the valley and plains. The main soil type is red clay formed by carbonate after its solution in hot, humid, and rainy climate conditions (Zhou et al., 2012). In such climate conditions, soils typically have high water content and high clay content (over 10%) (Zhou F I G U R E 2 Two test profiles and their digitalized features of the soil-rock interface and fractures for Prof-ID1 and Prof-ID2. The vertical and horizontal axes indicate distance in meters et al., 2012). Given these properties, the resistivity of soils in the study area is relatively low.
In this study, we selected two profiles (Prof-ID1 and Prof-ID2 in Figure 2), representing two typical carbonate rocks (limestone and dolomite, respectively), in the region. These two profiles are located on exposed excavations adjacent to roadsides, allowing the features of soil thickness and fracture distributions to be manually measured and then digitized as shown in Figure 2. Prof-ID1 is located in Maguan town, Puding County (Figure 1). The field survey shows that the subsurface zone is unconsolidated and concentrated in the shallow part of the profile (about 5 m thick). The uppermost soils are irregularly distributed, within about 1-m depth from the ground surface. The underlying fracture zone consists of horizontal and vertical fractures. The horizontal fractures are produced in the bedded limestones, whereas vertical fracture structures represent reduced dissolution kinetics as the widths decrease from the upper funnel-shaped dolines or grikes to the deeper fractures ( Figure 2). The bedrock layer contains micropores during the genesis of the carbonate rock belonging to the middle section of Guanling Formation (T 2 g 2 ).
The Prof-ID2 profile is dominated by dolomite ( Figure 1). The soils are covered with weeds and are relatively thick, mostly filled in the three funnel-shaped dolines. The underlying rock contains numerous sloping fractures with dip direction of 160˚and angle of 20˚.
Both Prof-ID1 and Prof-ID2 were measured by ERT and GPR in mid-May, 2017, when the soils were relatively dry after a non-rainfall period lasting more than a week. Air temperature was about 25˚C during the survey period. The GPR measurements were executed immediately after the ERT measurements.

The ERT method
Electrical resistivity tomography is an active source geophysical method, using two pairs of electrodes (dipoles) in contact with the ground: one is used to create an electrical field, and the other pair is used to measure the voltage difference from the source. By carrying out such measurements with different geometrical configurations, it is possible to assess the resistivity of the subsurface (Binley & Slater, 2020). The field measurements were carried out using a Syscal Pro 96 (Iris Instruments) (Figure 3a), which can connect to 96 electrodes and collect 10 measurements on dipoles simultaneously.
In this study, we adopted a typical 2D imaging configuration in the field. We installed the electrode array a short distance (0.6 m) away from the exposed face in order to ensure that the observed face is a reasonable match to the image zone (Prof-ID1 and Prof-ID2) ( Figure 2) A total of 48 electrodes were spaced at 0.3 m for Prof-ID1 and 0.5 m for Prof-ID2. However, as the dipoles are separated in an ERT survey, the footprint of the measurement increases and thus there is some inevitable impact of resistivity variation orthogonal to the line (i.e., away from the face). Different electrode configurations (the geometry of the quadrupole) are possible (Binley, 2015). The dipole-dipole mode (Binley, 2015) is most effective for assessing lateral variation in resistivity and was, therefore, adopted here.
Although the measurements were collected in a 2D collinear electrode array, a 3D finite element mesh is needed for forward modeling in order to account for the close proximity of the electrode array to the exposed face. The unstructured 3D finite-element mesh of Prof-ID1 model is shown in Figure 4, and it is similar for Prof-ID2. The measurement errors were estimated based on the difference between forward and reciprocal measurements. According to the previous study by Cheng et al. (2019) using ERT, the relative errors of the measurements are estimated to be 1.6% for Prof-ID1 and Prof-ID2.

3.2
The GPR method

Data collection and processing
In this study, the MALA Ground Explorer (GX) high-dynamic range (HDR) system with shielded antennas and unshielded antennas was applied for measurements ( Figure 3b-d). The MALA system used has antennas with three different frequencies (100, 500, and 1,200 MHz). Fractures can generally be envisaged as layers embedded in a homogeneous rock formation. This gives rise to two signals with opposite polarities reflected by the two sides of a fracture. The time elapsed between the transmission of the signal and its return back to the receiving antenna after reflecting in the subsurface (twoway travel time) is measured and later converted to depth.
High precision x-y-z geolocated data were acquired simultaneously with the free run data using a Trimble GNSS R8 differential GPS, allowing the corrections of GPR profiles for topography and XY positions. The ReflexW GPR and seismic data processing software (Sandmeier, 2015) was used for GPR data processing. The typical sequence of processing steps applied to the collected data includes static corrections (move start time) to adjust time zero at the ground surface, subtract mean (dewow) to remove signal drift or direct current (DC) shift caused by very low frequencies, gain functions, and filter functions. Migration is not applied in the data processing if obvious stratification characteristics and diffraction phenomenon cannot be identified after the above data processing steps. Additionally, topographic correction was applied to replace all traces in their exact location by using the data collected with the  (Chlaib et al., 2014;Di Prinzio et al., 2010;Reynolds, 2011;Telford et al., 1990)  global navigation satellite system (GNSS), and to visualize the topography on the final processed profile.

Physical properties of materials
The most important physical properties of the materials for GPR and resistivity data interpretation are the relative dielectric permittivity (or dielectric constant) and the electrical resistivity, respectively (Diallo et al., 2019). Table 1 shows examples of these properties for some common materials from the investigation of the other sites. Table 1 displays large variations of these properties since the parameters depend on the porosity of the medium, fractures, water content, and the conductivity of the water, mineralogy, pressure, dissolved ion content, temperature, and clay content (Lopez & Gonzalez, 1993;Sbartaï et al., 2007). In order to reduce uncertainty, it is better to measure the physical properties of the relevant materials. For proper conversion of the time axis into the depth axis to evaluate the depth of GPR reflectors, the wave velocity also needs to be assessed. This average velocity can be estimated by many ways, such as performing common mid-point (CMP) surveys, using time-domain reflectometry, measuring travel time between two boreholes, or using buried objects of known depth (Neal, 2004). In this study, the method of known burial depth was used to determine average GPR wave velocity in four typical karst profiles (Prof-a, Prof-b, Prof-c, and Prof-d in Figure 5). The profiles Prof-a to Prof-c are located nearby Puding city, and Prof-d is located upstream of the Puding catchment ( Figure 1). In each profile, boreholes were drilled horizontally at a certain depth and an iron bit was lowered in the borehole as a target body. The high-frequency 1,200-MHz shielded antennas were used to detect the propagation velocity of EM waves with and without the iron bit in the boreholes. A tape measure was used to determine the vertical height (H) between the iron bit and the profile surface (Table 2) for each borehole. The post-processing software (ReflexW) was used to process the radar detection data measured before and after drilling and the two results were compared ( Figure 5). The result of the comparison shows that the yellow point (the iron bit in the boreholes) was located at the highest point of the hyperbolic image of the isolated target body (the red circle). The travel time of radar EM wave at the target body can be obtained from the radargrams shown in Figure 5. The propagation velocity of EM waves in each of these four profiles can be inferred from where V is the EM wave velocity and T is the two-way traveltime.
The estimated propagation velocity, V, ranges from 0.099 to 0.159 m ns −1 with an average of 0.12 m ns −1 (Table 2). This average value was used in GPR data interpretation of the Prof-ID1 and Prof-ID2. The propagation velocity of EM wave is related to the medium's relative permittivity, ε r , according to where c is propagation velocity of the EM wave in a vacuum (∼0.3 m ns −1 ). According to the measured V, the relative permittivity ε r is 1 for air, 18 for wet soils filled in grikes and fractures, and 8.3 for limestone and dolomite ( Figure 5).

THEORETICAL MODELING OF TYPICAL SUBSURFACE PROFILES USING ERT AND GPR
A controlled study was performed prior to the field survey for the purpose of identifying reflection patterns in models representing characteristics of the study sites. Three model profiles were developed representing soil-filled grikes, inclined fractures, and layered structures (see Figures 6 and 7). Inversion and modeling of ERT and forward modeling of GPR were carried out for detection of the near-surface features.

Theoretical inversion by using ERT
For the ERT modeling, an intuitive open source software for complex geoelectrical inversion and modeling (ResIPy) (Blanchy et al., 2020) was used for 3D and 2D inverse modeling. The codes use an unstructured finite element mesh, allowing modeling of the theoretical apparent resistivity and representation of complex geometries. The electrode configuration chosen in this study was dipole-dipole configuration, with an electrode spacing of 0.3 m (Figure 6a-c). The 3D forward model data were perturbed with 2% Gaussian noise and then inverted using the 2D and 3D inverse codes. The 2D inverse modeling was applied as it is much simpler and more conventional in use but does not recognize the rock-air interface adjacent to the electrode array. The resistivity of soil and bedrock was set to 20 and 1,000 Ωm, respectively, following the range of the values in Table 1. The inverted synthetic models for the thin soil layer overlying grikes are shown in Figure 6d and e show the inversed resistivity distribution using the 3D and 2D functions of the ResIPy software, respectively. Both images clearly show the irregular distribution of the soil-rock interface. The 30and 40-Ωm contours (the white lines in Figure 6d and e, F I G U R E 5 Calibration of the relative permittivity and velocity of radar electromagnetic wave and the corresponding radar interpretation results of Prof-a-d without and with the iron bit (the yellow dot inside the red circle) respectively) capture well the irregular soil-rock interface of 3D and 2D inversion results, respectively. Three-dimensional ERT inversion provides a clearer interface between soils and rocks, indicated by a narrow band between the blue to red color for the image in Figure 6e.
A synthetic model consisting of a thin soil layer overlying a sloping fracture is shown in Figure 6b. The ERT tomography can capture the upper soil-rock interface, but it cannot clearly resolve the sloping fracture (Figure 6g). The layered fractures of limestone imbedded with rock fragment are shown in F I G U R E 6 Synthetic study for the electrical resistivity tomography (ERT) method. Modeling of three typical subsurface features that represent: (a) soil-filled grikes; (b) a thin soil layer overlying bedrock with an inclined fracture; (c) layered fractures. Inverted resistivity is shown in Panels d-g; Panels d and e use a three-(3D) and two-dimensional (2D) ERT inversion, respectively, for the synthetic model (a); Panels f and g show that for the synthetic model (b and c, respectively), using a 2D ERT inversion. The white line in Panels d and e shows 30 and 40 Ωm contour, respectively; the black line shows the true interface Figure 6c. As the fragments were surrounded by air, the resistivity was set to a high value (in this case, 2,000 Ωm). The inversion result in Figure 6g reflects the high resistivity zone due to the fragments, but it cannot identify the layered fractures.

Theoretical modeling by using GPR
The forward model of GPR used is based on the time domain simulation software GPRSIM developed by LAUREL (Goodman, 1994). GPRSIM is a diffraction model based on

F I G U R E 7
Forward modeling of three typical subsurface features that represent a thin soil layer overlying grikes (Prof-e), a steep sloping fracture (Prof-f), and layered fractures (Prof-g), respectively. Note: the relative dielectric constants of the corresponding materials are based on Table 1; the bottom of the figure represents the selected electromagnetic wave propagation type, respectively physical optics and ray tracing method based on geometric optics. The software uses the finite-difference time-domain method to obtain the numerical solution of Maxwell's equations. The software allows the use of a customized geological model, and simple setting of model parameters, which displays ray paths and other results. The modeled space was discretized into a grid with a resolution of 0.01 m. A time step of 0.0195 ns was used. The physical properties of the materials, such as the relative permittivity value ε r , were selected according to the measured values in Table 1.
The EM waves traveling through the subsurface encounter a buried discontinuity separating materials of a different physical properties, there various combinations of wave transmission (T) and reflection (R), depending on the properties and shape of the deposit off which they are reflected. As shown in Figure 7, for two horizontal layers, the "TRT" represents transmission (T) into the first layer, reflection (R) off the firstsecond boundary, and transmission (T) back to the surface. For further details, refer to Goodman (1994).
According to the forward modeling results in Figure 7a, the grikes can be identified in terms of the strong reflection signal of the soil-rock interface. The width of the grikes' surface can be identified as the signal segments between two adjacent reflection signals. The bottom of the grike is identified as the reflection segments below the grikes' surface. However, the grikes' side wall reflection is extremely weak, because the side wall is in the vertical direction. Therefore, the GPR image does not reveal any signals of the grikes' side wall, as shown in Figure 7a, but the connection lines of the signal segment terminating points between two adjacent reflection signals perfectly match the side walls.
The main characteristics of individual fracture identification (Figure 7b) are the rapid and regular variation of radar reflection waveform frequency, inconsecutive lineups, and strong signal amplitude. The reflection wave of the sloping fracture is clear and continuous, and its amplitude is obviously stronger than that in the solid rock area. As expected, GPR detects well the slightly slanted stratification, such as the three identified layers with the strong reflected signal of the EM wave (Figure 7c).

5.1
Prof-ID1 Figure 8a shows the 3D resistivity model  interpreted from inversion of the ERT Prof-ID1 profile data. The images based on ERT data interpretation represent ground surface elevation, elevation of top of rock, and soil thickness in the study area. The ERT can survey about 4-m depth below the profile surface using the configura-  Figure 8a shows clear demarcation of the soil-rock interface according to resistivity variations. The much lower resistivity areas (e.g., < 190 Ωm, in blue color, bounded by the solid black contour line in Figure 8a) corresponds to soil due to the presence of moisture and high clay content, whereas the high-resistivity areas (>700 Ωm demarked by the dotted white contour line in Figure 8a) most likely represents the intact rock. We interpret resistivity values between these thresholds (>190 and <700 Ωm) to infer rock that is intensely fractured. These resistivity thresholds are comparable with the resistivity values typically reported for limestone rocks (Table 1).
Although ERT detects the fractured rock characterized by low values of resistivity as a result of moisture presence, it cannot identify distributions of horizontal and vertical fractures shown in in situ measurements. For the GPR survey (Figure 8b, c), 500 and 100 MHz antennas were used. As the vertical resolution of GPR is a quarter of the wavelength of the radar wave, and the wavelength is inversely proportional to the frequency of measurement, the 500 MHz antenna gives a high-resolution image in the vertical layers, whereas the 100 MHz antennas gives a low-resolution image. The high-resolution image (Figure 8b) shows fractured rock properties in the depth less than 4 m, which is much shallower than the identified depth (8 m) from the low-resolution image ( Figure 8c). Unlike ERT and the forward simulation results from GPR (Figure 7), the GPR image of Prof-ID1 (Figure 8b, c) is ambiguous when used to interpret the grike structure filled with soils as the GPR image does not reveal any signals of the side walls of the grikes (Figure 7a). Moreover, many soils containing high gravel content (i.e., fragments of limestone and dolostone) observed in the study profiles produce interference signals (Wang et al., 2015), which can mask signals related to the soil-bedrock interface . Nevertheless, fractures and layered structures can be identified from the amplitude intensity, frequency variation, and phase continuity of GPR. For the joint fractures, the reflection wave represents inconsecutive lineups and obviously stronger amplitude than the intact rock area. Thus, we can decipher the joint fracture distributions in Prof-ID1 shown in Figure 8b and c. For the layered structure, the reflected wave represents the continuous in-phase axis of the signal, the uniform waveform distribution and the strong signal amplitude. A sketch of the layered structure is shown in Figure 8b and c.
The GPR results highlight many fractures within the limestone which are undetected by ERT, as reported by Carrière et al. (2013). However, the GPR results are not useful for detecting the soil-rock interfaces as in the case of ERT. Combining the advantages of the ERT and GPR interpretations, we can depict the structural feature diagram of the profile Prof-ID1 shown in Figure 8d, which is generally consistent with the digitized the structural feature of Prof-ID1 in Figure 2.

Prof-ID2
The ERT and GPR results for Prof-ID2 are shown in Figure 9. Figure 9a shows the resistivity model  interpreted from inversion of the ERT Prof-ID2 profile data. The low resistivity areas (<190 Ωm, the black contour line in Figure 9a) represents the upper funnel-shaped dolines filled by soils. Whereas the high resistivity (>190 Ωm) most likely represents the fracture zone. With the chosen ERT array and inter-electrode space array, it is again to resolve the soil-rock interface but is not possible to detect thin fractures. The GPR results highlight many sloping fractures (the yellow lines in Figure 9b) within the dolomite which are undetected by ERT. Based on the above interpretations of ERT and GPR methods, the structural feature diagram of the profile Prof-ID2 can be developed, as shown in Figure 9c, which is generally consistent with the digitized features of Prof-ID2 in Figure 2.

DISCUSSION
The synthetic and field examples shown above reveal how GPR and ERT may be effective in mapping shallow subsurface features of karst. However, the methods still have some limitations. One of the uncertainties from GRP and ERT modeling arises from the large ranges of resistivity and EM propagation velocity of the soils, fractures, and solid rocks (e.g., several orders of magnitude, as shown in Table 1). The wide range of resistivity values for the detected materials can lead to uncertainty in identifying the interface between unconsolidated materials and bedrock. Concerning GPR, in our study area, the measured propagation velocity within limestone ranges between 0.099 and 0.159 m ns −1 (Table 2) for the four selected rock profiles ( Figure 5). Use of the average velocity (0.12 m ns −1 ) to derive GPR images of the two test profiles ( Figure 2) could result in errors in delineating the subsurface structure. As shown in Figures 8d and 9c, the sloping fractures cannot be exactly captured by the GPR images. For ERT, there are large uncertainties of resistivity at the interface between soils and bedrock, as the porosity, saturation, and gravel may vary differently in the soils-bedrock interface. As shown in Figures 8d and 9c, although ERT can generally capture the soils-bedrock interface at our study sites using a value of 190 Ωm, departures between the in situ measured interface and the inverted interface still exist. Furthermore, as demonstrated by the synthetic models, even though the resistivity of soil and bedrock is known (e.g., 20 and 1,000 Ωm, respectively), the inverted interface resistivity from ERT is 30 and 40 Ωm for 3D and 2D inversion, respectively (Figure 6d, e), which is larger than the assigned soil resistivity because of the inherent smoothing. Concerning GPR, rock fractures in the subsurface typically have apertures less than a wavelength (λ) of the dominant frequency of the GPR signal. When the thickness of thin beds is smaller than the resolution limit, distinguishable anomalies may be lost and the "resolvable limit" is reached (Hosseini, 2014). For example, the Rayleigh resolution limit is λ/4. As shown in Figures 8d and 9c, the imaging cannot capture some vertical and inclined fractures. The study by Markovaara-Koivisto (2017) has shown that it is possible to estimate the fracture aperture when the aperture is wider than the vertical resolution of the antenna, For example, the resolution of a 800 MHz antenna allows detection of a 1-cm-wide waterfilled openings of crystalline rock fractures.
The ERT and GPR interpretation can be constrained with a priori information, such as from borehole measurements, which may help in reducing uncertainty and provide accurate and high-resolution interpretations (Hosseini, 2014;Kana, 2016;Obi, 2012). Prior knowledge on the nature of the rock under investigation, especially propagation velocity, will also help improve GPR modeling (Kana, 2016).

CONCLUSION
In situ explanations of the surveyed results are challenging because of the high heterogeneity in karst weathered medium and limited direct observations. The existence of exposed faces as field laboratories and theoretical modeling reveal how resistivity imaging may be effective in revealing localized infill of soil in karstic environments and how radar reflection imaging may be effective in characterizing of fracture distribution and stratified structure. All geophysical methods produce uncertainty in data interpretation. This can be a result of the ambiguity inherent in data inversion, the nature of signals generated in the subsurface using the particular method, variation in measurement support volume, and ambiguity between inferred geophysical properties (e.g., electrical conductivity and permittivity) and the quantity of interest. In particular, the presence of multiple sources of noise from materials which are not dominant or inhomogeneous (e.g., soils containing high gravel content) can obscure GPR and ERT signals. Inaccurate identification also arises from limitations caused by the resolution of the antenna for GPR and the electrode spacing for ERT, as well as limitations of the forward and inverse modeling of GPR and ERT data, respectively. As shown in this study, the 500-MHz antenna of GPR gives a high-resolution image that can detect detail fractures in the shallow layers (e.g., < 4 m) (Figure 8b). By contrast, a low-resolution image from the 100-MHz antennas can only detect the fractures and layered structures in the deep layers that interference signals represent obviously stronger amplitude than the intact rock area (Figure 8c). The theoretical modeling and inversion by using GPR and ERT also suggest that GPR signals cannot be directly used to visualize a vertical fracture wall and ERT cannot identify individual fractures.
The joint use of GPR and ERT is effective for providing an enhanced characterization of geological features in karst media. In this study, ERT is effective for detecting the shallow funnel-shaped dolines or enlarged fractures filled by soils since the ERT image provides a large contrast in resistivity of soils with respect to that of the rock. Joint use of different frequencies of GPR antenna (e.g., 100 and 500 MHz in this study) can be used to detect effectively most fractures underlying the soil and determine fracture features including joints and fractured rocks with specific inclinations. Therefore, a combination of ERT and GPR can fully delineate the soil-bedrock interface and identify fracture features.

A C K N O W L E D G M E N T S
This research was supported by the National Natural Science Foundation of China (42030506), and the UK Natural Environment Council (NERC) Grant NE/N007409/1 awarded to Lancaster University. We are grateful for the comments received from the Associate Editor and two anonymous reviews on an earlier version of the manuscript.

C O N F L I C T O F I N T E R E S T
We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.