Multivariate interrelatedness of geotechnical and petrophysical properties towards developing near-surface lithology clusters in a sedimentary terrain

The study reports the approach of interrelating geotechnical and petrophysical parameters to develop near-surface lithology clusters in a sedimentary terrain. Field-based and laboratory-based geophysical, geotechnical, and petrophysical measurements were carried out with a view to finding correlations among related parameters, grouping the parameters into depth clusters, and determining trend correlations between field-acquired and laboratory-obtained parameters. Four agglomerate hierarchical cluster signatures were identified at depth 0–30 cm, three at 30–60 cm, and five at 0–90 cm. Principal component 1 retained 74% of data variation and differentiated the petrophysical parameters according to Poisson’s ratio (v) and velocity ratio (Vp/Vs). Principal component 2 explained another 14% of variability in the original responses and separated parameters based on elastic modulus (E) and clay fraction (C). Degree of correlation of Vp/Vs with porosity (Φ), sand fraction (S), C, unconfined compressive strength (U), and California bearing ratio (R) followed order Φ > S > C, and U > R. Trend correlation of shear velocity at 1 m showed good agreement (R2 = 0.992) between field and laboratory data with an amplification factor of 4.13 separating the corresponding values. Lithology determined from soil analysis agreed with that inferred from measured electrical resistivity while geotechnical, petrophysical, and geophysical spatial analysis supported site lithology configuration patterns. Thus, multivariate association of geotechnical and petrophysical parameters exhibited potentials for rapid and exploratory near-surface lithology mapping even at a regional scale.


Introduction
Incessant failures of buildings and roads have been reported to cause colossal loss of lives and properties around the world, hence increasingly attracting the attention of investigators from related fields. This unpleasant situation is inimical to achieving sustainable cities and communities as it adversely affects the economy, environment, and safety of the citizens. Building collapse is a common phenomenon all over the world, but the trend in Africa is fast attaining worrisome dimensions. In Nigeria particularly, frequency of building failures and collapses has recently escalated with its associated colossal losses of lives and properties such that there were fifty-four cases of collapsed buildings across Nigeria between 2012 and 2016 (Ayininuola and Olalusi 2004). Out of the reported cases of collapsed building and failed structure in Nigeria, Lagos had the largest share of 51.6%. Statistical data showed that out of the reported cases of collapsed building and failed structure in Nigeria, Lagos, had the largest share of 51.6%, and out of the 105 cases of building collapses reported between 1978 and 2007 in Lagos, Island and Mainland recorded 30% and 11%, respectively (Oni 2010). The most recent one was a total collapse of a 21-story building that was still under construction on the Island on 1st November, 2021 where many were buried under the rubbles.
Factors reported to have been associated with failures of building and other engineering structures mostly included natural occurrences such as earth tremor, earth quake, land slide, erosion and flood; as well as technical and engineering defects such as overload, improper designs, and corrosion of reinforcement (Oni 2010).
Clustering and display instruments check for similarities between multiple objects and unrelatedness in a dataset thereby making correlations between variables (Beebe et al. 1998). Principal component analysis (PCA) explains variability in a dataset (Liu et al. 2018) while agglomerative hierarchical clustering (AHC) groups samples and interrelates groups to form patterns (Lee and Yang 2009). While PCA is usually deployed to identify underlying variables, AHC on the other hand helps to identify relatively homogeneous variables using selected characteristics (Hussain et al. 2008;Massart and Kaufman 1983;Schmitz et al. 2004). Soil stability investigations have always required both field-based and laboratorybased methods-geotechnical, geophysical, or geochemical approach-or an integration of either two or three of the approaches (Adebisi et al. 2016;Adeogun et al. 2020;Adeoti et al., 2016;Ayolabi et al. 2010;Barker 1997;Imai et al. 1976;Ishola et al. 2014;Oyedele et al. 2015;Park et al. 1999;Soupios et al. 2007;Stumpel et al. 1964), but these have always been very expensive especially if conducted at a regional scale. Attempts have also been made by previous investigators to characterize some selected rocks and soils geotechnically and petrophysically (Abatan et al. 2016;Boadu 2000;Boadu and Owusi-Nimo 2010;Kuforiji et al. 2021) for the purpose of determining their identities, and to develop correlations between some of the identities (Abd El-Rahman 1989; Bowles 2012; Moayed 2012; Roy and Bhalla 2017; Schmitz et al. 2004;Rezaei et al. 2018;Tribedi 2013). However, to the best of our knowledge, addressing ground stability by interrelating laboratory-obtained geotechnical and petrophysical identities towards developing unique field signatures at a regional scale is yet to attract significant attention. In addition, previous interests on laboratory-based seismic velocities have mostly focused on rocks and other solidframe geo-materials, owing to relative ease of measurement. There is therefore the need to establish rapid and explorative technique of assessing near-surface lithology that will be reasonably indicative of ground structure at a regional scale. Hence, this research was aimed at (i) obtaining relevant geotechnical and petrophysical parameters in the first 90 cm depth in the study area, (ii) interrelating obtained identities to develop near-surface lithology clusters using multivariate tools, (iii) spatially analyzing the study area using measured parameters and validating with acquired geophysical parameters, and (iv) correlating relevant field-based and laboratory-based parameters.

