Optimizing Subsurface Geotechnical Data Integration for Sustainable Building Infrastructure

Abstract

Sustainable building construction encounters challenges stemming from escalating expenses and time delays associated with geotechnical assessments. Developing and optimizing geotechnical soil maps (SMs) using existing data across heterogeneous geotechnical formations offer strategic and dynamic solutions. This strategic approach facilitates economical and prompt site evaluations, and offers preliminary ground models, enhancing efficient and sustainable building foundation design. In this framework, this paper aimed to develop SMs for the first time in the rapidly growing district of Gujrat using the optimal interpolation technique (OIT). The subsurface conditions were evaluated using the standard penetration test (SPT) N-values and soil classification including seismic wave velocity to account for seismic effects. Among the different geostatistical and geospatial models, the inverse distance weighting (IDW) model based on an optimized spatial analyst approach yielded the minimum error and a higher association with the field data for the understudy region. Overall, the optimized IDW technique yielded root mean square error (RMSE), mean absolute error (MAE), and correlation coefficient (CC) ranges between 0.57 and 0.98. Furthermore, analytical depth-dependent models were developed using SPT-N values to assess the bearing capacity, demonstrating the association of R² > 0.95. Moreover, the study area was divided into three geotechnical zones based on the average SPT-N values. Comprehensive validation of different strata evaluation based on the optimal IDW for the SPT-N and soil type-based SMs revealed that the RMSE and MAE ranged between 0.36–1.65 and 0.30–0.59, while the CC ranged between 0.93 and 0.98 at multiple depths. The allowable bearing capacity (ABC) for spread footings was determined by evaluating the shear, settlement, and seismic factors. The study offers insights into regional variations in geotechnical formations along with shallow foundation design guidelines for practitioners and researchers working with similar soil conditions.

1. Introduction

The safety and durability of any building structure are intrinsically linked to the geotechnical characteristics of the underlying soil [1,2]. In the absence of a thorough geotechnical investigation, the stability and load-bearing capacity of subsoil layers is indeterminate, presenting considerable risks to structural integrity, serviceability, and overall building performance, which may result in impaired functionality or structural failure [3]. However, such investigation procedures are cumbersome and require plenty of money and time, considerably escalating the project cost. Conversely, despite swift technological progress in the construction sector, urban subterranean regions remain predominantly unexplored and unutilized [4]. Additionally, in many urban and developing areas, plenty of projects have been executed in the past, and their geotechnical investigation reports (GIRs) ultimately became part of the documents rather than as guidelines for future projects. Such GIRs and subsequent structural design are based on real-time data, tedious laboratory and field testing, and decades of experience, which provide an in-depth understanding of soil behavior for a proposed project. Moreover, for projects of this magnitude, the preliminary feasibility evaluations are typically founded on fragmented information from multiple geotechnical reports instead of systematic data [2]. Conversely, in numerous small-scale construction projects, particularly within the urban residential expansion, geotechnical investigations are frequently bypassed due to economic constraints and the absence of comprehensive foundation design guidelines tailored to site-specific spatial characteristics [5]. Consequently, incidents of structural failure in small-scale buildings are prevalent in developing South Asian countries, often due to inaccurate estimations of the foundation design parameters [6]. To overcome these problems, the development of optimized SMs that account for varying geological formations and geotechnical settings can provide critical insights for site-specific foundation design. Such maps provide a more systematic approach by integrating data from earlier geotechnical studies and customizing it for future usage, improving the dependability of both small and large-scale building projects. This study aimed to create optimized geotechnical soil maps (SMs) for the rapidly expanding district of Gujrat as well as a geotechnical zonation map and design recommendations for shallow foundations.

In the past few years, the district of Gujrat has undergone a rapid development of infrastructural work as the population has been increasing at a tremendous rate. This district is pivotal in terms of agriculture and industrial trade; moreover, this area is planned to be connected to the China Pakistan Economic Corridor (CPEC) in the near future [7]. To meet the demand of the growing population, development works have been expedited, and various housing projects have been initiated by the Government of Pakistan under public–private partnerships in recent years. Present and future development work in the regions have given attention to the significance of SMs as these maps aid designers in delineating soil strength and qualities, offering essential pre-assessment and appraisal of the allowable bearing capacity (ABC). Moreover, they can also provide guidelines for local and international practitioners as well as researchers dealing with similar soil conditions for various infrastructure development work, where the project cost is insufficient to independently set up a geotechnical investigation.

Many researchers around the world have worked on soil mapping [8,9,10,11]. Nevertheless, there has been a scarcity of studies undertaken in the past that concentrated on several regions of Pakistan as research areas. The research conducted on the spatial interpolation of the district of Multan, Pakistan by using a large set of geotechnical data led to the development of shallow foundation design guidelines for the subjected area [12]. Likewise, the geotechnical zonification maps of the district of Faisalabad split the area into three zones based on the SPT-N ranges of soil (i.e., soft, medium, and stiff consistency soil) [13]. These maps provide reliable information on the geometry and engineering properties of underground layers that would make projects safer and more economical, thus helping to reduce the cost of geotechnical investigations, or at the very least, to have an initial understanding of the soil properties at the proposed project site. Consequently, the authors underscore the need to build such SMs, especially in quickly evolving urban areas or regions of considerable agricultural and economic significance such as the district of Gujrat. It is important to mention that no previous studies have been undertaken in this region to formulate foundation design recommendations for researchers and practitioners in the district of Gujrat in Pakistan. Such investigations are essential as they allow local and international practitioners to use data from adjacent or geologically similar places to create preliminary ground models. Accurate data on the subsurface geometry and qualities will improve safety, lower the project expenses, and offer a preliminary comprehension of the soil characteristics for new construction initiatives.

