Supervised Statistical Optimization of Interfaced Geophysical–Geotechnical Datasets for Holistic Surface–Subsurface Modeling

1. Introduction

In tropical crystalline basement terrains, characterizing the surface and underlying soil-rock characteristics poses significant challenges due to the complex lithologies and the absence of distinctive stratigraphic markers for quick delineation of soil-rock conditions (Kalantary et al., 2009).. Detailed spatial variability characterization of soil-rock profiles is crucial for geotechnical, geoengineering, and hydrogeophysical evaluations. Various geophysical (electrical and seismic methods) and soil test techniques have been developed for this purpose (Priya and Dodagoudar, 2015; Abudeif et al., 2019; An et al., 2022; Sujitapan et al., 2024).

Advancements in geophysical imaging have led to improvements in the aforementioned approaches. The most widely used near-surface geophysical methods for characterizing soil-rock profiles include electrical resistivity tomography (ERT) and seismic refraction tomography (SRT). SRT measures seismic P-wave velocity (Vp), while ERT determines changes in the georesistivity of subsurface properties (Chambers et al., 2006; Quigley, 2006; Akingboye and Ogunyele, 2019; Loke et al., 2022). Both techniques have proven valuable in numerous near-surface characterization projects compared to intrusive sampling techniques, which provide discrete point information (Gao et al., 2018; Hasan et al., 2022a; Akingboye and Bery, 2023a; Husainy et al., 2023). They are effective for small- and large-scale measurements due to variations in soil-rock properties and geohydrodynamics. To fully optimize the applicability of both methods, it is best to treat them as joint velocity–resistivity relationships, utilizing statistical correlations and analytical modeling (Meju et al., 2003; Zeng et al., 2018). However, these techniques have some weaknesses that affect their performance.

Over the years, important soil-rock engineering parameters such as rock mass quality (RMQ) and soil penetration test (SPT N-values) have become standard for detailed evaluations of in-situ rock mass integrity and soil strength (El-Naqa, 1996; Hasancebi and Ulusay, 2007; Priya and Dodagoudar, 2015; Akingboye, 2023; Akingboye and Bery, 2023b). These parameters can be used to address the limitations of SRT and ERT models. Establishing empirical relationships is a significant method for providing strong theoretical correlations between in-situ soil-rock geotechnical parameters and geophysical data. Standardized statistical analyses of geotechnical parameters (RQD and N-values) with velocity and resistivity data can effectively yield the required results (Sivrikaya and Toğrol, 2006; Priya and Dodagoudar, 2015; Kumar et al., 2016). Such estimates and derivations are crucial in infrastructure design, groundwater exploration, and soil/rock excavation to predict associated physical, geomechanical, and hydrogeological properties in areas with limited borehole information (Sivrikaya and Toğrol, 2006; Kalantary et al., 2009; Al-Heety et al., 2021; Akingboye and Bery, 2022). Additionally, detailed characterization of soil-rock profiles and their interfaces can provide two- and three-dimensional (2-D and 3-D) variability images of the near-surface, essential for mitigating risks like slope failure, soil liquefaction, and ground subsidence (Kalantary et al., 2009; Priya and Dodagoudar, 2015; Sujitapan et al., 2024). This approach can help prevent foundation failure. Similarly, delineating planes of weakness or soil-rock structures from such evaluations is crucial for groundwater prospecting (Al-Heety et al., 2021; Akingboye et al., 2022; Akingboye, 2024).

Moreover, the techniques mentioned above hold significant promise for the detailed characterization of the intricate soil-rock profiles and structures beneath Penang Island, Peninsular Malaysia, which is the study area. This is crucial due to the region's annual torrential downpours, highly diverse residual soil profiles, and varied topographic relief (Bery and Saad, 2012a; Akingboye and Bery, 2021; Mohamad et al., 2022). Additionally, with the yearly growth in tourism and population in the area, there is a pressing need for thorough investigations to support sustainable infrastructure design and groundwater development. Therefore, an innovative integrated approach is proposed, combining velocity-resistivity analysis, RMQ assessment, SPT N-value determination, and regression analytical modeling to comprehensively elucidate the surface and subsurface soil-rock characteristics of the study area. This study represents a significant advancement over the piecemeal information presented in earlier studies on the subject and its environs (Bery and Saad, 2012; Salleh et al., 2021; Akingboye and Bery, 2022).

This study used both 2-D and 3-D (fence diagram) visualizations to delineate and resolve uncertainties in surface–subsurface lithological conditions in the study area. Indirect subsurface Vp and resistivity models were validated using direct borehole measurements, from which the RQD and N-values were estimated. This combined approach yielded satisfactory results. The objectives of the study were to: (i) assess the nature, forms, and variability of soil-rock profiles, including their depths and architectures; (ii) define the geometries and apertures of bedrock weathered troughs and fractures; (iii) determine the surface–subsurface soil-rock quality and conditions; (iv) develop a suite of statistical (empirical) relations between velocity–resistivity data and RQD and N-values for rapid and accurate prediction of these soil engineering parameters, aiming to reduce borehole operational costs; and (v) evaluate the implications of objectives (i–iii) on infrastructure design and groundwater exploitation in the research area.

2. Regional and Local Geology of the Study Area

The Western Sibumasu Terrane and the Eastern Sukhothai Arc (referred to as the East Malaya Block) constitute two conjugated Late Triassic tectonic terranes comprising Peninsular Malaysia (Metcalfe, 1998, 2013; Pour and Hashim, 2017), Figure 1a. Evidence of these tectonic terranes is found along the southeast perimeter of the Eurasian continent. They originated during the Permian-Triassic period through the amalgamation of the Sibumasu and Indochina Blocks at the Bentong-Raub Suture Zone (BRSZ) (Metcalfe, 2001, 2000; Cao et al., 2020). The BRSZ, a deeply rooted, 13 km wide, N-S trending tectonic zone, divides East Malaya and the Sibumasu Block (Metcalfe, 2000; Pour and Hashim, 2015, 2017). This zone is a remnant of the main Devonian to Middle Triassic Palaeo-Tethys Ocean. Moreover, it delineates Peninsular Malaysia into the Western Belt, Central Belt, and Eastern Belt (Ng et al., 2015a, 2015b; Cao et al., 2020), while also giving rise to numerous roughly N-S trending faults (Pour and Hashim, 2015). In Malaysia, granitoids, mainly granite, constitute the major portion of the plutonic rocks. These granitoids predominantly comprise Permian to Late Triassic I- and S-type granites. The Eastern Belt of the Indochina region hosts I-type granites. Meanwhile, the Western and Central Belts of the Sibumasu domain are composed of Late Triassic S-type granites, typically biotite granites and granodiorites (Ng et al., 2015a).