Study area
Lagos, which happens to be the second largest city in Africa, is located at 6° 34′ 60″ N, 3° 19′ 59″ E along the West African coast. The metropolitan city covers the substantial part of Lagos State and comprises of lagoons, sandbars, and islands connected to each other and to the adjacent mainland by bridges. The land mass is estimated to be of about 3500 km 2 with an estimated population of about 15.5 million people. The average annual rainfall is about 1629.4 mm with peaks in June, and between September and October. Mean minimum temperature ranges between 22 and 26.6 °C, highest in March and lowest in August, while the mean maximum temperature ranges between 27 and 34 °C. The surface geology of Lagos metropolis is made up of thick bodies of coastal plain sands known as the Benin formation (Miocene to Recent), and the recent alluvial deposit consisting of unconsolidated beach type sands with varying percentages of clay and mud (Fig. 1). On the other hand, the sub-surface comprises of varying lithological patterns of clay and sand deposits, and sedimentary rocks (Jones and Hockey 1964;Longe et al. 1987).

Methods and applications
This section contains explanations on both laboratory and field data acquisition methodologies, as well as the details of statistical analyses. Petrophysical study is expected to provide insights into the seismic and mechanical properties of the earth materials while geotechnical survey is expected to explain mechanics of soil and rock properties and with direct applications to design and construction of engineering structures. Geophysical investigations, on the other hand, have been designed to reveal subsurface lithologies as well as underlying physical processes with a view to providing guidance on placement of engineering structures.

Sample collection
Bulk samples were collected at depths 30 cm, 60 cm, and 90 cm with the aid of cylindrical core samplers (12 cm by 8 cm) and sealed with a polythene bag for onward transportation to the laboratories. Samples were cored in triplicate around same place and especially in the vicinity of buildings that have been pre-identified as distressed by the relevant agencies and pre-confirmed during reconnaissance survey, and roads that have failed. Incidentally most of the tarred roads adjoining such building structures were also observed to have failed. Locations where samples were collected are shown in Fig. 2 but it was not possible to achieve bulk samples at some particular depths for some locations.

Petrophysical measurements
CT 130 Pundit 200 UPV instrument containing two p-wave 54 kHz transducers and two s-wave 40 kHz transducers was deployed to measure the compressive and shear wave velocities (V p and V s ) of the samples (Fig. 3). Before each pair of transducers was affixed to the end of a core specimen, a zero adjustment reading was performed in order to account for the several microseconds between the time that the electrical pulse is received and then converted into a seismic pulse. To take a zero adjustment reading, a calibration rod was placed on a regular basis in between each pair of transducers when the transducer frequency is changed or the cables are changed with the expected calibration value (μs) marked on the calibration rod. The transducers were coupled to the calibration rod or the core sample by applying couplant to the transducers and both  (Adepelumi and Fayemi 2012) ends of the rod or sample and pressing firmly together with constant 50-MPa pressure applied across the core samples. Apart from the velocities, other parameters of concern determined were Poisson's ratio (deformation in direction perpendicular to the applied force i.e. ratio of the transversal elongation to the axial compression), and the elastic modulus (measure of the ability of a material to withstand changes in length under compression). Velocities of compressional and shear waves propagating through a rock sample were obtained when the sample width, d, is supplied to the ultrasonic machine using the respective transducers. The machine measures the pulse time, t, and automatically determines the wave velocity (v p and v s ). Using the measured densities, Poisson's ratio (v) was determined using Eq. (1). The ultrasonic generator is capable of automatically generating the Poisson's ratio with density and velocity as input data (Khalil and Hanafy 2016;Proceq 2011).
With v known as well as the velocity ratio (V p /V s ), elastic modulus (E) was then calculated from Eq. (2) (Adams 1951). (1)