Numerous researchers have utilized diverse methodologies for the advancement of SMs, with geographic information systems (GISs) emerging as a formidable instrument for the collection, storage, retrieval, transformation, and visualization of spatial data from the physical space [14,15]. Studies show three diverse but overlapping views regarding the database, spatial analysis, and maps. Recent research has investigated the application of GISs in geotechnical engineering, employing spatial interpolation techniques to analyze geotechnical data [16,17]. GIS software encompasses various techniques of interpolation, but perspectives of its selection vary among different researchers. Previous research revealed that the ordinary kriging is less accurate in comparison to the inverse distance weighting (IDW) technique since it requires larger numbers of uniform spatial distributed data, which are rarely met in geotechnical data interpolation. Moreover, the efficacy of the IDW technique becomes evident in situations where the spatial distribution of data cannot be effectively modeled using variogram functions, as required in the ordinary kriging technique. Furthermore, GIS interpolated zonation maps of Surfers Paradise in Australia were explored via the IDW technique, revealing the varying spatial ranges of the geotechnical database without consideration of the efficiency of the prediction at unknown locations. Additionally, other researchers are of the view that the IDW method is relatively suitable for computing the results because of its straightforward interpretation as well as adjustment of the diminishing strength of the relationship over the study area via constant power or a distance decay parameter [17,18]. Many researchers have also developed geological, geophysical, and geotechnical maps using GIS-based coding and the analysis of soil investigation data [19,20]. Such maps are essential for guiding design principles, construction practices, and building regulations. Providing accessible data on the subsoil characteristics reduces the soil exploration costs and offers solutions to anticipated geotechnical challenges. Additionally, the development of these maps contributes significantly to regional and global geotechnical data repositories.

With the rapid expansion in infrastructure development in Pakistan, there is a strong urge to develop SMs for its major cities and districts to promote sustainable construction and geotechnical planned urban expansion in the region. The current study is an effort to develop such multi-stratified SMs for the unexplored district of Gujrat based on a field and laboratory investigated database entailing crucial parameters, such as the soil compressive resistance against the loading and soil type, to geospatially identify and comprehend the feasible and non-feasible zones for the prospective planned urban expansion. Furthermore, the study incorporates an extensive assessment of various interpolations with an optimization approach to identify the optimal interpolation technique (OIT). In addition, cross-validation was carried out based on critical performance metrics (CPMs) to ensure a quantification assessment of the accuracy of the generated SMs to ensure on-field applications. Moreover, the design considerations for shallow foundations based on ABC under the settlement, shear, and seismic criteria were also analyzed, which act as a guideline for local and international practitioners and researchers working on lightly loaded structures.

2. Description and Methods

2.1. General Description of the Study Area

The district of Gujrat is the tenth-major district of Punjab, Pakistan, comprising an area of about 3192 km², with a highly diverse environment. The district of Gujrat is located between 32°35′ north latitude and 73°45′ east longitude. The district is positioned on average about 233 m above sea level surrounding the two main rivers, the Chenab and the Jhelum. The trunk road network has been extended all over the district, which will be linked with the CPEC to create a trade course for the area. Figure 1 presents the location of the district of Gujrat on a map of Pakistan and its administratively divided tehsils, along with the average climate of the study area. Table 1 presents the distribution of land cover of Gujrat in km² carried out via survey by ALTAS Pakistan [21].

The district of Gujrat is underlain by a thick alluvial soil stratum comprised of clays, silts, and sand accumulated by the Chenab and Jhelum Rivers, which currently stream in the direction of the northern and southern sides of the study area, respectively. As per the geological survey of Pakistan, the project area lies in sub-piedmont deposits of the Indus flood plain, which is predominantly comprised of fine-textured unconsolidated soil deposits as deep as 900 m [22]. As per the United States Geological Survey (USGS), the district of Gujrat contains Neogene and Quaternary sediment deposits. The Neogene sediments are predominantly characterized by weathered clayey facies alternated by sand layers [23]. On the other hand, the quaternary deposit encompasses the environmental deposits (i.e., alluvial, fluvial, aeolian, lacustrine, and peat) as per the soil genesis. The United Nations Educational, Scientific, and Cultural Organization (UNESCO), in partnership with the Food and Agriculture Organization (FAO), has created the World Soil Map, which classifies soils according to their chemical, physical, and biological properties [24]. According to the UNESCO/FAO world soil maps, the district of Gujrat is governed by eutric cambisols, lithosols, calcaric fluvisols, orthic luvisols, and haplic xerosols representing weak weathered rocks, hard weathered rocks, flood plain alluvial deposits, illuvial clay deposits, and semi-desert soil, respectively [25]. The nomenclature of different soil groups along with its distribution within the boundary of the district of Gujrat is presented in Figure 2b. As per the UNESCO/FAO, the validity of the soil map was limited to the top 50 cm below the existing ground level. It should be noted that the network of lined and unlined irrigation channels is the main source of groundwater recharge [26]. Figure 2a represents the geology of the district of Gujrat.

2.2. Description of the Seismicity

As per the revised Building Code of Pakistan (BCP) (2007), the study area falls in seismic zone-2B. The peak ground acceleration (PGA) of zone-2B has a range of 0.16 g to 0.24 g [27,28]. The study area in the Punjab Plain shows the medium-intensity stage of seismicity. It is therefore recommended that the project structures should be intended to serve the constraints of zone-2B after giving due consideration to the foundation soil material. The characteristics of earthquake motion depend upon shear wave velocity, vs (m/s), which is an important parameter to calculate the seismic allowable bearing capacity.

2.3. Description of the Database

