Constructing a comprehensive three-dimensional (3D) geological model using borehole data is crucial for establishing a cohesive comprehension of complex subsurface geology, which is essential for groundwater studies. This study aims to introduce a novel geostatistical method–interval kriging–to efficiently conduct 3D lithologic modeling with multiple indicators using borehole data. The interval kriging is a best linear unbiased estimator utilizing irregular interval supports. The interval kriging can model 3D zonal or geometric anisotropies between a horizontal plane and a vertically orthogonal axis to the plane. To address the nonconvexity of the estimation variance, an additional regularization term is incorporated into minimizing the estimation variance. The minimization problem is solved by a global-local genetic algorithm embedded with gradient methods to obtain kriging weights and kriging length. The numerical and real-world case studies demonstrate that the interval kriging is more computationally efficient than the traditional 3D kriging method because the covariance matrix is largely reduced without sacrificing borehole data. Moreover, the interval kriging produces more realistic geologic characteristics than the 2.5D kriging, while conditional to spatial borehole data. Besides, the interval kriging produces more certainty in the classification problem than 2.5D and 3D kriging methods because the regularization term maximizes the margin of the hyperplane for classification. To conclude, the interval kriging is an effective and efficient 3D modeling technique capable of capturing the 3D structural complexity while significantly reducing computational time.