1 Introduction

The purpose for taking measurements on or beneath the earth’s surface is to gain an understanding of the subsurface conditions especially for the detection of subsurface voids such as washouts or settlement gaps, floor slabs, tunnels, vaults, caves which are natural limestone (karst) or solution cavities, and sinkholes. These voids which are naturally occurring or man-made are of geological and engineering interests. This is because they are potentially dangerous geological hazards [5]. Their presence might cause the weakening of foundation materials and in some cases lead to structural defects/failures/collapses. There are a number of methods/techniques available for voids detection [5, 8, 31, 32, 42, 54]. A combination of gravimetry and electrical resistivity imaging was used to identify the subsurface voids over the cave area [27]. Microgravity surveys were carried out for the investigations of natural cavities at mines [12, 39]; monitor of collapse of salt mines [7]. Elsewhere, surface wave diffractions were used for voids and faults detection [60]. The capability of using ground-penetrating radar (GPR) for mapping underground discontinuities due to faults, boulders, or pipes has been undertaken to gain a better understanding of the subsurface conditions for environmental, archaeological, hydrogeological, and engineering applications [5, 9].

The choice of appropriate geophysical tools for the detection of subsurface voids is not very easy to carry out due to the unpredictable nature and variations in the characteristics of the subsurface anomalies [57]. Also, getting to the depth of investigation and the resolution of the targets are others significant factors too [42]. The utilization of geophysical methods particularly the electrical resistivity in mapping buried objects (anomalies) arising from good electrical contrasts between the buried objects and the host surroundings. Electrical methods are vital geophysical tools for subsurface imaging and detection of structures such as sinkholes, cavities, and voids (natural underground chambers). These structures can either be air-filled (empty) or saturated with fluids [31]. A single geophysical method might not be adequate for subsurface voids detection due to their nature and expected depth of probe [57]. Furthermore, the ambiguities associated with the interpretation of most geophysical methods make a single technique inappropriate. To address the ambiguities, the application of multiple techniques or methods such as satellite imageries (remotely sensed data), borehole log data (ground truth) to constrain interpretation has been well acknowledged in the literature [22, 63].

Geotechnical investigations are essentially invasive methods which are mostly destructive and point specific. They are carried out to understand the dynamic physical properties of the subsurface structures [30, 45]. The applications of geophysical and geotechnical techniques for subsurface investigations provide vital information with respect to the geometrics (i.e., depth, size, and shape) of the subsurface geological structures, particularly in mapping cavities. Several studies have been carried out on subsurface mapping and characterization using both geophysical and geotechnical techniques [23, 24]. However, with the exception of [33] very limited studies have been reported on the use of geotechnical techniques for voids mapping. Against this backdrop, the geotechnical techniques are used to complement the geophysical techniques.

In this study, the two methods were used to map the presence cavity in the subsurface. First, the area investigated was delineated into a number of layers using variations in resistivity with depth from the VES measurement. Then, the VES data were laterally constrained to give two-dimensional (2D) resistivity models. Second, drilling tests were carried out from which resistance to penetration into the different soil strata was recorded in terms of SPT N-values. Also, sieve tests on representative soil samples were carried out. Also, accuracy assessment using uncertainties and coverage factor arising from the measured parameters was determined.

2 Location and geology of the study area

The study area, Klebang Putra at Kinta Valley District Perak in Peninsular Malaysia, lies between Latitude 4°41′8.20″N–4°41′10.70″N and Longitude 101°6′57.70″E–101°6′47.20″E. Geologically, the study area as shown in Fig. 1 is underlain by two Formations; the Carboniferous with Kinta limestone/sedimentary rocks and acidic undifferentiated igneous rocks which are mostly granitic [51]. Since the Kinta District is underlained by a sequence of sedimentary rocks ranging in age from Silirian to Permian, there are evidences of granitoids intrusion of Jurassic to Triassic ages [28]. Also, alluvium has been reported to cover the entire valley part of the district [28, 43]. In a related work, Tan [53] showed that Kinta District is underlain by an extensive limestone bedrock Formation called Kinta limestone. In situ measurements through boreholes or excavations revealed that subsurface geology features and soil profiles consist mainly of limestone bedrock and its associated karstic with the depth to the limestone bedrock put at less than 20 m [52]. The prevalence of limestone in this area shows the possibility of cavities in the area. According to Yassin et al. [61, 62] could be encountered at depth between 1.0 and 15 m. The alluvium which covers the floor of the Kinta Valley consists of Recent Clastic sediments obtained from weathering and erosion of the neighboring rocks. Karstification in the valley arises from the dissolution of soluble rock bodies, leading to formation of ridges, fissures, and drainage system [59].