Geotechnical investigation reports collected from various certified public–private design firms, which are regulated by engineering governing bodies, contain information regarding the geographical, geological, or geotechnical data in both numerical and alphanumerical forms. Table 2 elaborates on the detailed information retrieved from GIRs. In this study, the GIRs from 130 different construction projects that fell within the jurisdiction of the district of Gujrat were extracted, and the average information of the borehole was retrieved and compiled for further analysis. The borehole coordinates were marked based on the information acquired from the GIRs, as seen in Figure 3. The subsoil data obtained from each borehole comprised the thickness and position of each stratum as well as the SPT-N at different depths. It is important to note that the SPT-N values used for the development of SMs were corrected for hammer energy, tendon length, water table, and overburden pressure as part of the GIRs. These corrections ensured the accuracy and reliability of the SPT-N data, effectively standardizing and normalizing the values to facilitate comparability across various datasets. Furthermore, major data from 103 sites were employed for the development of SMs, and the remaining data were designated for validation purposes, which were strategically distributed uniformly throughout the study region to ensure comprehensive coverage. The choice of the selection of 80% of the data for SM development, with the remaining 20% for validation purposes, was guided by the data availability and aligns with common practices to prioritize model development, ensuring the efficacy of the resulting SMs [29]. To further evaluate the data distribution, the haversine approach was used to analyze the minimum and maximum distance between the two closest points. The minimum and maximum distance between the two closest points between boreholes were found to be 1.1 km, and 5.5 km, respectively (Equation (1)). Furthermore, the model’s performance was evaluated and quantified via CPMs such as RMSE, MAE, and CC to confirm the suitability of the data distribution. This distribution was found to provide an optimal balance between model training and validation accuracy.

(1)$d=2 Rarcsin (\sqrt{si{n}^2\left({\displaystyle \frac{\Delta lat}{2}}\right)+\cos \left(la{t}_1\right)\cos (la{t}_2)si{n}^2\left({\displaystyle \frac{\Delta lon}{2}}\right))}$

Here, R = 6371 km is the Earth’s radius, while $\Delta lat$= $la{t}_2$-$la{t}_1$ and $\Delta lon$ = $lo{n}_1$-$lo{n}_2$, respectively.

Table 2

Data extracted from the geotechnical investigation reports.

Borehole ID	Contains general details about the borehole including identification number, project name, point, depth, location, contractor, and other relevant information.
Groundwater table	Records the water table depth and its variations over the monitoring period.
Lithology	Provides a detailed description of soil layers encountered including thickness, consistency, and color. For rock layers, additional information such as aperture, roughness, discontinuities, and weathering effects is included.
In situ tests	Contains data from various tests conducted within the borehole, offering reliable insights. Established correlations further aid in understanding the soil’s mechanical properties.
Lab tests	Includes test results from laboratory analyses on the soil and rock samples. Details on the sample quality, depth, sampling method, and physical and mechanical parameters are also provided.

Figure 3

Study area with the locations of the borehole points.

[Figure omitted. See PDF]

3. Results and Discussion

3.1. Statistical Analysis of the SPT-N

A statistical analysis of the SPT-N dataset was performed to evaluate its homogeneity, variability, and distribution trends, as presented in Table 3. Descriptive statistics including the mean and standard deviation (S.D) revealed a direct relationship with depth, indicating stiffer strata beyond 3 m below the existing ground level (EGL). This increased heterogeneity reflects transitions in geological strata and the presence of denser interbedded layers. Furthermore, thee analysis of variance (ANOVA) confirmed statistically significant differences in the SPT-N values across depth intervals, providing strong insights into depth-dependent soil behavior, as detailed in Table 4. It was evident from the results of the analysis that the values of variance and S.D increased down the depth due to the higher dispersion of SPT-N data. The majority of SPT-N values at shallow levels (i.e., 1–3 m) fell under the soft to medium consistency limits of soil, which depicted lower values of S.D and variance. Beyond a 3 m depth, soil from soft to stiff consistency ranges was observed, indicating that a larger gap among the lower and higher extremities of SPT-N values resulted in higher deviation.

The frequency distribution of the SPT-N values at various depths is plotted in Figure 4a–i. These graphs ensured that the dataset accurately represented the geotechnical characteristics, while the statistical analysis validated the interpolation assumptions. Trends in cumulative frequency provide key insights for soil consistency and foundation design. The results showed that the frequency of the SPT-N values for 1 m, 2 m, and 3 m depths were less scattered over the range of the x-axis in contrast to those at the depths of the SPT-N values beyond 3 m. Similarly, Figure 4i shows that the frequency of SPT-N values ranging between 0 and 4 became more pronounced at shallow strata (i.e., 1–3 m depth below the EGL), while at a higher depth, the frequency distribution was relatively scattered over a wide range of SPT-N values, which indicates that the study area adheres to the trend of high variance beyond 3 m depth. Figure 4j shows the maximum and minimum deviation that can occur from the average mean SPT-N values at a particular depth. In addition to various linear regression models, which were established and documented in Table 5, geotechnical data were utilized to forecast the SPT-N values by incorporating the depth factor (see Figure 4j). Additionally, the regression models demonstrated a high coefficient of determination, suggesting that it can be reliably used to estimate the SPT-N values at different depths during the early planning and design stages of future projects in the study area.

3.2. Selection of Optimum Interpolation Technique (OIT)