Geophysical measurements
Geophysical investigations involving electrical resistivity and seismic refraction surveys were carried out on traverses connecting previously geotechnically sampled points in the vicinities of buildings pre-identified to be exhibiting foundation failures. Resistivity profile lines were conducted with the aid of PASI 16GL-N resistivity meter using Wenner and Schlumberger arrays to obtain 2D and 1D models, respectively (Koefoe 1979). Vertical electrical sounding (VES) was achieved using Schlumberger configuration to assess the vertical trends of the geoelectric parameters on traverses at depths. Apparent resistivity was measured based on: Current was introduced to the ground through a pair of electrodes separated by a distance AB, and the resulting potential was measured through another pair of electrodes separated by a distance MN. Half current electrodes separation, AB/2, was varied from 1 to 100 m at inter-VES station spacing ranging between 10 and 70 m. 2D resistivity survey was achieved using Wenner configuration to delineate lateral and vertical changes in the apparent resistivity values of the study location subsurface. Inter-electrode spacing was varied from 5 to 60 m in steps of either 5 or 10 m depending on availability of spaces (Telford et al. 1990). Inverted models of the subsurface were obtained from apparent resistivity data using iterative smoothness-constrained least-square inversion. Seismic profiling was done with a 24-channel ABEM Terralock Mark 6 seismogram using 5 kg sledge hammer as energy source based on the split-spread configuration. The refraction survey was carried out using 4.5-Hz geophones on 30 traverses to develop s-wave velocity models (Park et al. 1999). Acoustic energy was generated by vertically striking the metal plate with hammer shots. Geophone spacing of 2.5, 3, 3.5, and 4 m were used along the profile with a minimum offset of 1 m was employed for the investigation (Jorshi and Bhardwaj 2018). The longest source-to-the-last-geophone interval was 94 m on a maximum of two spread per traverse. VES curves plotted based on field data were processed using partial curve matching approach in order to obtain an initial model as input into an inversion procedure using WinResist 1.0 numerical package. RMS error on smoothed data from VES ranged from 1.1 to 2.9, with an average of 2.0. Word package was then used to process the 1D subsurface lithologies for all traverses. In addition, acquired 2D resistivity dataset were inverted to obtain a 2D image of the subsurface using the finite difference modeling (FDM) suite of the DiProWin4.0 numerical code which solves the forward problem of electrical resistivity for the purpose of generating theoretical dataset models. An integral step in the inversion procedure was to compare theoretical dataset to the observed resistivity data with the aim of achieving a good fit between the two, using a smoothness-constrained least-squares inversion algorithm and an active constraint balancing (ACB) (3) a = KV ∕ I to achieve stable results. The inversion process would normally be accepted to have achieved its aim, and hence consequently terminated, when the difference between the theoretical and observed resistivity models have been so iteratively minimized with an output of less than 5% mismatch, thereby resulting in the inverted subsurface 2D resistivity structure. In a similar vein, using SeisImager 4.2.0.0 suite applications, dispersion curves were calculated by converting to frequency domain and subsequently inverted based on the models initially generated. Inverted curves were then used to develop 2D s-wave velocity models and consequently the velocity profiles.

Scanning electron microscopic imaging and X-ray diffraction tests
The type and quantity of clay in any soil has much implication on the foundation design structure. Swelling clay minerals expand when moist and shrink when dry, leaving large voids in the soil. Samples obtained at 90 cm depth with high (IFK3) and low (LMA3) porosities were further processed to extract the clay content for scanning electron microscopy and x-ray diffraction mineralogy. SEM was used to obtain photomicrographs of the identified samples at 20 and 1000 µm. Images of the microstructures were picked at 100 × and 2500 × magnitudes of display, respectively. The images of the microstructures were later analyzed using the combination of ImageJ and OriginPro 2021 graphing and analysis tools (Fredrich et al. 1995;Mikrajuddin and Khairirrijhal 2009;Ziel et al. 2008) to determine the surface porosity of the samples. Since mineral composition play some defining roles, samples were further prepared for X-ray diffraction analyses using standard procedures (Kittrick and Hope 1963). Extracted clays from samples were therefore subjected to XRD tests (DX-27 SSC 40 kV/30 mA Slit:1 deg&1d) to determine mineralogy (Chauhan and Chauhan 2014).