Fig. 1
figure 1

A simplified geology map of Perak with Districts (modified after [61])

Full size image

3 Theoretical background

3.1 Basic concepts of Electrical Resistivity Technique

In vertical electrical soundings (VES) technique, the subsurface is assumed to consist of a sequence of layers of finite thicknesses; each of these layers is assumed to be electrically homogeneous, isotropic, and the boundaries between subsequent layers are horizontal. The method is based on the application of a linear filter to determine the apparent resistivity curve from the kernel function [16]. The process includes calculating resistivity transforms and apparent resistivity values which facilitate the interpretation of resistivity soundings. A typical VES model modified after Auken et al. [3] as part of initial model in the design of two-dimensional (2D) subsurface model is shown in Fig. 2. It can be described with model parameters (i.e., layer resistivities, ρ, and depths to layer boundaries, d). The models are connected with their neighbors by laterally constraining them with respect to their model parameters (Wisén et al. 2005). Figure 2a shows a 2D resistivity pseudo-section obtained from inversion of continuous vertical electrical soundings (CVES) data using any of the 2D inversion programs [13, 56]. As shown, the resistivity model is divided into cells (represented by different colors) by a process called model discretization. To achieve the laterally constrained inversion (LCI), it is necessary that all the VES models (Fig. 2b) should have the same model parameters and are simultaneously constrained using appropriate algorithm. Then, inversion is performed to find an optimal model by minimizing an objective function ([2]; Auken et al. 2005).

Fig. 2
figure 2

Laterally Constrained VES models (modified after [3])

Full size image

The apparent resistivity and depth of each layer across the boundary were used to simulate the resistivity responses of various geological conditions. The process of resistivity modeling is expressed in terms of set of equations that characterize each layer through the forward modeling.

The procedure for the resistivity modeling can be estimated by an integral equation of the form in Eq. (1).

$$d_{i} = \int_{0}^{z} {K_{i} (z) p(z){\text{d}}z}$$
(1)

where di is the measurable response due to electric current injection into the ground, p(z) is a resistivity function related to the Earth’s structure as a function of depth in a laterally homogeneous Earth. It is referred to as the model parameter which is a continuous function of position. Ki are called the Kernels function of the data. The Kernels describe the relationship between the data and the earth model function p(z). Forward modeling is used to predict the data or responses that would be recorded over a hypothetical earth structure. For a two-layer Earth model using the Schlumberger electrode array, the apparent resistivity according to [38] is given by Eq. (2).

$$\rho_{2 } = \rho_{1} (1 + 2L^{2} )\int_{0}^{\infty } {K(\lambda ) J_{i} (\lambda L)\lambda {\text{d}}\lambda }$$
(2)

where

  • L is the half electrode separation,

  • Ji is the Bessel function of ith order,

  • K(λ) is a function transforming resistivity model parameters,

  • (λL)λ is the lambda coefficient which determines the electrode parameter,

  • ρ1 and ρ2 are the resistivities of the first and second layers.

The K(λ) is given by Eq. (3).

$$K(\lambda ) = \frac{{ - K_{1,2 } \exp ( - 2\lambda t)}}{{1 + K_{1,2} \exp ( - 2\lambda t)}}$$
(3)

where

  • t-factor is depth of layer

  • $$K_{1,2} = \frac{{\rho_{1} - \rho_{2} }}{{\rho_{1} + \rho_{2} }}$$

  • K1,2 is the resistivity contrast between the two layers. It ranges between − 1 and + 1 as the resistivity ratio

  • (ρ2/ρ1) for an infinitely conducting and resistive bottom layer, respectively, varies between 0 and infinity.

Since it is not easy to express Eq. (3) in the form d = Gm (i.e., forward model as a Kernel function that describes the mapping of the model parameters into the dataset), the resistivity–depth relationship is regarded as a nonlinear problem [26].

Inverse modeling on the other hand is the process of looking for the best geological model of the earth from a given measured field data. The inversion process starts with datasets which contain the measured resistivity for a given electrode configuration and ends with an electrical image of the subsurface. Geophysical inversion suffers from non-uniqueness because there are several models that can explain the same datasets. One of the ways to reduce non-uniqueness is to use additional information from other sources to constrain the inversion process [50]. To find the model that best fits the original measurements in terms of model misfit, a regularized optimization technique is introduced into the inversion process [55].

To find the model that best fits the original measurements in terms of model misfit (i.e., least root means square error), a regularized optimization technique is introduced in the inversion process. The objective function which is minimized in solving the inverse problem by optimization is given in Eq. (4).