The various interpolation techniques were evaluated to assess the accuracy of prediction, and integrated with a geotechnical dataset based on SPT-N dispersed across the unexplored study region. The geotechnical database was joined via a linkage algorithm to associate it with the spatial location of the borehole for the precise formulation of spatial autocorrelations. The SPT-N database at a 1.0 m of depth below the EGL was evaluated via various geostatistical and spatial analyst approaches to access and quantify the accuracy for the subsequent development of SMs to ensure its practical implementation in the field. The spatial analyst approach entails the assessment of the IDW, spline, and radial basis function, while the kriging, diffusion, and polynomial interpolations under the geostatistical analyst approach were evaluated in terms of CPMs that involve the evaluation of the root mean square error (RMSE), mean absolute error (MAE), and correlation coefficient (CC), respectively. The optimal interpolation method was chosen for the further development of SMs based on the SPT-N values and soil types at various depth intervals, specifically from 1 to 10 m below EGL at 1 m intervals. During the extensive iterations of various input variables, the spatial analyst-based IDW found, with a power parameter of ‘2’, integrated with five neighboring boreholes, accounting for the radius of influence, and was identified to be optimal with the RMSE, MAE, and CC ranges between 0.57 and 0.98 (Figure 5). On the other hand, the RMSE, MAE, and CC were found to be in the range between 0.73 and 2.5 for the spline, radial basis function, kriging, and diffusion, with the least desired techniques identified, to the polynomial interpolation technique with a higher magnitude of error and the least association with the actual field database. Based on the CPM values, the IDW interpolation technique based on the spatial analyst approach was subsequently selected for the development of geotechnical soil maps for the unexplored study region.

3.3. Development of Soil Maps

Geotechnical data were gathered from various locations across the district, with sampling points uniformly distributed across all three administratively controlled tehsils (Figure 3). The site coordinates, elevation above mean sea level, SPT-N values, and soil type at different depths were numerically recorded from the GIRs and imported into ArcGIS as input data. ArcMap was then utilized to develop soil maps using the Spatial Analyst tool with the IDW interpolation technique. Based on a detailed evaluation of different interpolation methods, the IDW spatial interpolation approach was selected, which estimates the value at an unknown point as a weighted average of values from nearby points within a specified cut-off distance or from a set number of closest points. For this analysis, an IDW power parameter of “2” was applied.

3.4. Development of SPT-N Maps

Figure 6 presents the soil resistance against compressive loading at depths of 1 m, 2 m, 3 m, 4 m, 5 m, 6 m, 8 m, and 10 m under the EGL. These SMs indicate the soil resistive capacity and soil stability against the loading at different intervals of the stratum. The classification of SPT-N values was based on the soil stiffness as suggested for soft, medium, and stiff soil. Eight maps were developed to analyze the stratigraphical changes of SPT-N across the shallow and deep intervals of depth. The observed ranges of SPT-N values, namely 1–16 for 0–3 m, 1–34 for 3–6 m, and 5–39 for 6–10 m, revealed increasing resistance with depth, consistent with stratigraphical densification. At a depth of 1.0 m, the majority of the area consists of soft zones with SPT-N values ranging from 1 to 4, indicating weak strata unsuitable for foundation placement. As the depth increases, these soft zones diminish, giving way to stiffer strata. For instance, at a 2 m depth, the SPT-N values increased to a range of 5 to 8, reflecting improved soil strength. This transition from soft to stiff behavior became more pronounced with depth; by 10 m, the majority of the area was characterized by significantly higher SPT-N values ranging from 9 to 16 and exceeding 16.

Furthermore, as the study area is prone to peak ground acceleration in the range pf 0.16–0.24 as per the building code of Pakistan, the analysis incorporated the shear wave velocity in the soil maps. The predominant vs values ranged from 83 to 272 m/s for shallow depths (0–3 m), while depths beyond 3 m exhibited vs values ranging from 165 to over 272 m/s. This variation highlights the differences in the seismic wave propagation characteristics, which are crucial for assessing the seismic site response. It is pertinent to mention that shallow soil layers exhibit weak zones with low compressive resistance to loading, necessitating careful consideration during the early stages of construction to ensure the stability and reliability of foundation designs for future building infrastructure projects.

3.5. Soil Maps Based on Geotechnical Soil Type

SMs were created based on soil types identified at different stratigraphical layers of the study area. Soil types were classified according to the Unified Soil Classification System (USCS) and assigned numerical codes as follows: (1) CL, lean clay; (2) CL-ML, silty clay; (3) ML, silt; (4) SP-SM, poorly graded sand with silt; (5) SP, poorly graded sand. Based on the observed pattern of soil variation together with depth, nine maps were created, as shown in Figure 7. The depth intervals of the soil maps were 0 m, 1 m, 2 m, 3 m, 4 m, 5 m, 6 m, 8 m, and 10 m, respectively. The developed SMs demonstrated the predomination of CL and CL-ML soil in the shallow stratification (i.e., up to 3 m), which is consistent with the alluvial deposits formed by the sedimentation from the Chenab River and its tributaries, contributing to the fine-grained interbedded stratigraphical layers. Between 4 m and 6 m, the stratigraphical composition shifted significantly, with SP and SP-SM soils becoming more prevalent, indicating the presence of coarser-grained alluvial layers, possibly reflecting historical fluctuations in sediment deposition regimes.

Moreover, the deeper layers exhibited a transitional behavior, where the SP-SM soils were gradually predominated by CL or CL-ML soils, demonstrating a mix of fluvial and lacustrine depositional environments influenced by historical river migration and flooding events. Figure 7 illustrates this transition, highlighting the complexity of soil layering and the influence of the region’s geological history on its geotechnical characteristics.

3.6. Validation of Soil Maps

To assess the accuracy of the SMs at specific depths and locations, actual SPT-N values and soil types were compared with interpolated predictions generated via the optimized IDW technique. The borehole points reserved for validation were well-distributed across the study region to ensure comprehensive coverage, revealing only minor differences between the observed and predicted values, confirming the reliability and robustness of the IDW method for the development of SMs. Such SMs are critical to providing regional variation considering the regional variation in geological and geotechnical variation other than local anomalies. Figure 8 demonstrates the scatter plot of the field observations and predicted values, showing a more concentrated alignment at shallow depths (1–3 m below the EGL) that gradually became more dispersed at depths beyond 3.0 m. It is pertinent to mention here that the maximum variability between the actual and predicted SPT-N values in the first 3 m ranged between ±3. In order to further analyze the strength of the predicted SPT-N values and soil type, the RMSE, MAE, and CC were computed for various depths, as can be seen in Table 6.