Underlying Penang Island are typical granitic rocks, primarily granites, characterized by an alkali to total feldspar content (Ong, 1993). These granites are categorized into various groups: Tanjung Bungah, Paya Terubong, Batu Ferringhi, and the Muka Head micro granite, collectively forming the North Penang Pluton (NPP). The Tanjung Bungah group, composed of medium- to coarse-grained biotite granite, was formed during the Early Jurassic. The Paya Terubong group's medium to coarse-grained biotite granite originated during the Early Permian to Late Carboniferous period. The Batu Ferringhi group, dating back to the Early Jurassic, consists of medium- to coarse-grained biotite granite. The South Penang Pluton (SPP) comprises the Sungai Ara and Batu Maung granites (Ong, 1993; Mohamad et al., 2022). The Sungai Ara granite, representing 20% of the area, is predominantly microcline granite with fine-grained biotite granite characteristics. On the other hand, medium- to coarse-grained biotite granite constitutes approximately 80% of the microcline granite found in the Batu Maung granite (Mohamad et al., 2022).

3. Methods

The study area encompasses three locations within Penang Island, Malaysia: the Universiti Sains Malaysia (USM) main campus as Site 1, Batu Maung as Site 2, and Jelutong as Site 3 (see Figure 1 and Figure 2). Eleven geophysical traverses (Figure 2 and Table 1) were conducted in this study to analyze the spatial variability of soil-rock conditions and structural architectures beneath the specified sites. Furthermore, in situ soil-rock conditions and integrity were assessed using borehole (BH)-based RQD and SPT N-values. The obtained research findings serve as typical representations of the near-surface soil-rock characteristics underlying Penang Island. The methodological framework of the research is depicted in Figure 3.

3.1. Geophysical Survey Traverses and Field Data Measurements

Sites 1, 2, and 3 are situated in the eastern, southern, and northern parts of the study area, respectively. Sites 1 and 2 lie within the Sungai Ara granitic formation in the South Penang Pluton (SPP). Site 1 is positioned closer to the North Penang Pluton (NPP) as it is situated almost at the northeastern end of the SPP. Meanwhile, Site 3 is within the NPP.

At Site 1, two ERT traverses (TR1 and TR2) were conducted in the Maiden area of USM. The third traverse (TR3), involving both SRT and ERT measurements, was established at the Restu students’ residence in the southeastern part of USM. TR1 and TR2 had ERT traverse lengths of 200 m and 100 m, respectively, with spacings of 5 m and 2.5 m. TR3 consisted of an ERT profile of 100 m, while its SRT line was positioned between stations 17 m and 87 m, covering a profile length of 70 m (see Figure 2a). This configuration enabled the resolution of complex soil-rock structures beneath the line. The lithologic log obtained from boreholes (BH) at 85 m station distance along TR1 at Site 1 was correlated with the ERT models of TR1 and TR2. Similarly, the lithologic log from a nearby existing borehole was correlated with the ERT and SRT models of TR3 to ascertain the true nature of soil-rock profiles and resolve layer thickness.

The second location (Site 2) is situated a few kilometers from Site 1. A total of four traverses, comprising four SRT and two ERT traverses, were conducted. Lines 1 and 2 were established generally in an SE–NW direction, while Lines 3 and 4 ran in a SW–NE direction (refer to Figure 2b). SRT Line 1 had a length of 80.5 m with a geophone separation of 3.5 m. In contrast, SRT Lines 2–4 were each 115 m long with a geophone spacing of 5 m. Both geophone intervals were utilized to capture detailed near-surface soil-rock features. The soil-rock profiles along Lines 3 and 4, generated by SRT, were further analyzed using ERT. Each ERT line measured 160 m in length and had an electrode spacing of 4 m (see Figure 2b). Lines 1 and 2 posed significant challenges for conducting ERT surveys due to the rugged terrain and steepness. To achieve a balance between the dual spacings (i.e., 5 m and 3.5 m) used for SRT at Site 2, a 4 m electrode spacing for ERT was deemed suitable. A total of six boreholes (BH1–BH6) were drilled at Site 2 to facilitate detailed interpretation of the subsurface lithologic units.

At Site 3, four traverses (L1 to L4) were conducted for both ERT and SRT measurements (see Figure 2c). Shorter traverses and station intervals were utilized due to the site's location in a developed area. L1 and L2 had a profile length of 80 m each, oriented approximately in an east-west direction. L3 and L4 had spread lengths of 74 m and 40 m, respectively, oriented in a north-south direction. Except for L4, which had a 1.0 m station spacing, all traverses employed a 2.0 m station spacing (refer to Table 1). One of the two existing boreholes at this site, BH1, intersected L1 and L4 at stations 38 m and 36.5 m, respectively. The second borehole (BH2) is located at a station point of 44 m on L3.

3.2. Geophysical Field ERT and SRT Data Measurements

The ERT surveys utilized the ABEM Lund Imaging Resistivity Meter, comprising the ABEM SAS Terrameter 4000 and the ES 64-10C Electrode Selector, along with their respective accessories. The Wenner-Schlumberger array was employed to measure surface–subsurface resistivity distributions across the study area. For the SRT surveys, the ABEM Terraloc Mk8 was utilized, accompanied by a 12-pound hammer and a channel of 24 geophones operating at 14 Hz. This system facilitated the measurement of seismic P-wave velocity (Vp) values, providing insight into the spatial variability of soil-rock conditions. To generate high-resolution and detailed Vp models, clustered shot points were utilized, with 5–7 shots each at shot points between station intervals and offsets. During data inversions, surface topographic variations along the surveyed lines were extensively corrected using recorded topographic elevations for station positions.

3.3. ERT Data Processing, Inversion, and Modeling