$$\varphi (m) = \varphi_{d} (m) + \alpha \varphi_{m} (m)$$
(4)

where

  • φ(m) is data model misfit;

  • φd(m) is Chi-squared measure of the data misfit;

  • φm(m) is model objective function containing the desired model characteristics;

  • α is regularization parameter;

  • m is variable model parameters.

$$\varphi_{d } (m) = \left\| {W_{d} \left[ {d - f(m)} \right]} \right\| ^{2}$$
(5)

The first term on the right-hand side of Eq. (5) is a measure of the data misfit, with f denoting the forward operator, Wd represents the data weighing matrix associated with the individual (uncorrelated data errors). In addition, φ(m) contains a stabilizing model objective function usually expressed as Eq. (6)

$$\varphi (m) = \left\| {W_{m} [m - m_{f} ]} \right\|^{2}$$
(6)

It is used to incorporate certain model constraints relative to a reference model (mref) by suitable choice of a model weighting matrix, Wm The regularization parameter, α, in Eq. (4) controls the tradeoff between the data misfit and the model objective function.

3.2 Geotechnical investigations (in situ field tests)

In geotechnical investigations, soil properties are determined by both laboratory and in situ testing since they provide complementary information about the subsurface conditions. The standard penetration test (SPT) is an in situ measurement designed to provide information on the geotechnical engineering properties (i.e., both static and dynamic) of the soils. The dynamic soil properties can be determined using an active source of energy to either excite the soils or induce a measurable seismic wave [46]. Examples of the dynamic soil properties are soil moduli [i.e., shear modulus (Go), Young modulus (E) or bulk modulus (B)], damping ratio (D), and Poisson’s ratio (υ). For some of these soil properties, the effect of disturbance on the soils is very important. As a result, correlations using empirical relations or regression models have been established to determine these properties particularly from in situ penetration tests [19, 40, 47]. The SPT value indicates the density of ground in terms of number of blows or N-values. SPT due to its simplicity and low cost, is a popular in situ tests for measuring resistance to penetration of soils in engineering applications [14, 36, 48]. ASTM has become a standard for SPT classification of soils for subsurface geotechnical investigations. Although SPT was initially designed for coarse grained soil, it can also be used for estimating the properties of fine grained soil such as undrained compressive strength, undrained shear strength, and coefficient of volume compressibility [48] as well as soil classification [41]. SPT were carried out through boring or in situ measurements of the number of blows (N-values).

The particle size distribution (PSD) involves the measurement of the different sizes of particles of subsoils such as clay, silt, and sand as determined by their abilities to pass through sieves of various mesh size or by their rates of settling in water. The various sizes of the subsoils were used as diagnostic characteristics in the classification schemes as they provide information on the soil’s PSD [6]. In geotechnical engineering, knowledge of PSD of soils can be significant in providing an initial rough estimates of soil engineering behaviors like compressibility, shear strength, permeability, expansibility, and soil classification [21, 35]. Also, the PSD can be useful in understanding soils physical and chemical properties [10]. A variety of methods are available for measuring individual particles in the terms of particle size, size distribution, and shapes. Examples of common methods include sieve analysis; sedimentation [11, 34]; microscopic sizing [25, 58]; image analysis [4, 18]; and laser diffraction methods [15, 49]. The sieve technique which is one of the oldest methods is adopted in this study because of its simplicity, accuracy, cheapness, and easy of interpretation [20].

3.3 Overview of uncertainty and coverage factor

Uncertainty at certain degree of confidence is usually associated with every measurement. The result of a measurement is only complete provided the uncertainty that may arise due to a number of factors from the measurement is stated. Many of the techniques for estimating uncertainty are based on analyses resulting from the perturbation of the original data sets [37]. Statistical analysis, however, has proven to be a vital tool in estimating uncertainties in measurements. When uncertainties in measurements are evaluated, a good judgment with respect to fitness for the purpose of the measurement is arrived at. In order to estimate combined standard uncertainty in measurements, rescaling of the components of uncertainty becomes necessary. Rescaling, however, can be carried out via coverage factor (K). The standard uncertainty demarcates an interval where it is possible to estimate say about two-third of the results obtained. Thus, the coverage factor (k) is always given together with the uncertainty. Although a combined standard uncertainty (uc) is used to express the uncertainty in several measurements, however, in geophysical surveys parameters what is often required is a measure of uncertainty that defines an interval about the measurement result y within which the value of the measurand Y can be confidently asserted to lie [17]. The expanded uncertainty (U) is obtained as product of coverage factor (k) and the combined uncertainty, uc(y) as in Eq. (7) [44].