3.7. Generalized Geotechnical Zoning for Lightly-Loaded Structures

Laying the foundation of lightly-loaded structures within the first 3 m of the EGL is a common practice in Pakistan. Therefore, it is necessary to develop generalized geotechnical zones for the SPT-N value for the average of the first 3 m overburden, which provides a quick estimation of the soil parameters and bearing capacity for the local practitioners dealing with lightly-loaded structures, keeping in mind the economical constraints. In this context, an average SPT-N-based SM was incorporated based on the average SPT-N values of the first 3 m below the EGL, as seen in Figure 9. The stratification reflects the influence of the depositional environment and overconsolidation ratio, which directly impact the soil stiffness and strength properties. A total of three zones were proposed for soft, medium, and stiff consistency soils. For the top 3 m, the average SPT-N value in zone-I was 1–4; for zone-II, it was 5–8; and for zone-III, 9–16. Table 7, Table 8 and Table 9 present the average subsurface soil parameters based on field and laboratory investigations for the design of shallow foundations.

Spread footings such as square footings are more appropriate for lightly-loaded structures. In addition, for these types of footing, the influence zone falls within 2B, with an average footing width ranging between 0.5 and 2.5 m, which is a common practice in Pakistan. Therefore, the developed SMs up to 10 m depth well incorporated the influence zones of spread footing that fells within the above-mentioned range. Furthermore, the ABC for spread footing was evaluated based on the shear and settlement criteria. For the shear criteria, the ABC was calculated against the factor of safety (FOS) 3.0, while for the settlement criteria, the ABC against 25 mm settlement was assessed. The bearing capacities were computed by varying the foundation width between 0.5 and 2.5 m by keeping the depth of footing to 1 m. The Terzaghi equation was used to evaluate the safe bearing capacity in terms of the shear strength of the soil whereas the settlement was computed by incorporating the finite element method software known as the Plaxis 3D Foundation. A typical illustration of the simulation results on the footing model is presented in Figure 10. The governing values of bearing capacity are the lowest of two criteria (i.e., shear and settlement). Moreover, as the study area falls under the medium-intensity seismic zone, a third criterion based on seismicity was also incorporated into the current study. The governing equation used to calculate the seismic allowable bearing capacity is presented in Equation (2) [30]. Furthermore, a comparison was drawn between the ABC computed based on the SPT-N values and shear wave velocity.

(2)$q_a=0.025\gamma V_s\beta .$

where $q_a$ is the allowable bearing capacity, $\gamma$ is the unit weight of soil, vs is the shear wave velocity, and $\beta$ is the correction factor for footing width (B), as shown in Equations (3)–(5). The suggested values of $\beta$ for various footing widths (B) are as follows:

(3)$\beta =1.0 \mathrm{for} 0.0< \mathrm{B} < 1.2$

(4)$\beta =1.13-0.11B. \mathrm{for} 1.2< \mathrm{B} < 3$

(5)$\beta =0.83-0.01B \mathrm{for} 3.0 < \mathrm{B} < 12.0$

Figure 10 presents the ABC curves for the spread footing based on the SPT-N values and shear wave velocity for all three zones. In general, the ABC for spread footing based on the SPT-N values was higher than the seismic ABC for all three zones (i.e., I, II, and III). Furthermore, zone-I exhibited a lower ABC for spread footing due to the presence of soft strata, while zones-II and -III presented reasonably good strata for spread footing. Among all of the zones, zone-III presented the highest ABC based on SPT-N values. In addition, the ABC of spread footing in zones-I and -II was controlled by the shear criteria at a footing width <2 m; conversely, beyond 2 m, the settlement criteria prevailed. In the case of zone-III, the ABC of spread footing was controlled by the settlement criteria. On the other hand, the ABC based on the shear wave velocity for zones-I, -II, and -III was found to be in the range of 28–40%, 45–51%, and 38–48%, respectively, lower than the ABC based on the SPT-N values. This implies that seismicity has a major influence on the bearing capacity calculations and must be incorporated into the design. For comparison purposes, the ABC was further compared with the study pertaining to the Faisalabad region. It was found that the ABC curves for the spread footing of Faisalabad district were relatively higher compared to the district of Gujrat and predominantly comprised of cohesionless soil with shear criteria that significantly controlled the ABC of the district, as can be seen in Figure 11a–c. It is pertinent to mention here that in medium- to high-rise buildings where structural loads are significant and settlement control is critical, mat footings could be a more appropriate choice to ensure uniform load distribution and reduce the risk of differential settlement.

4. Conclusions

This article focused on the development and optimization of geotechnical SMs of an unexplored region to analyze the spatial variability in geotechnical settings, which were further used to establish recommendations and suggest potentially appropriate areas for the construction of lightly-loaded structures with shallow footings. The key outcomes of this study are as follows:

• The extensive iteration, the spatial analyst approach-based IDW, was found to be the most effective interpolation technique (OIT) with an RMSE, MAE, and CC ranging from 0.57 to 0.98. The least desirable method was polynomial interpolation, with an RMSE, MAE, and CC ranging from 0.73 to 2.5.

• The SMs showed that the top 3 m deposit was mostly composed of lean clay, clayey silt, and/or silt with an average SPT-N value of 1–16 and a shear wave velocity of 83–272 m/s, resulting in soft to hard strata. For construction in soft strata regions, the thickness, depth, and seismicity factors should be considered for a safe and cost-effective design. The strata between 3 and 6 m were mostly comprised of poorly graded sand with silt, based on the soil types. Beyond 6 m, this trend changed to lean clay, clayey silt, or silt, with typical mean SPT-N values ranging between 9–16 and >16, which is suitable ground support for most civil engineering constructions.