RES2D inversion software was employed to process and iteratively invert the resistivity values measured in the research area, following the methodology outlined by Loke and Barker (1996). The primary objective of data inversion in geophysical interpretation is to generate a detailed model for the subsurface resistivity distribution, ensuring that the calculated model response closely matches the measured apparent resistivity values. However, it's important to note that the solution obtained through inversion is typically non-unique in practical field surveys. Therefore, the standard least-squares inversion technique was utilized, as described by DeGroot-Hedlin and Constable (1990), Loke et al. (2022), and Akingboye and Bery (2023a).

The four-node L₂-norm finite element approach formed the basis for the field data inversion constraints, aimed at minimizing the discrepancy between the calculated and observed resistivities. Additionally, the accuracy and resolution of the inversion model were enhanced by employing a damping factor of 0.05, with a minimum value set at 0.01. The inversion process resulted in resistivity sections for the surveyed traverses converging by the 5th iteration, with inverse convergence root mean square errors (RMSEs) below 10%. The fact that none of the data points were eliminated during inversion serves as evidence of the high quality of the field datasets. Figure 4 illustrates a composite resistivity section showing the observed and calculated pseudosections, along with the inverted resistivity section for TR1 at Site 1.

3.4. SRT Data Processing and Modeling

To process and invert the measured field Vp datasets, observed, recorded, and source data files were compiled, incorporating distances of the shot points and elevations of station points for the measured Vp data points. FirstPix software was utilized to enhance weak seismic traces using the auto gain control method. This technique amplified the first arrival travel times, facilitating easier identification, tracing, and picking. A representative example of the raypath generated for the picked seismic first arrival traces for L3 at Site 3 (Jelutong) is depicted in Figure 5. Subsequently, the resulting data files were loaded into SeisOpt@2D to generate the final Vp model for each of the investigated traverses.

The SeisOpt software utilizes a nonlinear optimization technique for the inversion of the picked first arrival traveltimes. In its forward modeling process, the software quickly solves the eikonal problem using finite differences, as described by Akingboye and Bery (2023b). The data inversion was conducted using the Refraction Inversion and Optimization Technique (RIOT) in high-resolution mode. Compared to medium- and low-resolution modes, the high-resolution mode produces detailed subsurface Vp inversion models with minimal structural artifacts. Subsequently, the final Vp models generated for all SRT traverses were refined using Surfer^TM software. This approach was employed to create Vp models with topographic reliefs and a uniform velocity color scale. Finally, the resulting Vp models were correlated with the borehole lithological logs to facilitate detailed interpretation of surface–subsurface soil-rock conditions, bedrock structures, water-saturated fills, and other relevant features.

It is important to highlight that all analyzed locations were subjected to the same range of Vp values, reflecting a common Vp color scale. This consistent range of the Vp color scale is suitable for all SRT models across the study area, as variations in overburden thicknesses, lithologic unit ages, and densities result in similar distinctive velocity features (Akingboye, 2023). However, it is worth noting that not all models produced by ERT, particularly those for Site 2, are compatible with a single resistivity scale. When a lower resistivity color scale was applied to the ERT models of Site 2, similar to that used for Sites 1 and 3, highly resistive materials without signs of low-resistive bodies were eliminated. Consequently, information about weathered granite, resistant stiff-to-hard soil bodies, and resistive boulders present in the ERT models of Sites 1 and 3 was excluded. Conversely, using a high-end resistivity scale for Site 2 resulted in the exclusion of features that may be indicative of weathered granite, stiff-to-hard soil bodies, and boulders. Therefore, to achieve a comprehensive interpretation, it is advisable to adopt two different resistivity color scales, ensuring synergy between borehole-derived logs and ERT models. Variations in soil-rock minerals, water saturation, porosities, etc., are attributed to observed resistivity variations (Rucker et al., 2021; Akingboye, 2024; Sujitapan et al., 2024).

3.5. Evaluation and Modeling of RMQ of Subsurface Lithological Conditions

To acquire accurate data on the RMQ of lithologic units within the study area, the borehole-based RQD of subsurface lithologic units was assessed (refer to Figure 3). This evaluation was exclusively conducted in the newly drilled boreholes at Site 2, as the existing boreholes at Sites 1 and 3 lacked RQD data. This approach minimizes drilling costs for the research. Interestingly, the RQD values derived from six strategically located boreholes (BH1–BH6) along Lines 1–4 at Site 2 provide representative RQD indices indicative of the RMQ for the entire study area.

For geoengineering purposes, lithologic units (or rock masses) are categorized using the RMQ system, which assesses the integrity or quality of the rocks, the presence of significant fractures or planes of weakness, and the orientation of these planes of weakness. RMQ is vital for determining various factors such as the depth of the foundation, the capacity of the foundation rocks to bear weight, and the potential for foundation sliding (Hoek and Brown, 1997; Barton, 2006; Hoek and Diederichs, 2006; Hasan et al., 2023;).

To evaluate the RQD of the study area’s lithologic units, forty-five core samples (i.e., $n=45$) from BH1–BH6 at Site 2 were estimated and used. A core drilling machine with a double-layer core tube and a diamond 75 mm drilling bit was employed to recover the core samples. In the past, RQD was calculated using only the 54.7 mm (NX-size core) core diameters. Topsoil and completely weathered materials constituting the residual soil profiles do not meet the requirements for RMQ estimation; hence, they were not considered (Barton, 2006; Pells et al., 2017). It is also worth noting that the water table varied between 2 m and 5 m during the drilling operations. The RQD values of each recovered core sample in a borehole, with respect to the core run length, were evaluated using Eq. 1 (Abzalov, 2016). RQD represents the percentage of sound intact core pieces longer than 4 inches (i.e., >100 mm) to the rock’s total core run length.

$$RQD (\%)=\frac{\sum length of sound pieces>100 mm}{Total core run length}\times 100$$

(1)

The resulting borehole-RQD values at various subsurface depths were correlated with their corresponding SRT models to determine their Vp values. Utilizing Surfer software, the RMQ models of the investigated four lines at Site 2 were derived based on the combined results of the borehole-RQD and their Vp values. These RMQ models provided information on the weathering conditions of the foundation rocks, jointing/fracturing patterns, load-bearing capacities, and water saturation state. To optimize the performance of the derived RMQ models, the simple linear regression (SLR) technique was applied to further analyze the borehole-RQD and their Vp values. This process aimed to establish a new empirical relationship specifically tailored for granitic terrains with shallow overburden depths of <50 m. Such an empirical relationship proves especially valuable in regions with limited borehole data and steep terrains where traditional drilling methods are challenging.