Results and discussion
Micrographs and diffractograms were obtained for select samples and interpreted. Experimental data were subjected to multivariate analysis and spatial correlations for development of interrelated clusters and proper classification.

Multivariate statistics
Correlation matrix was first explored to reveal the hidden relationship between each pair of parameters. PCA and AHC were subsequently deployed on the data using XLSTAT, 2020, and OriginPro, 2020, to analyze for geographical origin and evolution of geophysical parameters. Mechanical properties (Table 1) were subjected to investigation on the basis of depths and locations (Granato et al. 2018). Column 1 represents sampling IDs, where suffices 1, 2, and 3 indicate sampling points at depths 30 cm, 60 cm, and 90 cm, respectively. It was possible to draw inferences on the mineralogical compositions especially for samples having complete information on compressive and plastic parameters (Casagrande 1932;White 1942) as listed on Table 2. From the analysis of results in Table 2 in combination with Table 1: SG gave an indication that samples are a combination of organic sands, ferrous sands, inorganic clay, and silty sands with values ranging from 1.89 to 3.89; R indicated a largely sandy clay structure with values ranging from 1.34 to 24.42%; SS and U showed medium to minimally stiff structure with values ranging from 42.41 to 72.6 kN m −2 , and 81.82 to 158.19 kN m −2 , respectively; LL indicated a soil structure low/medium plasticity of silt clay structure with range 17.85 to 55.82%; and PI showed a structure of 7% as sand, 23% as clay (C), 39% as clay loam, and 31% as silt clay (SC), with values ranging from 1.24 to 49.35%. Hence, the samples were either clay, silt clay or sandy clay thus supporting the instability indication for most of the sites. OC and OM indicated medium presence of clay from 1.88 to 5.87%, and 3.25 to 10.1% range of values respectively. Organic content is the ratio, expressed as a percentage, of the mass of organic matter in a given mass of soil to the mass of the dry soil solids. Organic matter influences many of the physical, chemical and biological properties of soils such as soil structure, soil compressibility and shear strength, water holding capacity, nutrient contributions, biological activity, and water and air infiltration rates. Equivalent basal spacing (EBS) was determined from LL using Eq. (4) (White 1942) for the purpose of inferring clay mineral compositions.
Though clay is very plastic and cohesive when moist, its extreme behavior of expansion and contraction, subject to moisture content or otherwise, puts pressure on foundation thus inhibiting its stability. Sand/gravel drains very easily because of the large openings and large particle sizes, and when compacted and moist, it holds together fairly well, thus providing good stability support though vulnerable to loss of friction and cohesiveness with heavy water movement. Since a PI above 25% gives an indication of expansive/shrinking clay when moist/dry, and given the EBS values ranging from 5.92 to 6.86 Å which will ordinarily infer kaolinite mineral, hence the studied clay samples can be said to be mixtures of minerals with inferred mineral composition (IMC) as kaolinite, illite, smectite, and quartz in different proportions.

Correlation and regression analysis
As a way of providing preliminary understanding of the processes controlling stability parameters, Pearson's correlation coefficients were calculated using XLSTAT, 2020, as shown in Table 3. Pearson's correlation coefficient is a statistical test that measures relationship, or association, between two continuous variables on the same interval based on the method of covariance. Coefficient with value −1, 0, or 1 is an indication of a total negative linear correlation, no correlation, or a total positive correlation, respectively. For samples collected at 30 cm depth, perfect positive correlation was expectedly observed between U and SS, and ρ and SG. U was equally observed to correlate positively with R, R being a penetration test which evaluates mechanical strength, and U being a test which evaluates the strength of the soil when it is crushed uniaxially without lateral constraint. For samples collected at 60 cm depth (Table 4), similar trends were observed between U and SS, and between ρ and SG. However, unlike at 30 cm depths, positive correlation was between U and K. Poisson's ratio is usually positive given that most common materials become thinner in cross section when stretched. A zero Poisson's ratio is an indication that there is no transverse deformation resulting from an axial strain. The strength and stability of soil depend on its physical properties, and it is important that soil is stable through wetting and drying cycles to avoid cracks on roads and foundations. Expansive soil poses significant threats to foundations through its capacity to swell thereby causing uplift with moisture increase and resulting in cracks. Hence, the investigated soils did not have the capacity to provide needed support for foundations because soil that demonstrates stability across wetting and drying cycles is most appropriate for structural foundations (Firoozi et al. 2017;Mitchell and Soga 2005). As a matter of fact, an appropriate combination of particle types is best for structural stability. It was important to investigate how properties correlate across depths. Correlation of parameters in Table 5 was carried out, exclusive of mechanical properties due to incompleteness of data. Correlation of the velocity ratio with Poisson's ratio, and elastic modulus followed expectation since Poisson's ratio was calculated with known values of V p and V s , and elastic modulus was calculated with known values of bulk density and Poisson's ratio. Velocity ratio correlated significantly, though negatively, with density. Correlation of velocity ratio with porosity, sand fraction, and clay fraction was according to the order Φ > S > C (Table 6). Hence, regression equation was developed for the velocity ratio and the three parameters as expressed in Eq. (5).
Correlation of available mechanical properties with corresponding velocity ratio was investigated. Table 7 presents the correlation matrix between velocity ratio, California bearing ratio, and unconfined compressive strength. V p /V s correlated positively with U and R according to the order U > R. Hence, regression equation was developed for the velocity ratio and U as expressed in Eq. (6).  LGA