• The comprehensive validation of the developed SMs at each interval, quantified through CPMs, demonstrated RMSE and MAE ranges of 0.57–1.65 and 0.30–0.59, respectively, for SPT-N and soil type. Additionally, the CC ranged from 0.93 to 0.98, ensuring the reliability and applicability of these models for preliminary design phases in construction projects.

• Based on the spatial variation, three geotechnical zones were identified based on SPT-N values up to a 3 m depth, with averages of 1–4, 5–8, and 9–16 in zones-I, -II, and -III, respectively. The shallow footing designs in zones-I and-II were governed by shear criteria, while the settlement criteria dominated in the zone-III designs. The ABC with seismic consideration was found to be 28–51% lower than the ABC based on the shear and settlement criteria across all zones, emphasizing the significance of dynamic soil properties in the design.

• The SMs and geotechnical zone classifications provide valuable guidelines for local practitioners and researchers offering insights into key geotechnical parameters along with geotechnical variability across the study region, which is critical for the initial phase of project planning, site selection, and the design of a preemptive response system.

Author Contributions

Conceptualization, N.I.; Methodology, N.I.; Software, N.I.; Data curation, N.I., Investigation, N.I.; Project administration, N.I.; Writing—original draft preparation, N.I.; Conceptualization, Z.I.; Methodology, Z.I.; Validation, Z.I.; Software, Z.I.; Investigation, Z.I.; Review and editing, Z.I.; Visualization, Z.I.; Supervision, N.Z.; Formal analysis, N.Z.; Writing—review and editing, N.Z.; Formal analysis, Z.u.R.; Visualization, Z.u.R.; Writing—review and editing, Z.u.R.; Validation, H.I.; Formal analysis, H.I.; Visualization, H.I.; Writing—review and editing, H.I.; Formal analysis, A.I.; Writing—review and editing, A.I.; Formal analysis, M.H.; Writing—review and editing, M.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

The authors would like to express their sincere gratitude to Quanzhou University of Information Engineering, Tongji University, the Communication and Works Department, the University of Portsmouth, UK, the University of Engineering and Technology Lahore, and Shenzhen University for their invaluable support and contributions to this study.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Spatial details of the district of Gujrat. (a) Location of the district of Gujrat on a map of Pakistan, (b) district of Gujrat map in Punjab province, (c) district of Gujrat map with the administrative controlled tehsils Sarai Alamgir, Kharian, and Gujrat, and (d) average climatic conditions.

Figure 2. (a) Geological map of the study area, (b) UNESCO/FAO soil classification.

Figure 4. (a–f) SPT-N histogram at corresponding depths. (i) Comparison of SPT-N distribution at various depths. (j) Descriptive changes of SPT-N with depth.

Figure 4. (a–f) SPT-N histogram at corresponding depths. (i) Comparison of SPT-N distribution at various depths. (j) Descriptive changes of SPT-N with depth.

Figure 5. CPMs of the various interpolation techniques.

Figure 6. SPT-N variation up to a 10 m depth at each 1 m interval of the stratum.

Figure 7. Geotechnical properties of soil based on soil type.

Figure 8. Validation of SPT-N and soil type.

Figure 9. Average SPT-N zonation map.

Figure 10. Illustration of the simulation results computed from Plaxis 3D. (a) Isometric view of the model after the application of the load, (b) cross-sectional view of the model, (c) cumulative settlement.

Figure 11. Allowable bearing capacity curves for a spread foundation for zones-I, -II, and -III.

Table 1

Distribution of land covers of the district of Gujrat.

Distributed Land Cover for the District of Gujrat
	km2	%
Crop irrigated	0.83	0.026
Crop marginal	2134	66.85
Flood plains	0.09	0.002
Crop rainfed	284	8.89
Forest	292	9.152
Vegetation	18.08	0.59
Rangelands	160	5.01
Built-up areas	212	6.64
Bare areas	20	0.62
Wet areas	71	2.22
Total	3192 km2	100

Table 3

Statistical descriptors of the SPT-N data.

Statistical Descriptor of SPT-N Data
Depth (m)	1	2	3	4	5	6	8	10
Mean value	4.66	6.06	7.5	9.57	11.62	13.023	15.65	20.81
St. deviation	2.34	2.63	3.37	4.704	5.10	5.66	5.67	5.76
Minimum	1	1	1	1	3	4	5	11
Maximum	13	15	16	31	34	36	39	39
Data count	130	130	130	130	130	130	128	128

Table 4

ANOVA analysis of SPT-N data.

Analysis of Variance of SPT-N Data
Groups/Depth		Count	Sum	Average	Variance
1 m		130	470.65	4.66	5.519
2 m		130	612.43	6.063	6.93
3 m		130	758.28	7.5	11.37
4 m		130	966.57	9.57	22.14
5 m		130	1174.1	11.62	26.09
6 m		130	1315.33	13.02	32.09
8 m		128	1424.42	15.65	32.14
10 m		128	1706.9	20.81	33.21
Results of ANOVA analysis for the district of Gujrat
Source of Variation	SS	dF	MS	F-stat	p-value
Between groups	18,270.04	7	2610.006	124.6194	0
Within groups	16,168.62	772	20.9438
Total	34,438.67	779

Table 5

Formulation regression correlations with depth.

Regression Formulations of SPT-N Values with Depth
Profile	Correlation	R2
Average	N = 1.753 (D) + 2.565	0.992
Average − S.D	N = 2.172(D) + 4.926	0.989
Average + S.D	N= 1.334(D) + 0.203	0.956

where N = SPT-N value, D = depth.

Table 6

Statistical assessment of the predicted and field values.