3.6. Evaluation and Modeling of N-Values of Subsurface Soil Profiles

SPT serves as a cost-effective geotechnical and geoengineering testing technique commonly employed in infrastructure design (Bery and Saad, 2012; Syukri, 2020; Al-Heety et al., 2021). While traditionally utilized for engineering site assessments, SPT results also offer valuable insights for groundwater exploration. In this study, in situ, borehole SPT N-values were derived for Sites 1 and 3 (refer to Figure 3) to assess the strengths, consistency, and penetration resistance of the soil profiles beneath these sites. Given the overburden thickness and weathering profiles at these locations, Sites 1 and 3 were deemed most suitable for SPT testing, as drilling boreholes for RQD evaluation in these areas would be both challenging and expensive.

To determine the SPT N-values at the studied sites, a total of 35 N-values were derived from four existing boreholes drilled at Sites 1 and 3. These N-values, along with their corresponding depths, were modeled in Microsoft Excel to clarify the spatial variability across the evaluated sites. Subsequently, the 35 N-values were correlated with the corresponding resistivity and Vp models of TR1 and TR3 (Site 1), and L1 and L3 (Site 3), at respective depths. The resulting N-values, resistivity, and Vp values were averaged to obtain 16 datasets based on comparable depths. These datasets were then utilized to develop empirical relationships for N-values with Vp and resistivity using both SLR and multiple linear regression (MLR) approaches.

4. Results

4.1. Borehole Lithologic Logs of the Study Area

The borehole lithologic units from the three studied sites within Penang Island primarily consist of clay/silt to sandy topsoil, silty-to-sandy weathered units, weathered/fractured granite, and integral/fresh granitic bedrock (Figure 6a–e). Site 2 is predominantly characterized by silty sand to sand (Figure 6c), while Site 1 (Figure 6a–b) and Site 3 (Figure 6d–e) comprise sandy silt with minor clay content. Consequently, silty sand emerges as the most prevalent surficial and subsurface soil profile in the study area. The identified soil profiles in these areas range from medium stiff to hard. Interestingly, bedrock was delineated at Sites 1 and 2 in the borehole logs (Figure 6a–c). However, even at a drilling depth of approximately 18 m at Site 3, bedrock penetration was not achieved (Figure 6d–e), indicating severe weathering of the granitic bedrock at Site 3. The thickness of the overburden varies across the study area, with Site 1 (eastern part) exhibiting a greater overburden thickness compared to Site 3 (northern part). Beneath TR1 (Site 1), the overburden, including weathered/fractured zones, extends to a depth of >41 m (Figure 6a).

4.2. Spatial Variability of Surface–Subsurface Soil-Rock Conditions in the Study Area

According to the correlated borehole lithologic logs with Vp and resistivity models of the three sites, the study area exhibits distinct lithologic units comprising topsoil (sandy silt to sand), medium stiff to hard silty-to-sandy weathered units, weathered/fractured granite, and integral/fresh granitic bedrock. The separation of clay/silt from the sand and water-saturated columns in the ERT models was achieved by carefully incorporating the borehole lithologic logs of each analyzed traverse.

4.2.1. Characterization of Soil-Rock Conditions and Architecture at Site 1

The borehole lithologic log (Figure 6a) and the ERT models of TR1 and TR2 (Figure 7a–b) provide insights into the underlying lithologic units at Site 1, along with their geometries, interfaces, and depths. Dry sand dominates the near-surface between stations 17 m and 58 m in Figure 7a, with depths reaching 14 m and the highest soil resistivity around 1800 Ωm. A similar feature can be observed between stations 75 m and 110 m, at depths ranging from 5 m to 17 m, where a central portion divides two low-resistive sandy silt bodies into an oval form, likely representing water pathways (fractures). Figure 7b depicts the sandy silt bodies beneath TR2 (Site 1) with distinct geometries and dipping directions. However, a thin layer of silty sand characterizes depths of 30–31.5 m beneath TR1 (Figure 7c), likely representing deeply weathered feldspar-rich sections of the granitic bedrock, significantly enhancing the residual soil’s thickness.

The lithological units of TR2 (in N–S) are comparable to those beneath TR1, except for the water-saturated area below the topsoil (Figure 7b). Similar characteristics are observed beneath TR3 (Figure 8a–b), suggesting extensive granitic bedrock weathering across the area. The extensive low-resistive body (sand) in the central part of TR2, with <200 Ωm, is saturated and extends to the southern part, likely representing a deep-weathered/fractured conduit for water. It is underlain by resistive dry sand or relatively weathered granite. In general, sandy silt and sand constitute approximately equal proportions of the subsurface soil constituents beneath TR2.

Furthermore, the borehole lithologic log of TR3 primarily consists of silt/clay to silty sand, extending to depths of 20 m (Figure 8c), with low resistivity values of <200 Ωm (Figure 8a). Moderately to relatively high resistivities were observed in the western end and central sections. The top layer to a depth of 10 m exhibits low Vp and resistivity values of <800 m/s and <1000 Ωm, respectively, indicative of the residual soil cover (Figure 8a–b). This section is underlain by the weathered granitic unit with moderate Vp (800–2000 m/s) and resistivity values (<1000 Ωm). Additionally, potential fractures are indicated by the SRT model of TR3 at station positions of 42–47 m and 72–77 m (Figure 8b), although deeper mapping of this structure was impeded by the propagation of weak seismic waves.

4.2.2. Characterization of Soil-Rock Conditions and Architecture at Site 2

Figure 9a–d illustrates the SRT inverted models for Lines 1–4 at Site 2 in Batu Maung. Figure 9e and 9f, respectively, show the 3-D SRT model of the surveyed lines and the borehole logs of BH1-BH6. The borehole lithologic units exhibit a strong correlation with the SRT models of Lines 1–4. The surface–subsurface soil-rock profiles are characterized by residual soils (<800 m/s), weathered and weathered/fractured granitic units (800 to <2000 m/s), and integral (hard)/fresh bedrock (2000–4500 m/s), Figure 9a–e. The ERT models of Lines 3 and 4 at this site are depicted in Figure 10a–b, showing highly varied resistivity values compared to other sites. The delineated soil-rock profiles, including residual soils, weathered/fractured granitic units, and integral/fresh bedrock, exhibit resistivity values ranging from 100 to >1000 Ωm, 900–2000 Ωm, and >2000 Ωm (Figure 10a–b). The presence of stiff to hard sand bodies, interspersed/outcropped granite, and boulders contributes to the high resistivity values of residual soils (such as >1000 Ωm). Except for Line 1, the topsoil of silty sand is thin across the site. Beneath Site 2, bedrock fractures with various axial directions, apertures, and depths are also delineated.

