Spatial variability (randomness, correlation, and singularity) within the geotechnical parameters of complicated geological movements influences the estimation quality that depends on how well mathematical tools can account for variability through limited observations of a spatial field. Classical statistical methods depict randomness well, but cannot account for the problems associated with spatial correlations. Geostatistical methods such as ordinary Kriging (OK), universal Kriging (UK), and co-Kriging (CK) can produce predictions based on spatial auto-correlation and cross-correlation, but are always accompanied by average smoothing effects; a local singularity created by nonlinear geo-processes, therefore, requires special methods to be properly evaluated. In this study, a shallow load-bearing stratum of silt clay (length = 525 m, width = 80 m) at the former 2010 Expo Park in Shanghai was explored by performing 42 borehole laboratory experiments, which provided the key geotechnical parameters: the cohesion coefficient (\( C \), in kPa), the friction angle (\( \varphi \), in o), and the compression modulus (\( E_{\text{S}} \), in MPa). First, Kriging methods such as OK, UK, and CK estimated these geotechnical parameters, then a multi-fractal analysis was employed to measure the local singularity. Cross-validation illustrates that multi-fractal analysis has the ability to depict a local anomaly, and further that the auxiliary information utilized in CK improves spatial estimation accuracy.