Depth	SPT-N			Soil Type
m	RMSE	MAE	CC	RMSE	MAE	CC
1	0.57	0.36	0.98	0.41	0.33	0.98
2	0.79	0.60	0.98	0.44	0.42	0.98
3	0.86	0.61	0.98	0.54	0.46	0.98
4	0.92	0.80	0.97	0.57	0.49	0.98
5	1.10	0.90	0.97	0.59	0.51	0.97
6	1.16	0.92	0.96	0.48	0.49	0.97
8	1.50	1.28	0.94	0.3	0.32	0.94
10	1.65	1.32	0.93	0.3	0.38	0.94

Table 7

Foundation design parameters for zone-I.

Soil Identification	Depth	Bulk Density (γ)	Undrained Cohesion (c)	Coefficient of Compressibility (mv)
	(m)	(kg.cm3)	(kg/cm3)	(cm2/kg)
CL	0–1	0.0016	0.2	0.02
CL-ML	1–3	0.00165	0.3	0.017
CL-ML	3–4.5	0.0017	0.45	0.015
CL-ML	4.5–10	0.00175	0.6	0.013

Table 8

Foundation design parameters for zone-II.

Zone #2
Soil Identification	Depth	Bulk Density (γ)	Undrained Cohesion (c)	Coefficient of Compressibility (mv)	Angle of Internal Friction (φ)
	(m)	(kg.cm3)	(kg/cm3)	(cm2/kg)	(Degrees)
CL-ML	0–3	0.00175	-	0.0018	-
SP-SM/SP	3–5.5	0.00185	0.7	-	28
CL/CL-ML	5.5–10	0.0019	-	0.0012	-

Table 9

Foundation design parameters for zone-III.

Zone # 3
Soil Identification	Depth	Bulk Density (γ)	Undrained Cohesion (c)	Coefficient of Compressibility (mv)	Angle of Internal Friction (φ)
	(m)	(kg.cm3)	(kg/cm3)	(cm2/kg)	(Degrees)
SP-SM	0–4.5	0.0018	-	-	30
ML/CL-ML	4.5–7.5	0.00185	0.75	0.0012	-
CL/CL-ML	7.5–10	0.0019	0.85	0.001	-

References