The residual soil’s thickness, as depicted in Figure 9a, ranges from 8 to 17 m. From the northwest to the southeast, the layer dips slightly, reaching a steep slope between station positions of 38 m and 58 m along Line 1. In contrast, the weathered granitic unit is about 6–12 m thick. Between station distances of 45 and 70 m, the thickness of the bulged zone decreases, likely due to soil erosion from the steep sections or fractures. The steep slopes, particularly Lines 1 and 2 (Figure 9a–b), exhibit high soil erodibility, leading to significant sediment load downslope toward the site’s southern half (Figure 9e). Figure 9b shows similar Vp trends and subsurface features to Line 1 (Figure 9a). However, the overburden of Line 2 is thinner and underlain by rugose bedrock. Its residual soil and weathered layers are also approximately 1–12 m and 5–15 m thick, respectively.

Moreover, the weathered layers and residual soils of Lines 3 and 4 have thicknesses ranging from <1 to 11 m and 6 to 18 m, respectively (Figure 9c–d and Figure 10a–b). The granitic bedrock has suffered significant weathering to depths >20 m, as shown in Figure 10, contrasting with the underlying bedrock characteristics in Figure 9c–d. The central section of Site 2 depicts massive bedrock (Figure 9e–f and Figure 10), which plunges farther toward the southeasterly portion (Figure 9e). Interestingly, new fractures were delineated in Figure 10 along with prominent penetrative fractures that conform to those seen in the SRT models.

4.2.3. Characterization of Soil-Rock Conditions and Architecture at Site 3

The soil-rock conditions at Site 3 in the northernmost region of the research area (Jelutong, Penang Island) (Figure 11 and Figure 12) have several distinguishing characteristics from those identified beneath Site 1 (Figure 7 and Figure 8). The soil-rock profiles beneath this part of the study area are composed largely of silty sand with little to trace clay (Figure 11c, Figure 12c,f). However, Site 3 is characterized by generally low resistivity signatures as opposed to Sites 1 and 2 due to the high water retention rate of the residual soils. L1 and L2, in the approximately E–W direction, are characterized entirely by low resistivity values (<200 Ωm). L1 depicts low to moderate Vp values (50–1900 m/s), extending to a shallow depth. Between station positions of 42 m and 62 m, a resistive body that is indicated as fresh granitic bedrock plunges from a depth of 6 m westward (Figure 11a). This feature has high resistivity and moderate Vp values of <2000 m/s (Figure 11b). A similar feature stretching from station position of 33 m eastward was also delineated with moderate resistivity values of 400–800 Ωm at a depth of 8 m.

Furthermore, L2 (Figure 11d–e) is characterized by lithological features that are comparable to those beneath L1. At a depth of 6 m between station positions 30 m and 45 m, a weathered granitic body with moderate to high resistivity values (400–1000 Ωm) was identified in Figure 11d. Contrary to L1 (Figure 11a), where it was separated by a deep-weathered trough with water fills, the feature is still relatively intact. The SRT model (Figure 11e) depicts well-correlated lithologic features to Figure 11d. From the eastern to station point of 44 m, Figure 11e depicts generally flat weathered bodies. From there through the end of the model, the layer, however, undulates. In general, Vp values of >1000–4500 m/s, which range from moderate to extremely high, define L2. Station 42 m from 6 m depth to the western end of the model is characterized by soil-rock materials of high Vp values. This anomalous feature, however, depicts moderate resistivity values, as shown in Figure 11d. This can be categorized as partially weathered to integral granitic bedrock, with a depth comparable to the feature delineated beneath L1.

Figure 12a–c displays the ERT and SRT models (in N–S) and the BH log for L3 at Site 3, while Figure 12d–f displays comparable results for L4 in the same direction. Similar features are shown in the resistivity model of L3 (Figure 12a) as those beneath L1 and L2. A weathered granitic mass with reasonably high resistivity (400–1300 Ωm) that stretches from 40 to 56 m deep was also identified (Figure 12a). This body was not delineated by the borehole along this traverse (Figure 12c). However, it identified stiff to hard silty sand produced from weathered granitic bedrock. The SRT model depicts low to moderate Vp values (Figure 12b). A moderate Vp layer similar to the features in Figure 12a was also mapped at a depth of 7 m in Figure 12b. L4 is also characterized by low to moderate resistivity and Vp values (Figure 12d–e). Despite the penetration depth of the SRT model (Figure 12e), the weathered layer at the depth of about 2 m was delineated. This exact body was likewise penetrated by the borehole (Figure 12f) in that location. Interestingly, the actual orientations of the delineated subsurface soil-rock profiles, as well as their geometries and depths, are shown in the 3-D resistivity and Vp models of Site 3 (Figure 13a–b). The delineated resistive granite at the depths of 6–8 m, reaching the basal section of the model, cut across Site 3. This feature may have been tectonically uplifted at the central part because of the depth of the weathered granitic unit. The deep-weathered/fractured zones at Site 3 include the segmented granitic bedrock in L1 and other traverses.

In summary, the eastern part of Site 1 (Figure 7 and Figure 8) to the northern part of Site 3 (Figure 11, Figure 12 and Figure 13) are characterized by comparable surface–subsurface soil-rock characteristics and architecture. Deep weathering, fractures, and high water retention are features of these areas. However, there are considerable differences in the constituents of the soil-rock profiles beneath Batu Maung (Site 2) in the southern section of Penang Island (Figure 9 and Figure 10). The site is characterized by sandy-rich residual soils, multiple penetrative fractures, and outcropping bedrock/boulders. Particularly at Site 2, the fractures/faults extend to depths >20 m, while the deep-weathered zones at Sites 1 and 3 are more pronounced. Notably, in Figure 1b, linear high hills dominate the study area's central region. The summarized spatial variability of surface–subsurface soil-rock conditions, as identified in borehole-correlated logs with ERT and SRT models of the studied sites in Penang Island, are presented in Table 2. The soil-rock parameters such as compactness, uniaxial strength, porosity, density, degree of void ratio, moisture content, fractures/faults, etc., are thought to be responsible for the observed differences in resistivity and Vp values (Barton, 2006; Priya and Dodagoudar, 2015; Akingboye and Bery, 2022, 2023a).

