High-quality mesh generation is critical in the finite element analysis of displacements and stabilities of geological bodies. In this paper, we propose a clustering-based bubble method for generating high-quality tetrahedral meshes of geological models. The proposed bubble method is conducted based on the spatial distribution of the point set of given surface meshes using the clustering method. First, the inputted geological models consisting of triangulated surface meshes are divided into several parts based on spatial distribution of point set, which can be used for the determination of the positions and radii of initial bubbles. Second, a procedure based on distance of nearby bubbles is used to obtain the initial size of bubbles. Third, by enforcing the forces acting on bubbles, all bubbles inside the 3D domain reach an equilibrium state by the motion control equations. Finally, the center nodes of the bubbles can form a high-quality node distribution in the domain, and then the required tetrahedral mesh is generated. Comparative benchmarks are presented to demonstrate that the proposed method is capable of generating highly well-shaped tetrahedral meshes of geological models.