1. Otake, Y.; Honjo, Y. Challenges in geotechnical design revealed by reliability assessment: Review and future perspectives. Soils Found.; 2022; 62, 101129. [DOI: https://dx.doi.org/10.1016/j.sandf.2022.101129]

2. Pinheiro, M.; Vallejos, J.; Miranda, T.; Emery, X. Geostatistical simulation to map the spatial heterogeneity of geomechanical parameters: A case study with rock mass rating. Eng. Geol.; 2016; 205, pp. 93-103. [DOI: https://dx.doi.org/10.1016/j.enggeo.2016.03.003]

3. Skinner, H.; Toll, D.G.; Higgins, K.; Brown, M.; Burland, J. ICE Manual of Geotechnical Engineering Volume 1: Geotechnical Engineering Principles, Problematic Soils and Site Investigation; Emerald Group Publishing: Bingley, UK, 2023; Volume 1.

4. Von der Tann, L.; Sterling, R.; Zhou, Y.; Metje, N. Systems approaches to urban underground space planning and management–A review. Undergr. Space; 2020; 5, pp. 144-166. [DOI: https://dx.doi.org/10.1016/j.undsp.2019.03.003]

5. Tougwa, F.N. Some major challenges faced by civil engineering professionals in the execution of their profession and the impact of the challenges to the environment, society and economy of developing countries. Society and Economy of Developing Countries. Curr. Trends Civ. Struct. Eng.; 2020; 5, 622. [DOI: https://dx.doi.org/10.33552/CTCSE.2020.05.000622]

6. Hakro, M.R.; Kumar, A.; Ali, M.; Habib, A.F.; de Azevedo, A.R.G.; Fediuk, R.; Sabri, M.M.S.; Salmi, A.; Awad, Y.A. Numerical analysis of shallow foundations with varying loading and soil conditions. Buildings; 2022; 12, 693. [DOI: https://dx.doi.org/10.3390/buildings12050693]

7. Imran, M.; Ali, M.; Janjua, M.S. Urban development challenges under the China-Pakistan Economic Corridor (CPEC). Int. Perspect. Belt Road Initiat.; 2021; pp. 81-99.

8. Lundström, K.; Larsson, R.; Dahlin, T. Mapping of quick clay formations using geotechnical and geophysical methods. Landslides; 2009; 6, pp. 1-15. [DOI: https://dx.doi.org/10.1007/s10346-009-0144-9]

9. Ijaz, Z.; Zhao, C.; Ijaz, N.; Rehman, Z.U.; Ijaz, A. Spatial mapping of geotechnical soil properties at multiple depths in Sialkot region, Pakistan. Environ. Earth Sci.; 2021; 80, 787. [DOI: https://dx.doi.org/10.1007/s12665-021-10084-z]

10. Kim, H.S.; Ji, Y. Three-dimensional geotechnical-layer mapping in Seoul using borehole database and deep neural network-based model. Eng. Geol.; 2022; 297, 106489. [DOI: https://dx.doi.org/10.1016/j.enggeo.2021.106489]

11. Ijaz, Z.; Zhao, C.; Ijaz, N.; Rehman, Z.U.; Ijaz, A.; Junaid, M.F. Geospatial modeling of heterogeneous geotechnical data using conventional and enhanced conception of modified Shepard method-based IDW algorithms: Application and appraisal. Bull. Eng. Geol. Environ.; 2023; 82, 428. [DOI: https://dx.doi.org/10.1007/s10064-023-03435-6]

12. Aziz, M.; Khan, T.A.; Ahmed, T. Spatial interpolation of geotechnical data: A case study for Multan City, Pakistan. Geomech. Eng.; 2017; 13, pp. 475-488.

13. Kamal, M.A.; Arshad, M.U.; Khan, S.A.; Zaidi, B.A. Appraisal of geotechnical characteristics of soil for different zones of Faisalabad (Pakistan). Pak. J. Eng. Appl. Sci.; 2015.

14. Ijaz, Z.; Zhao, C.; Ijaz, N.; Rehman, Z.U.; Ijaz, A. Development and optimization of geotechnical soil maps using various geostatistical and spatial interpolation techniques: A comprehensive study. Bull. Eng. Geol. Environ.; 2023; 82, 215. [DOI: https://dx.doi.org/10.1007/s10064-023-03244-x]

15. Islam, M.A.; Hasan, M.A.; Farukh, M.A. Application of GIS in general soil mapping of Bangladesh. J. Geogr. Inf. Syst.; 2017; 9, pp. 604-621. [DOI: https://dx.doi.org/10.4236/jgis.2017.95038]

16. Chen, S.; Arrouays, D.; Mulder, V.L.; Poggio, L.; Minasny, B.; Roudier, P.; Libohova, Z.; Lagacherie, P.; Shi, Z.; Hannam, J. et al. Digital mapping of GlobalSoilMap soil properties at a broad scale: A review. Geoderma; 2022; 409, 115567. [DOI: https://dx.doi.org/10.1016/j.geoderma.2021.115567]

17. Vaysse, K.; Lagacherie, P. Evaluating digital soil mapping approaches for mapping GlobalSoilMap soil properties from legacy data in Languedoc-Roussillon (France). Geoderma Reg.; 2015; 4, pp. 20-30. [DOI: https://dx.doi.org/10.1016/j.geodrs.2014.11.003]

18. Almasi, A.; Jalalia, A.; Toomanian, N. Using OK and IDW methods for prediction the spatial variability of a horizon depth and OM in soils of Shahrekord, Iran. J. Environ. Earth Sci.; 2014; 4, pp. 17-27.

19. Hasterok, D.; Halpin, J.A.; Collins, A.S.; Hand, M.; Kreemer, C.; Gard, M.G.; Glorie, S. New maps of global geological provinces and tectonic plates. Earth-Sci. Rev.; 2022; 231, 104069. [DOI: https://dx.doi.org/10.1016/j.earscirev.2022.104069]

20. Qian, Y.; Xiao, L.; Wang, Q.; Head, J.W.; Yang, R.; Kang, Y.; van der Bogert, C.H.; Hiesinger, H.; Lai, X.; Wang, G. et al. China’s Chang’e-5 landing site: Geology, stratigraphy, and provenance of materials. Earth Planet. Sci. Lett.; 2021; 561, 116855. [DOI: https://dx.doi.org/10.1016/j.epsl.2021.116855]

21. Aksoy, T.; Dabanli, A.; Cetin, M.; Senyel Kurkcuoglu, M.A.; Cengiz, A.E.; Cabuk, S.N.; Agacsapan, B.; Cabuk, A. Evaluation of comparing urban area land use change with Urban Atlas and CORINE data. Environ. Sci. Pollut. Res.; 2022; 29, pp. 28995-29015. [DOI: https://dx.doi.org/10.1007/s11356-021-17766-y]

22. Raza, M.; Farooqi, A.; Niazi, N.K.; Ahmad, A. Geochemical control on spatial variability of fluoride concentrations in groundwater from rural areas of Gujrat in Punjab, Pakistan. Environ. Earth Sci.; 2016; 75, pp. 1-16. [DOI: https://dx.doi.org/10.1007/s12665-016-6155-7]

23. van der Elst, N.J.; Hardebeck, J.L.; Michael, A.J.; McBride, S.K.; Vanacore, E. Prospective and retrospective evaluation of the US Geological Survey Public aftershock forecast for the 2019–2021 Southwest Puerto Rico Earthquake and aftershocks. Seismol. Soc. Am.; 2022; 93, pp. 620-640.

24. Mantel, S.; Dondeyne, S.; Deckers, S. World reference base for soil resources (WRB). Goss, Margaret Oliver Encyclopedia of Soils in the Environment; 2nd ed. Michael, J. Elsevier: Amsterdam, The Netherlands, 2023; pp. 206-217.

25. Buenemann, M.; Coetzee, M.E.; Kutuahupira, J.; Maynard, J.J.; Herrick, J.E. Errors in soil maps: The need for better on-site estimates and soil map predictions. PLoS ONE; 2023; 18, e0270176. [DOI: https://dx.doi.org/10.1371/journal.pone.0270176] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36630410]

26. Arrouays, D.; Mulder, V.L.; Richer-de-Forges, A.C. Soil mapping, digital soil mapping and soil monitoring over large areas and the dimensions of soil security–A review. Soil Secur.; 2021; 5, 100018. [DOI: https://dx.doi.org/10.1016/j.soisec.2021.100018]

27. Siddique, M.S.; Schwarz, J. Elaboration of multi-hazard zoning and qualitative risk maps of Pakistan. Earthq. Spectra; 2015; 31, pp. 1371-1395. [DOI: https://dx.doi.org/10.1193/042913EQS114M]

28. Khan, S.A.; Pilakoutas, K.; Hajirasouliha, I.; Garcia, R.; Guadagnini, M. Seismic risk assessment for developing countries: Pakistan as a case study. Earthq. Eng. Eng. Vib.; 2018; 17, pp. 787-804. [DOI: https://dx.doi.org/10.1007/s11803-018-0476-3]

29. Vabalas, A.; Gowen, E.; Poliakoff, E.; Casson, A.J. Machine learning algorithm validation with a limited sample size. PLoS ONE; 2019; 14, e0224365. [DOI: https://dx.doi.org/10.1371/journal.pone.0224365]

30. Tezcan, S.S.; Ozdemir, Z.; Keceli, A. Seismic technique to determine the allowable bearing pressure for shallow foundations in soils and rocks. Acta Geophys.; 2009; 57, pp. 400-412. [DOI: https://dx.doi.org/10.2478/s11600-008-0077-z]