5. Discussion

5.1. Surface–Subsurface Lithologic Characteristics, RMQ, and Statistical Modeling

To wholly comprehend the RMQ of the lithologic units, we integrated the evaluated borehole-RQD and their Vp values from BH1–BH6 (Figure 6c) at Site 2 with the SRT models derived for Lines 1–4 in Figure 9a–d to develop their 2-D/3-D RMQ models. We statistically analyzed the evaluated borehole-RQD and their Vp values (45 data in total) using the SLR approach, treating borehole-RQD and Vp data as dependent and independent variables, respectively. The developed RMQ regressed plot via SLR modeling (Figure 14) illustrates the trends and the exact range of Vp values for the borehole-RQD data. Additionally, the model generated a concise linear lithology-based empirical relation for rapid and accurate prediction of RQD. Figure 14 provides detailed classifications for the highly weathered to integral/fresh granitic units, which are listed in Table 3.

Via this modeling technique, we inferred five major RMQ ratings based on the quality of the evaluated granitic rock mass. An RMQ rating with Vp of 800–950 m/s and an RQD of 0–25% defines the highly weathered granitic rock unit of very poor quality. The poor-quality weathered granite is rated with Vp and RQD values of 950–1500 m/s and 25–50%, respectively, while the relatively weathered granite of fair quality has Vp values of 1500–1900 m/s and RQD of 50–75%. The very poor to fair quality weathered granitic rock units are rippable in nature, constituting the weathered/fractured granitic units. Conversely, the Vp and RQD values of 1900–2400 m/s and 75–90% depict the integral (hard) granitic rock. The fresh (competent/non-rippable) granitic bedrock is characterized by Vp and RQD values of 2400–4500 m/s and 90–100%. The residual soil profiles across the area are characterized by <800 m/s and are also considered part of the rippable unit in the area. The derived lithology-based empirical equation between borehole-RQD and Vp datasets yielded a high regression coefficient of determination $R^2 value$ of ≈0.96 (Eq. 2). This lithology-based empirical equation can adequately resolve overlapping Vp values for different soil-rock layered conditions.

$$RQD=0.0394Vp-7.5136 (R^2=0.96)$$

(2)

Furthermore, to ensure the accuracy of the borehole-RQD and Vp values with their lithology-based empirical relation, we conducted an accuracy assessment using IBM SPSS (v26) to mitigate statistical bias. This SLR accuracy assessment involved significant statistical tests such as multicollinearity (tolerance and variance inflation factor [VIF]), multivariate normality, homoscedasticity, Durbin-Watson (D-W), and analysis of variance (ANOVA), rather than relying solely on typical statistical estimates like mean, median, mode, and standard deviation.

Table 4 presents the summarized results of the regression accuracy assessment. The correlation (R) value of 0.979 derived for the 45 RQD and Vp analyzed datasets indicates a strong positive relationship between the two parameters. Additionally, the $R^2 value$ of 0.959, with a Vp standard error of 0.001, indicates high prediction strength in the analysis. The $R^2$ result affirms the strong linear relationship between RQD and Vp, explaining approximately 96% of the data. This value translates to a prediction accuracy of approximately 96%. The SLR model also yielded a D-W statistic of 1.724, indicating no autocorrelation between RQD and Vp data. The variables appear to fit well to the regressed model based on the p-value of <0.05%, the large F-statistic of 1015.264, and the collinearity tolerance and VIF of 1.0, implying a model free from multicollinearity issues (Akingboye and Bery, 2022). Hence, the derived lithology-based empirical equation is significant for the direct estimation of RQD from Vp models, particularly in areas with few or no borehole data and inaccessible locations where drilling is practically impossible. This approach is also cost-efficient for use in large-scale geotechnical/engineering investigations. However, the derived empirical relation for RMQ estimation is site-specific, so it should be applied in granitic terrains with comparable lithological characteristics.

Based on the achieved statistical analytical modeling, with certified accurate assessment, the 2-D/3-D RMQ models of Lines 1–4 were developed (Figure 15) based on the RMQ chart produced by Caterpillar Incorporation (2010). The 2- and 3-D RMQ models (Figure 15a–e) of Site 2 distinctly depict the lateral and vertical strength of the penetrated rock masses with their range of Vp values, as well as the planes of weakness (fractures). The study area’s topsoil and completely weathered granite comprise the first layer (i.e., the residual soils), which is understood as the rippable rock unit. This layer’s values are below 800 m/s. The residual soil profile has a thickness of between 1 and 17 m. The highly weathered to relatively weathered granitic units, with Vp values of 800–1900 m/s, are identified as the second RMQ unit. This layer might be as thin as a few centimeters beneath the residual soil or as deep as 20 m. This layer is underlain by the integral (hard) granitic unit, which has Vp values of 1900 to 2400 m/s. With values >2400 m/s, the fresh granitic bedrock is identified as non-rippable granitic bedrock. Figure 15 also clearly displays the subsurface fractures.

The 3-D RMQ (Figure 15e) indicates an increase in the thickness of rippable residual soil towards the eastern and northeastern regions. There are variations in the thickness of the highly weathered to relatively weathered granitic units in the northeastern, eastern, and western parts. However, the fresh granitic bedrock unit extends southward from the north. Despite the thinness of the integral granitic unit, ranging from 2 to 5 m, it plays a crucial role in determining the appropriate section for infrastructure foundations. This layer serves as the boundary between the competent/non-rippable (fresh) granitic bedrock and the rippable soil-rock units.

5.2. Statistical Analytic Modeling of N-Value and Velocity-Resistivity Datasets