Principal component analysis
According to the Kaiser criterion, factor loading (FL) of at least 0.4 is to be considered important to extract sufficient variance from the variable (Nagaraju et al. 2016). Also, eigenvalues higher than 1 are considered as significant in the PCA analysis and are normally used for rotation which is also used to explain the scree plot. FL is the correlation coefficient for the variable and the factor, which reveals the variance explained by the variable on that particular factor. Hence, petrophysical parameters (Table 5) were subjected to PCA across the entire depths (0-90 cm). PC1 retained about 74% of data variation and differentiated the petrophysical parameters according to Poisson's ratio and velocity ratio (Table 8). Similarly, PC2 explained another 14% of variability in the original responses and separated the parameters based on elastic modulus and clay fraction. The signs indicated that V p /V s will be lower in the negative axis of the PC1 while v will be higher in the positive axis of PC1, and both E and C will be higher in the positive axis of PC2 (Fig. 4). Table 9 explains the score labels for PC1 and PC2 where scores represent the coordinates of the score labels on the biplot, while the score labels indicate the sampling IDs. For PC1, samples located on the left-hand side of the biplot (IKJ1 and IFK3) have lower mean values of velocity ratio than those located on the right-hand side, while samples located on the right-hand side of the biplot (IKJ2a, IKJ2b, LMA2, LMA3, IKT2, LIS1, and IKO2) have higher values of Poisson ratio than those on the left-hand side. For PC2, samples located on the upper side of the biplot (IKJ1, IKJ2b, LMA3, and IKT2) have higher values of elastic modulus and clay fraction than those on the lower side.

Hierarchical cluster analysis
AHC was used to group sample locations in clusters based on their similarity. AHC makes grouping decisions by considering local patterns or neighbor point rather than global data distribution. Group average was the cluster method used to build the cluster stages while the distance type was Euclidean. Each sample location was placed in its own cluster and then subsequently merged into larger clusters until a single cluster was achieved. Sum of distances method was used to determine the most and least representative observation. The final outcome was in form of a dendogram where the vertical axis represents the distance or dissimilarity between clusters while the horizontal axis represents the clustering sample locations. Hence, geotechnical parameters (Table 1) were subjected to AHC at two levels of depths (0-30; 30-60 cm), while petrophysical parameters (Table 5) were similarly subjected but across the entire depths (0-90 cm). Figure 5 illustrates sampling locations at 30 cm with unique properties associated to clusters which can be easily used as stability signatures. Sampling points can be grouped into four clusters: {KSF, BDG} belonging to cluster 1; {IKJ} is clearly as an outlier sampling point which can be classified as cluster 2; {IFK, APA, LMA} categorized as cluster 3; and {LIS, EPE} classified as cluster 4. The dissimilarity is in the order of cluster 3 < cluster1 < cluster 4 < cluster 2. Hence, clusters 1, 2, 3, and 4 are said to have unique geotechnical fingerprints at 30 cm depth. For samples locations at 60 cm as shown in Fig. 6, three clusters appeared very distinct: {EPE, IKJ} belonging to cluster 5; {ISO} was an outlier cluster 6; and {LMA, APA, BDG} categorized as cluster 6.The dissimilarity was in the order of cluster 7 < cluster 5 < cluster 6. Hence, clusters 5, 6, and 7 are said to have unique geotechnical fingerprints at 60 cm depth. Figure 7 presents petrophysical properties across depths with EPE1 as the clustroid, i.e., the sample location closest to other sample location in terms of geotechnical and petrophysical properties. Samples {BDG3, IFK1, IKJ2} fall into cluster 8 of sample locations with the least dissimilarity, followed by {EPE1, IKT2} belonging to cluster 9, {IFK2, IKT3} categorized as cluster 10, {LMA3} classified as cluster 11, and finally {IFK3} as clearly an outlier sampling location, with the widest dissimilarity compared to the rest points this presenting a markedly different fingerprint, and classified as cluster 12. Hence, in the order of dissimilarity, cluster 8 < cluster 9 < cluster 10 < cluster 11 < cluster 12, and are said to have unique petrophysical fingerprints across the depths.

