Voxel Interpolation of Geotechnical Properties and Soil Classification Based on Empirical Bayesian Kriging and Best-Fit Convergence Function

To support bearing capacity estimates, this study develops and tests a geoprocessing workflow for predicting soil properties using Empirical Bayesian Kriging 3D and a classification function. The model covers a 183 m × 185 m × 24 m site in Astana (Kazakhstan), based on 16 boreholes (15–24 m deep) and 77 samples. Eight geotechnical properties were mapped in 3D voxel models (812,520 voxels at 1 m × 1 m × 1 m resolution): cohesion (c), friction angle (φ), deformation modulus (E), plasticity index (PI), liquidity index (LI), porosity (e), particle size (PS), and particle size distribution (PSD). Stratification patterns were revealed with ~35% variability. Maximum φ (34.9°), E (36.6 MPa), and PS (1.29 mm) occurred at 8–16 m; c (33.1 kPa) and PSD peaked below 16 m, while PI and e were elevated in the upper and lower strata. Strong correlations emerged in pairs φ-E-PS (0.91) and PI-e (0.95). Classification identified 10 soil types, including one absent in borehole data, indicating the workflow’s capacity to detect hidden lithologies. Predicted fractions of loams (51.99%), sandy loams (22.24%), and sands (25.77%) matched borehole data (52%, 26%, 22%). Adjacency analysis of 2,394,873 voxel pairs showed homogeneous zones in gravel–sandy soils (28%) and stiff loams (21.75%). The workflow accounts for lateral and vertical heterogeneity, reduces subjectivity, and is recommended for digital subsurface 3D mapping and construction design optimization.