SPT N-values serve as indicators of dynamic strength, consistency, and penetration resistance of soils in BH1 and BH2 (TR1 and TR3 at Site 1), as well as BH1 (L1) and BH2 (L3) at Site 3. They are valuable for assessing the spatial variability of surface–subsurface soil conditions, particularly in geologic terrains with diverse and complex soil-rock compositions. Figure 16 illustrates the variation of the 35 evaluated N-values with depth for Sites 1 and 2. These N-value results exhibit a strong correlation with the interpreted Vp and resistivity models. Notably, at depths of 8 m and beyond, there is a significant increase in soil penetration resistance (N-value >35). This value reaches 50 at approximately 9 m depth and remains consistent above 12 m for all SPT boreholes. The depths where the most cohesive layer and the integral granitic rock, both exhibiting N-values of 50, are encountered, are correlated. Additionally, the topsoil layer is characterized by a thickness of less than 2 m and N-values ranging from 7 to 17. Remarkably, the depths of N-values, Vp and resistivity, and RMQ models are closely aligned, indicating the high precision of the models.

The SPT N-value classifications for cohesive soil conditions by Karol (1960) and Terzaghi & Peck (1967), as presented in Table 5, offer detailed insights into the correlation between N-values and soil penetration strength and conditions. Consequently, in tropical granitic terrains, the ranges of N-values for various soil conditions are shown in Table 6 along with their corresponding Vp values. The depths of the entire granitic rock unit and the most cohesive layer, characterized by an N-value of 50, exhibit a strong correlation (Figure 15; Table 6). The N-values obtained in this investigation generally align closely with the depths of the Vp, resistivity, and RMQ models. Although the observed N-values and their corresponding Vp values are slightly lower than those predicted for borehole-RQD, N-values of 50 produced a range of Vp values that correlated well. This variation may stem from the unstandardized relationship between N-values and soil characteristics (Sivrikaya and Toğrol, 2006). Overall, this study successfully demonstrates the relationships between N-values and Vp and/or resistivity data with soil-rock quality.

The statistical empirical relationships for the averaged 16 N-values with their corresponding Vp values, as determined via SLR, are presented in Figure 17. Most importantly, Figure 17 illustrates the distributions and trends of the analyzed N-values and their Vp values. The derived lithology-based empirical relation (Eq. 3) is critical for predicting N-values from Vp data in granitic environments without conducting SPT results. This will significantly reduce borehole operation costs for a similar purpose. Eq. 3 yielded an $R^2$ of 0.943 for N-values with Vp, indicating a 94.3% prediction accuracy. The observed correlation is attributed to increasing soil-rock velocities with depths.

$$Vp=44.392\left(N\right)+182.98 (R^2=0.943)$$

(3)

A unified empirical relation was also generated using MLR for the evaluated N-values utilizing the Vp and resistivity datasets. MLR is a crucial statistical instrument for comprehending indicator variables as a prediction model to examine geophysical data from an investigated site (Tabachnick and Fidell, 2019). To ascertain the empirical link for N-values with Vp and resistivity data, the average Vp and resistivity data acquired were translated into linear forms using MLR. $B_1$, and $B_2$ are the coefficients of $X_1$ (Vp), and $X_2$ ($\rho$), respectively, while $C$ is the equation constant, given in Eq. 4. These parameters were determined via MLR analysis, and their resulting values are presented in Table 6. Equation 5 represents the empirical relation that can be used to predict N-values in granitic terrains using Vp and resistivity data. Table 7 shows impressive statistical results from the MLR analysis, with an $R^2$ of 0.943, low standard error, a p-value of 0.0001, VIF of 1.692, and D-W of 1.938. Normal D-W levels range from 1.5 to 2.5. The violation of no autocorrelation, which can lead to incorrect standard errors for the regression coefficient estimations, perhaps makes D-W values outside of this range (particularly values 1 or >3) of major concern (Akingboye and Bery, 2023a). The estimated p-value for resistivity and N-values is higher than the p-value of 0.05 for regressed data. This explains the variability of the soil-rock resistivities with the determined N-values probably due to increasing water affinity. Generally, N-values and Vp data yielded a strong correlation of $R^2$ of 0.943 (Figure 17), indicating that soil penetration resistance increases with density and depth and that the weathered layer contains probably more sand than the topsoil.

$$N=B_1\left(X_1\right)+B_2\left(X_2\right)+C$$

(4)

$$N=0.021\left(Vp\right)-0.00001033\left(\rho \right)+5.583$$

(5)

5.3. Geoengineering Implications of Velocity–Resistivity, RMQ, and N-Value Results

It is crucial to thoroughly assess the impact of near-surface lithologic units on both pre- and post-infrastructure design to mitigate the risk of foundation failure or collapse. In addition to inadequate foundation design, infrastructure failure can be caused by soil liquefaction, soil failure, differential settlement, ground subsidence, and other factors. These challenges arise from variations in the physical, hydrogeological, and geomechanical properties of surface and subsurface soil-rock profiles (Hasan et al., 2023; Sujitapan et al., 2024) . Hence, a comprehensive evaluation of soil-rock suitability, incorporating Vp, resistivity, RMQ, and N-values models, along with associated regression analytical models (Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 and Figure 15 and Figure 16), is essential for designing resilient infrastructure.

The soil-rock profiles beneath Sites 1 and 3, located in the eastern and northern parts of the study area, predominantly consist of high clayey/silty soil content (<200 Ωm and <1900 m/s) (Figure 7 and Figure 8 and Figure 11, Figure 12 and Figure 13). Moreover, the presence of thick saturated silt/silt sandy weathered bodies and deep-weathered/fractured zones indicates a high affinity for water. These soil-rock conditions, along with delineated outcropping rock (or boulders), highly varied saturated soils, and multiple fractures at Site 2 (Figure 9 and Figure 10), pose significant threats to engineered infrastructure (El Shinawi et al., 2022; Akingboye, 2023; Hasan et al., 2023). Site 2, in particular, is susceptible to increased run-off due to its proximity to extensive mountainous hills. Consequently, excessive rainfall could lead to saturation of steep soils, potentially triggering soil failure.