$$U = k*u_{c} (y)$$
(7)

where

  • U = Expanded uncertainty

  • k = Coverage factor

  • uc(y) = Combined uncertainty

Traditionally, a particular value of coverage factor gives a particular confidence level for the expanded uncertainty with the overall uncertainty scaled by coverage factor k = 2, to give a 95% of confidence level.

4 Material and methodology

In this study, the geophysical method used for the electrical resistivity measurement was the VES technique. Also, the geotechnical investigation was carried out through field test in which boring with standard penetration tests (NSPT values) and laboratory tests in which particle size distribution (PSD) analysis were carried out on some of the soil samples collected. Then, the verification of the earth model reconstructed was carried out using statistical technique via coverage factor (K) and uncertainty values (U).

4.1 VES measurements

The resistivity measurements were carried out along five traverses established within the Klebang Putra district with five VES stations spaced at 10 m to 50 m from one another as shown in Fig. 3. For the field measurements, the half-spacing of current electrode (AB/2) used ranges from 1.0 to 50 m. The profiles crossed some of the available boreholes for comparison. The soundings measurements were carried out to have an insight of the variation of resistivity with depths in the subsurface of the study area. In principle, the VES measurements involve the injection of current into the ground through the two current electrodes A and B. Then, the resulting voltage (i.e., potential difference) was measured by the other two potential electrodes M and N. The Wenner-α electrode configuration was adopted for the resistivity measurements using ABEM SAS 4000 model resistivity terrameter. The interpretation of the field apparent resistivity data was obtained quantitatively using curve matching technique to generate the model parameters. Then, an inversion of data sets was carried out using IPI2WIN software version. The input model parameters in the IPI2WIN program were the electrode spacing (AB/2) m, the layer thicknesses, the number of layers, and the resistivity values of each layer. They were used to search for subsurface optimal resistivity model during the inversion process by an iterative process. The apparent resistivity curves were analyzed using the inversion technique to obtain the true resistivity. The principle of ambiguity in the inverse problem solution—there are many possible resistivity models having the same or similar solution. The reconstructed resistivity and thickness values were constrained with the borehole data to provide information on the subsurface geology.

Fig. 3
figure 3

Location map of the Klebang Putra in Kinta Valley showing the VES and Borehole points for SPT test

Full size image

The VES measurements along the profiles were used to generate the 2D LCI resistivity models through the concatenate option from RES2DINV [29]. This allows the resistivity variations with depth of investigations to be imaged with lateral constraints and sharp boundaries. The data sets are inverted as one system, producing layered solutions with laterally smooth transitions. The 2D LCI inversion was performed with equal size constraints on both the layer resistivities and depths to layers across the boundaries. For the resistivities and depth, constraint factor of 1.1 and 1.3, respectively, was used [2].

4.2 Geotechnical investigation (field test)

Five boreholes (BHs 1–5) were drilled to allow for the control of the basic inference made on the resistivity models. To this end, the SPT N-values and the PSD measurements were carried out. The SPT-N values were obtained in accordance with the ASTM D- [1] standard from which the number of blows made by the sampler was recorded. The samples were obtained at 1.5 m intervals into the ground surface by driving the sampler with a 50 mm diameter thick-walled tube called the split spoon sampler into the hole using a 63.5 kg hammer. For the PSD tests, nineteen samples were collected from BHs 1, 2, and 4, while eleven samples from BHs 3 and 5 for sieve analysis. The sieves test consists of a wire mesh or slotted metal plate that permits particles smaller than the mesh size to pass through. The various sizes of the clay, silt, and sand were used as diagnostic characteristics in the classification schemes.

4.3 Model validation

Statistical analysis was used to validate the field data measurements using the coverage factor, k = 2, and uncertainty values, U, at 95% confidence level. Furthermore, a tool bar applied in Dinver program was used to analyze data from 1D geo-electric measurements on a single piece automatically or semi-automatically to get the smallest error and vertical lithological variations along with their thicknesses, depths, and resistivity values.

5 Results and discussion