Morphological and mineralogical analysis
Visual interpretation of the micrographs at 1000 µm revealed clay minerals ranging from moderate porosity (LMA3) to low porosity (IFK3) which was consistent with the surface porosity trends (Table 10) theoretically determined from OriginPro. LMA3 (Fig. 8) and IFK3 (Fig. 9) have surface porosity of clay extracts separated by a wide margin of about 10% (Mikrajuddin and Khairirrijhal 2009). This further confirmed that clay types for the two locations have unique structural and morphological identities (Fredrich et al. 1995). X-ray diffraction (1.54056A, Cu/Kalpha1) run on select soil-clays revealed alluminosilicates as the prominent mineral (Kittrick and Hope 1963). Presence of silicates renders the soil very unstable because of the swell-shrink characteristics during wetting and drying cycles, and hence puts engineering loads at risk of failure.

Site classification using laboratory-obtained shear wave velocity
In order to classify the capacity of the soils in the 1 m depth, shear wave velocity was used to determine stability indicator parameters such as standard penetration resistance or N-value (N), bearing capacity (B r ), and reaction module (R m ) on samples with known V s . The parameters were calculated using the following expressions (Imai et al. 1976;Moayed 2012;Stumpel et al. 1964;Tribedi 2013):  N-value gives an expression on soil resistance to penetration or cohesiveness (Eq. (7)); B r is an indication of the maximum load required to produce soil shear failure (Eq. (8)), while R m is an efficient indicator for structural analysis of foundation stability measured as load intensity per unit of displacement (Eq. (9)), also known as modulus of subgrade reaction, and it expresses the relationship between soil pressure and deflection (Abd El-Rahman 1989;Birch 1966;Gassman 1973;Sheriff and Geldart 1986;Tatham 1982

Correlation between field-acquired and laboratory-obtained shear wave velocity
Field-based topsoil (V s 1) shear wave velocity was extracted for Ikorodu, Ebute Meta, Ikeja, Lagos Island, and Ifako Ijaye (Table 12) for the purpose of comparing with corresponding laboratory-based values. In similar vein, topsoil electrical resistivity was extracted for the same locations with corresponding inferred lithology for the purpose of comparison and validation. Trend correlation of shear velocity showed good agreement (R 2 = 0.992) between laboratory-based and field-based data though an amplification factor of 4.13 separates the corresponding values. Lithology determined from soil analysis showed good agreement with that inferred from the electrical resistivity, and trend correlation of shear velocity showed good agreement (R 2 = 0.992) between field and laboratory (Fig. 10).

Spatial analysis of laboratory-based petrophysical and field-based geophysical parameters
The ArcMap and the ArcScene Software of the ArcGIS 10.4.1 package were used for the interpolation of data in 2D, and the visualization of interpolated data in 3D respectively.
Soil properties data for 30 cm, 60 cm, 90 cm depths which were in excel (.csv) data was tailored to suit the GIS environment. This was done by finding the mean for certain parameters and adding coordinates to the same. This data cleaning and engineering was done for the 30 cm, 60 cm, and 90 cm, respectively, while all parameters were scaled from −10 to 10. Geophysical data for 30 cm depth which was in Microsoft excel (.csv) data was imported into Arc-Map and was plotted using the X and Y coordinates on the attribute table. Plotted points were exported as a layer in shapefile (.shp) data format. This process was repeated for the geotechnical and petrophysical properties of the samples for the 60 cm and 90 cm depths, respectively. Inverse distance weighting (IDW) interpolation tool under the Arc-Map geostatistical tool was used for the interpolation of 30 cm, 60 cm, and 90 cm depths, respectively. Exported shapefiles served as the input for the analysis. Geotechnical and petrophysical parameters to be interpolated were selected respectively while Lagos boundary shapefile layer was used as the extent of the analysis which after analysis generated an output which formed a new raster layer for each of the parameter that was interpolated. The interpolated parameter layers were imported to ArcScene to render depths of each parameter for ease of interpretation and visualization.  Laboratory V s Field V s Figures 11 and 12 represented spatial spread of carbon contents-organic carbon and organic matter, and spatial distribution of lithological properties-silt fraction, sand fraction, clay fraction, and silt and clay fraction respectively. Figures 13 and 14 respectively show spatial distribution of geotechnical plastic properties-Atterberg limits and plasticity index, and spatial spread of geotechnical compressive properties-shear strength, California bearing ratio, and unconfined compressive strength. Figures 15 and 16 show spatial distribution of geophysical compaction properties-bulk density and maximum dry density, and spatial distribution of geophysical flow properties-porosity, permeability, and void ratio, respectively. Figures 17 and 18 represent spatial spread of geophysical moisture retention properties-natural moisture content and optimum moisture content; and spatial distribution of petrophysical properties-p-wave velocity, s-wave velocity, elastic modulus, and Poisson's ratio, respectively. Hence, clusters were further developed based on similarities of signatures.
Trend correlation using visual pattern recognition revealed clusters of parameters with similar spread patterns as follows: cluster I-sand fraction, porosity, bulk density, p-wave velocity, and s-wave velocity; cluster II-shear strength, unconfined compressive strength, and California bearing ratio; cluster III-elastic modulus and Poisson's ratio; cluster IV-clay fraction, liquid limit, plastic limit, and plasticity index; cluster V-permeability and void ratio; and cluster VI-silt fraction, organic carbon and organic matter. The six clusters revealed innate interrelatedness among the parameters.
For the purpose of validation, spatial maps for fieldacquired parameters (shear wave velocity and electrical resistivity) were developed for traverses covering specific points where geotechnical analysed samples have been obtained. From Ikorodu traverse 1, s-wave velocity profile revealed three poorly competent clay layers on top of a competent layer of sand. 2D profile length of lateral length 180 m and depth of 50 m with resistivity values ranging from 14 to 263 Ωm revealed three incompetent layers sitting on top of a competent layer. Top layer was clayey sand with apparent resistivity values of 14-29 Ωm at 3 m depth along 20-170 m, second layer (sandy clay) was observed at 18 m depth with resistivity value ranging between 91 and 120 Ωm, third layer (sand) spreading    (Fig. 19). With respect to Ebute Meta traverse 2, s-wave velocity profile revealed three poorly competent layers (clay/sandy clay, clay, and peat) sitting on top of a moderately competent layer of sand. Electrical resistivity model showed first, second, third and fourth layers as clay/sandy clay (poorly competent), clay (poorly competent), peat (poorly competent), and very thin continuous peat (moderately competent) respectively. Geoelectric section combining VES stations 1 and 2 at 45 and  (clay), third layer (sandy clay), and fourth layer (sandy clay) across the VES stations. Resistivity ranged between 54 and 62, 10 and 14, 2, and 38 and 57 Ωm at depths 1, 11, 33, and > 33 m, respectively (Fig. 20). For Ikeja traverse 3, s-wave velocities revealed poorly competent to fairly competent, poorly competent, and moderately competent profiles. Inverted models showed top layer (clay/clayey sand), second (sandy clay), third (sand/clayey sand), and fourth (clayey sand) lithologic unit with model resistivity values of 17-20 Ωm to depth 2. Model resistivity values of 70-110 Ωm along the lateral distance 80-120 m at 17 m, 80-120 Ωm flowing across the profile beneath the second layer to a depth of 50 m, thus prone to structural failure because of overburden. First three layers of VES1 and VES 2 have resistivity ranging from 38 to 98, 29 to 60, and 99 to 169 Ωm at depths 0.7-1, 4-5, and 38-40 m, respectively. Clay was observed under sand layer of VES 2 between 40-50 m depths with resistivity value of 12 Ωm. Geoelectric section combining VES point 1 and 2 at 80 and 100 m respectively up to 50 m depth identified four subsurface layers comprising of incompetent top soil, incompetent second layer (sandy clay and clay), competent third (sand and clayey sand) and incompetent fourth layer (clay and sand) unevenly spread across the VES stations especially at the third and fourth subsurface layers with resistivity ranging between 38 and 98, 30 and 60, 99 and 178, and 13 and 155 Ωm at depths 1, 3, 38, and > 38 m, respectively (Fig. 21). With respect to Lagos Island traverse 1, s-wave velocity profile revealed three poorly competent clay/peat layers on top of a competent layer of sand. Electrical resistivity profile showed top layer sand/ clay lithologic unit with model resistivity values of 2-52 Ωm, along the lateral distance 0-180 m on the traverse, second layer was clay/peat with resistivity values ranging from 6 to 52 Ωm up to 30 m depth, and third was peat/ sandy clay with resistivity ranging from 10 to 62 Ωm up to 50 m. Geoelectric section combining VES stations 1 and 2 at 100 and 120 m respectively and covering a total depth of 40 m showed four incompetent subsurface layers comprising of top soil, second layer (clay) third layer (peat) and fourth layer (sandy clay), unevenly distributed across the VES stations. Resistivity ranged from 147-149, 9-18, 6-8, and 40-56 Ωm at depths 1, 3, 21, and > 21 m, respectively (Fig. 22). For Ifako Ijaye traverse 1, s-wave velocity profile indicated fairly competent, incompetent and competent lithologies. 2D ERT spreading over a lateral distance of 180 m and 50 m deep showed three layers comprising of sandy-clay, clay, and sand. Top layer lithologic unit revealed apparent resistivity values of  Fig. 18 Spatial map of laboratoryobtained V p , V s , E, and v 91-100 Ωm, second layer (sandy clay) was observed at distance 10-190 m to a 15 m depth with resistivity values ranging between 91 and 120 Ωm, and third layer was sand spreading across the profile with resistivity values of 101-534 Ωm up to 50 m depth. Geoelectric section identified uniform layers across the VES stations with three subsurface layers comprising of incompetent top soil (sandy clay), incompetent second (clay), and competent third (sand) with resistivity ranging from 40 to 47, 21 to 22, and 199 to 730 Ωm at depths 1, 3, and > 3 m, respectively (Fig. 23). Hence, field acquired s-wave velocity and electrical resistivity profiles revealed incompetent, poorly competent, fairly competent or moderately competent lithological classification of layers, and this is again consistent with the trends observed at the surface from laboratory-obtained geotechnical parameters.
However, the unavailability of s-wave velocity data beyond 20 m for most of the sampling traverses posed some limitations to accurate site classification; hence, s-wave velocity at 30 m depth was further required for proper classification of the studied sites. Vs30, the average seismic shear-wave velocity from the surface to 30 m depth, has been widely used by many investigators as a parameter for seismic site classification to characterize  site suitability for building codes.Vs12.5 and Vs20 were used to determine Vs30 using Eq. (10) (Boore 2004;Sheriff and Geldart 1986;Wang and Wang 2015).

Conclusions
In this work, multivariate statistical tools were used to interrelate measured geotechnical, petrophysical, and geophysical parameters towards developing stability signature on a sedimentary terrain. Correlations were developed for the geotechnical and petrophysical parameters, and for fieldacquired and laboratory-obtained geophysical parameters. Principal component analysis tool was used to explain the variance in the data while agglomerate hierarchical clustering tool was used to identify clusters according to dissimilarities and relatedness. Principal component 1 retained about 74% of data variation and differentiated the petrophysical parameters according to porosity and velocity ratio, while principal component 2 explained another 14% of variability in the original responses and separated the parameters based on elastic modulus and clay fraction. Four agglomerate hierarchical cluster lithologies were identified at depth 0-30 cm, three at 30-60 cm, and five at 0-90 cm. Lithology determined from soil analysis showed good agreement with that inferred from the electrical resistivity, and trend correlation of surface shear velocity showed good agreement (R 2 = 0.992) between field and laboratory. Spatial maps were developed for geotechnical, petrophysical, and geophysical parameters for the purpose of validating site configuration patterns. Field acquired s-wave velocity and electrical resistivity profiles revealed incompetent, poorly competent, fairly competent or moderately competent lithological classification of layers, and this is consistent with the trends observed at the surface from laboratoryobtained geotechnical parameters. Interpolated shear wave velocity at 30 m depth supported geophysically measured electrical resistivity profiles with largely incompetent site classification of 63% E (soft clay profile) and 33% D (stiff soil). Multivariate association of geotechnical and petrophysical parameters exhibited good potentials for rapid and explorative near-surface regional lithology mapping.

Availability of data and materials
Most of the data generated or analyzed during this study are included in this published article, while the remaining are available from the author on reasonable request.

Declarations
Ethics approval and consent to participate N/A.

Conflict of interest
The author declares no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.