The northern portion (Site 3) exhibits a higher water retention level compared to other studied locations, as depicted in Figure 11, Figure 12 and Figure 13. Consequently, meticulous technical planning and unique foundation designs are imperative for infrastructure development in this area. Identifying suitable foundation locations relies heavily on the competent (fresh) granitic unit, which varies significantly in depth across the study area, particularly beneath Sites 1 and 3. The RMQ (Figure 14 and Figure 15) and N-value models (Figure 16; Table 6) have greatly facilitated the assessment of foundation rock suitability. The integral/fresh granitic bedrock, characterized by RQD values of 75–100% (Site 2), and sections delineated with N-values >50 at Sites 1 and 3, are capable of supporting infrastructure foundations. For high-rise structures or buildings with continuous footings, reinforced foundations are recommended, preferably piled to rest on competent/fresh bedrock with resistivity, Vp, and RQD values exceeding 2000 Ωm, 2400 m/s, and 90%, respectively. This is particularly crucial for the eastern and northern parts, where penetrable competent bedrock is scarce at considerable depths. Pile foundations reaching the competent/fresh bedrock may extend to depths ranging from 12 m (less common) in the southern part (Site 2) to 40 m at Site 1 (except TR3) and Site 3 (Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15). Typical foundation depths for superstructures range from 20 to 40 m across the study area.

5.4. Implications of Modeled Soil-Rock Conditions on Groundwater Productivity

Understanding the surface–subsurface soil-rock conditions, as well as the structures facilitating fluid migration in exhumed bedrock, provides a comprehensive insight into the spatial distribution, channeling, and storage of groundwater (Gao et al., 2018; Kumar et al., 2021). This knowledge is essential for sustainable groundwater abstraction in complex geological terrains (Soupios et al., 2007; Dahlin and Wisén, 2018). Adequate overburden thickness (with minimal clay/silt content), well-connected fractures, deep-weathered troughs, and other factors contribute to enhancing the hydrodynamics of aquifer units, thereby increasing groundwater yield capacity. Therefore, integrating the Vp, resistivity, RMQ, and N-value models with their regression-derived counterparts is crucial for identifying productive water-bearing aquifer zones within the soil-rock profiles of the study area.

The 2-D/3-D Vp and resistivity models for the investigated sites (Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13) reveal intriguing surficial and subsurface soil-rock conditions and architecture, exhibiting contrasting geophysical characteristics. Notably, penetrative weathered/fractured zones serve as fluid migration systems within the study area. The boreholes in the region have water table limits ranging from 1.5 m to 5 m. By carefully analyzing Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13, the soil-rock conditions of the study area were assessed to pinpoint nearly silt-free sections, which are likely to offer optimal groundwater yield in boreholes.

Interestingly, Site 1, particularly TR1, exhibits high groundwater potential zones within its subsurface conditions. The topsoil thickness is less than 2 m, while the weathered unit exceeds 10 m, reaching over 41 m beneath TR1 (Figure 7 and Figure 8). Potential groundwater conduits include deeply weathered/fractured zones identified at station positions 85–125 m (TR1; Figure 7a) and 50–70 m (TR2; Figure 7b). Beneath TR3 (Figure 8), the western overburden soil constituent of USM is shallower than its eastern part. Here, weathered/fractured zones at station positions 40–50 m and 68–77 m may be viable drilling targets. However, due to the limited mapping of weak zone continuities, the certainty of high groundwater yield at these locations may be uncertain. At Site 2, residual soils and deeply weathered/fractured units exceed 35 m in thickness (Figure 9 and Figure 10). Utilizing SRT, ERT, and RMQ models (Figure 9, Figure 10a–b, and Figure 15), delineated multiple penetrative fractures may provide the necessary groundwater yield, particularly in the central to eastern parts. Site 3 shows a network of groundwater conduits and storage zones, with an overburden thickness exceeding 20 m (Figure 11, Figure 12 and Figure 13). Overall, these identified sections in the study area typically serve as sustainable water-bearing zones, with targeted drill depths above 40 m. However, due to significant silt content, especially at Sites 1 and 3, appropriate borehole completion procedures are proposed to ensure the long-term productivity of aquifers in the study area.

6. Conclusions

A holistic understanding of the spatial variability of surface–subsurface soil-rock conditions and architecture at strategic sites on Penang Island, Malaysia, is provided via the combination of multidimensional velocity–resistivity, RMQ, and SPT N-values with their optimized regression modeling. This study aims to determine the suitability of defined soil-rock conditions for the development of sustainable infrastructure and groundwater resources. Moreover, novel empirical relationships have been established to directly predict RQD and N-values from Vp and resistivity data, offering significant cost reductions for borehole drilling and associated research in granitic terrains on a large scale.

The Vp, resistivity, RMQ, and N-value models of the study area show distinct characteristics across the sites. Site 1 in the eastern region and Site 3 in the northern part exhibit thick, saturated, and loose silty to sandy materials, along with deep-weathered/fractured zones. Batu Maung (Site 2) primarily comprises sand, boulders, and resistive bedrock at shallow depths, accompanied by penetrative fractures. The N-value results align well with the resistivity and Vp models, indicating N-values of 7–17 for residual soils with a thickness extending to approximately 2 m. Moreover, N-values of 18 to 50 and >50 are estimated for the highly weathered to relatively weathered granites and the integral granitic bedrock, respectively.

Understanding the soil-rock characteristics, as well as the geometries and apertures of the weathered/fractured zones, provides valuable insights into selecting suitable sections for infrastructure and groundwater development in the area. It is important to address potential risks, such as soil failure in the highly steep soil segment at Site 2, through thorough soil excavation and reinforced retaining walls before infrastructure design. Generally, suitable sections for infrastructure placement are identified based on competent bedrock, characterized by RQD values exceeding 90% (Site 1) and N-values above 50 (Sites 1 and 3). Piling depths can vary, ranging from 12–25 m at TR3 (Site 1) and Site 2 to 40 m at Sites 1 and 3. The eastern to northern sections of the study area exhibit greater soil water-retention affinity compared to the south, attributed to deep-weathered/fractured zones with depths exceeding 40 m, which represent potential zones for sufficient water exploitation. Moreover, the intra-bedrock fractures beneath Site 2 have the potential to enhance groundwater productivity beyond other sections with high silt content.

This study effectively highlights the importance of integrating geophysical and geotechnical parameters for soil-rock characterization through optimized statistical modeling. The derived lithology-based empirical relationships for RQD and N-values with Vp and resistivity datasets offer valuable insights for geotechnical and civil engineers, as well as hydrogeologists. These empirical relationships enable direct predictions of RMQ and N-values, particularly in areas with limited or no borehole data, thereby facilitating informed decision-making in infrastructure and groundwater development projects.