The design and construction of underground structures are significantly affected by the distribution of geological formations. Prediction of the geological interfaces using limited data has been a difficult task. A multivariate adaptive regression spline (MARS) method capable of modeling nonlinearities automatically was used in this study to spatially predict the elevations of geological interfaces. Borehole data from two sites in Singapore were used to evaluate the capability of the MARS method for predicting geological interfaces. By comparing the predicted values with the borehole data, it is shown that the MARS method has a mean of root mean square error of 4.4 m for the predicted elevations of the Kallang Formation–Old Alluvium interface. In addition, the MARS method is able to produce reasonable prediction intervals in the sense that the percentage of testing data covered by 95% prediction intervals was close to the associated confidence level, 95%. More importantly, the prediction interval evaluated by the MARS method had a non-constant width that appropriately reflected the data density and geological complexity.