## 2 Materials and Methods
### The Mesh Generations for the HD Model  
Characterizing the interaction between fractures and rock matrix employs the tetrahedral mesh in the model to represent the rock matrix. The fractures are constructed by using triangle meshes in the simulation domain [16; 17; 18; 19]. In the mesh generation process, the fracture geometries must be identified first to define the inner boundaries. The fractures are groups of connected 2D plates in a 3D domain. Note that the identified fractures have excluded the isolated fractures for computational efficiency. Based on the identified fractures in the simulation domain, the triangle meshes for the fractures are generated, and the element sizes can vary with the complexity of the fracture connections. The 2D fracture meshes are the basis for generating the tetrahedral mesh for the 3D rock matrices.  
The challenge of the mesh generation process is to construct the 3D tetrahedral meshes, which share the existing triangular elements of the identified fractures (see Figure 1). The generated meshes are crucial to computational convergence and efficiency. This study proposes the following criteria to meet the requirements: (1) The integrity of the triangle meshes for fractures needs to be maintained, (2) The boundary recovery algorithm must be considered to account for the collinear or co-points in the processes of generating tetrahedral meshes [35], (3) The obtuse angles need to be avoided when generating triangles and tetrahedrons, and (4) The refinement of the meshes near interfaces between fractures and matrices is preferred to characterize the large discrepancy in hydraulic conductivity between fractures and rock matrices. These requirements have been issues in modeling flow and transport in fractured rocks [35; 36; 37; 38; 39].  
In this study, the generation of tetrahedral and triangular meshes is based on the concept of the Delaunay triangulation algorithm. The algorithm enables the generation of tetrahedral and triangular meshes by adding additional variables in three different directions. However, the tetrahedron might involve a four-point coplanar, a high acute angle, or a line with four points. In this study, the developed mesh generation program has a completed feature to loop over the bad connections of the elements and nodes. In general, procedures in the loops are to search fracture intersections, sequentially add nodes inside triangles to refine the mesh until the node distances reach the specified criteria, and smooth and adjust the triangles to close to acute triangles. Once an obtuse angle appears, the mesh generation model destroys the obtuse-angled triangular mesh with the adjacent triangles and regenerates the local meshes until all of the meshes are close to the acute triangles. The completed fracture meshes are the basis for generating 3D tetrahedral mesh in fractured rocks.  
Figure 1 is an example generated by the developed mesh generation model. In the example, there are two fractures with one collinear line in the center of the model domain (see Figure 1a). The fractures are composed of triangular meshes. Figure 1b shows the two fractures comprised of triangular meshes. In this case, the length of the triangular mesh is constant so that the sizes and areas of triangular elements are similar. With the available fracture meshes, the rest portions of the domain are filled with tetrahedral meshes. Figure 1c is the exterior view of the domain composed of fracture and rock matrix elements. Inside the simulation domain, the generated meshes consist of triangular and tetrahedral meshes (see Figure 1d). Note that the mesh generation concept is suitable to embed irregular objects inside the simulation domain when the volume surfaces of the objects are available. The surfaces of the objects could be treated as the inner boundaries, i.e., similar to the fractures. The generated meshes are then used in our HD model for simulating flow and particle tracking.