## 1 Introduction  
Numerical simulation has been applied to solve many problems in various fields, and mesh generation always is a vital topic in finite element analysis. Among the mesh generation algorithms, Delaunay triangulation [1; 2; 3; 4], the advancing front technique [5], the mapping method [6] have been widely developed and applied due to the high quality of the generated mesh and easy implementation of the algorithms. Many researchers have focused on tetrahedral mesh generation methods [7; 8]. For example, [PERSON] [9; 10] proposed a novel link-Delaunay-triangulation method to achieve geometric and topological consistency. Indirect methods [11; 12] were proposed to convert Delaunay triangular meshes into quadrilateral meshes by combining adjacent pairs of triangles. Constrained triangulated surfaces can be used to mesh a three-dimensional geological model by using a series of tetrahedral meshes [13; 14]. Other methods [15; 16] for generating three-dimensional meshes were also proposed. Moreover, [PERSON] et al. [17] and [PERSON] [18] first proposed the bubble method for mesh generation of the non-manifold geometry. This method generates point sets and the Delaunay triangulation technique is used for mesh generation. [PERSON] [19] proposed an algorithm for three-dimensional bubble mesh generation. In this method, many points were inserted by the mapping-based method. [PERSON] [20] used bubble method to discretize two-parameter surfaces into a high-quality mesh. [PERSON] [21; 22] proposed a fast local mesh generation algorithm based on the node-set generated by this method and proposed a subdivided dynamic bubble system to shorten the calculation time. [PERSON] [23] presented an automatic triangular mesh generation method, based on bubble dynamics simulation and a modified Delaunay method to generate high-quality mesh on complex surfaces efficiently. [PERSON] [24] implemented bubble method for adaptive mesh generation in two and three dimensions. The geological model is a three-dimensional model representing the ground surface, strata, and faults, which consists of multiple geological blocks describing the structure and geometry of the stratum. Each block in the model has a different shape. The thin layer and the fault are very flat geometry, and the thick layer is a very thick geometry. Therefore, when the geological model is transformed into a numerical model, the geological model needs to be divided into the high-quality mesh to meet a balance between calculation accuracy and efficiency. 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 point set of given surface meshes using the clustering method. There are four main procedures in the 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 points. 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 node of the bubble can form a high-quality node distribution in the domain, and then the required tetrahedral mesh is generated. The novelty of the proposed method can be explained as follows. As mentioned above, most of the above bubble methods focus on two-dimensional cases and three-dimensional parameter surface. However, since geological structures and conditions are generally quite complex, a three-dimensional geological model cannot be well represented by parametric surfaces but can be well represented of triangulated surfaces. The proposed bubble method is specifically proposed to generate the high-quality meshes of geological models consisting of triangulated meshes.