The VES results are presented as resistivity sounding curves in Fig. 4a–e. A summary of the inverted resistivity, depth, and inferred lithology units is presented in Table 1. The depth sounding curves are generally of K-type (ρ1 < ρ2 > ρ3) with resistivity rising to a maximum and then decreases, indicating that the middle (intermediate) layer has higher resistivity than the overlying (top) and underlying (bottom) layers. Three geo-electric layers are delineated. They are the topsoil, weathered bedrock, and sandy clay/clayey sand. The soundings curves reveal the presence of a relatively high resistive zone with resistivity of about 500 Ωm and low resistive substratum with resistivity less than 100 Ωm. The high resistivity layer is the weathered bedrock (limestone with cavity). This is underlain by a layer with low resistivity between 72 and 120 Ωm. It is presumptuous of weathered bedrock but filled with fluid or sand/clayey sand sediments. The 2D LCI models with the distribution of resistivity are shown in Fig. 5. The high resistivity region represents air-filled (cavity) weathered/fractured bedrock denoted by red color band at lateral distance of 110–130 m. The resistivity is about 430 Ωm and at a depth of 10 m from the ground surface. Also, a low resistive structure with resistivity between 40 and 164 Ωm (blue color) is prominent in all the sections. The inverted resistivity section also shows a possible weathered bedrock (cavity) which is characterized with lower resistivity of about 120 Ωm which could be due to clay infill or groundwater accumulation. The SPT profile for the five boreholes is shown in Fig. 6. Low N-values between 0 and 5 are obtained for layers of the earth classified as weak zones or unconsolidated topsoil, clay or clayey silt soils while relatively high SPT-N values between 10 and 20 show that dense materials are encountered at greater depth into the subsurface strata. The dense materials are suggestive of consolidated sediments typically sand layer or bedrock.

Fig. 4
figure 4

1D VES inversion curves showing variation in resistivity with depth

Full size image
Table 1 A summary of interpreted resistivity with thicknesses, depths, and inferred lithology
Full size table
Fig. 5
figure 5

2D inverted electrical resistivity (ERI) sections

Full size image
Fig. 6
figure 6

Subsurface profile based on SPT N-values and BHs at test site

Full size image

The result of the PSD tests examined on some representative soil samples plotted on a semi-logarithmic graph presented as curves of cumulative percentage as function of soil particle’s diameter for the five boreholes is shown in Fig. 7. It is observed that the PSD curves are consistent in all the boreholes both near the surface and at greater depths. The position and shape of the grading curves determine the soil classes to which the samples belong. The different soil grades obtained from the gradation sieve analysis are clay, silt, sand, and gravel. At BHs1, 2, and 4, the gradation analysis at depths between 3.5 and 6.5 m for the PSD curves indicates clay/silt sediments with sizes smaller than 0.075 mm. As a result, the clay/silt sediments are of little or no importance in geo-engineering applications. For BHs 3 and 5, at depths between 3.0 and 4.5 m, the PSD curves indicate sand/gravel sediments with sizes greater than 0.075 mm which is of greater value in examining the competence of the subsurface conditions to supports massive engineering structures. Thus, the PSD curves show that more than 80% of the soil samples passed the 75 μm sieve. Based on percentage of samples that passes through the sieve, the soils are classified as clay, silt, sand, and gravel.

Fig. 7
figure 7

Particle size distribution curves for samples from BH1-5

Full size image

The results of coverage factor (K) and uncertainty values (U) at 95% confidence level are shown in Fig. 8. The histograms with the uncertainty analysis are generated using the IP2WIN software. The difference in the models created on the histograms is due to the changes in the electrical properties of the earth layers delineated. The uncertainty values are as a result of the deviation in model parameters such as resistivities, thicknesses, and the upper boundary depths between the theoretical and observed field data. It is a reflection of the results obtained in field tests through SPT N-values, PSD, and BHs analyses. In all the boreholes, the histograms are characterized with low percentage errors between 0.14 and 0.33%. Figure 9 shows a correction of the geophysical and geotechnical results. The methods showed good agreement with respect to the soil layering information revealed.

Fig. 8
figure 8

Bar charts with uncertainty analysis along traverse a one for BH1 b two for BH2 c three for BH3 d four for BH4 e five for BH5

Full size image
Fig. 9
figure 9

Correlation of geophysical and geotechnical results

Full size image

6 Conclusions

A combination of geophysical and geotechnical methods was undertaken to gain knowledge of the subsurface conditions at Kinta Valley District, Malaysia. The VES technique delineated three geo-electrical layers; topsoil, sandy clay/sand, and weathered bedrock. The 2D resistivity sections from the laterally constrained VES models revealed the location and extent of a geological structure which is characterized with high resistivity (air-filled cavity) to about 10 m from the surface at lateral distance of 110–130 m. At this zone, the SPT N-values were relatively low, suggesting little or no resistance to penetration. While for deeper portion of the sections, fairly high N-values were obtained. This region is the limestone bedrock. The uncertainties associated with the techniques statistically lies within the intervals of 0.14 and 0.33% for coverage factor k = 2 and 95% confidence level.