Numerical Assessment of the Hybrid Approach for Simulating Three-Dimensional Flow and Advective Transport in Fractured Rocks

: This study presents a hybrid approach for simulating ﬂow and advective transport dynamics in fractured rocks. The developed hybrid domain (HD) model uses the two-dimensional (2D) triangular mesh for fractures and tetrahedral mesh for the three-dimensional (3D) rock matrix in a simulation domain and allows the system of equations to be solved simultaneously. This study also illustrates the HD model with two numerical cases that focus on the ﬂow and advective transport between the fractures and rock matrix. The quantitative assessments are conducted by comparing the HD results with those obtained from the discrete fracture network (DFN) and equivalent continuum porous medium (ECPM) models. Results show that the HD model reproduces the head solutions obtained from the ECPM model in the simulation domain and heads from the DFN model in the fractures in the ﬁrst case. The particle tracking results show that the mean particle velocity in the HD model can be 7.62 times higher than that obtained from the ECPM mode. In addition, the developed HD model enables detailed calculations of the ﬂuxes at intersections between fractures and cylinder objects in the case and obtains relatively accurate ﬂux along the intersections. The solutions are the key factors to evaluate the sources of contaminant released from the disposal facility. Contributions:


Introduction
Modeling groundwater flow and solute transport in fractured rocks is essential in engineering, geology, and hydrogeology studies. Many applications, including gas/oil transport within the fractured system in the shale formation, artificial fractures by hydraulic fracturing methods for increasing the connectivity of a fracture system, and enhanced geothermal system (EGS) in hot and dry formations, have been widely investigated based on the knowledge of flow and transport in the fractured rocks [1][2][3][4]. In the procedures of site characterization investigation for the spent nuclear fuel final disposal, the solute transport in the fractured systems plays an essential role because the critical issue has been focused on the release of radionuclides from the deposition holes (DHs) in the disposal facility to the human environment, i.e., the biosphere [5,6]. The concept of spent nuclear fuel final disposal is mainly for radioactive wastes. The radioactive wastes are encapsulated in tight, corrosion-resistant, and load-bearing copper canisters deposited in the DHs. In the DHs, the canisters are protected by buffers to prevent canister corrosion. The high flow velocity in fractures could erode the buffers, and the flow might contact the canister at relatively high flux. Such high-velocity water carries corrosive materials and could obtain the solutions of DFN and ECPM models for comparison purposes [25][26][27][28][29]. Based on the same flow conditions, two synthetic cases were designed to evaluate the developed HD model. The first simple case (Case I) was composed of three intersected fractures in a simulation domain. In the second case (Case II), a cylinder object was embedded in the fractured rock to assess the flow and advective transport in the complex flow system. In this study, the generated computation meshes for different cases needed to meet the stability requirement for flow simulation. A systematic comparison was conducted to quantify the differences in heads and particle paths obtained from three different modeling concepts.

Materials and Methods
This study develops the HD model for modeling flow and advective transport in fractured rocks. Many previous studies employed the FracMan and DarcyTools for different applications. FracMan is a commercial software developed by Golder Associates and widely used software for fracture generation, flow and advective transport simulation based on the DFN concept. However, DarcyTools is a non-commercial software developed by Svensk Kärnbränslehantering AB and an integrated software package used for the same purposes based on the ECPM concept. The detailed theoretical concept and the numerical algorithms can be obtained in the user's guide or programming guideline [25][26][27][28][29]. Here, we focus on presenting the governing equations and the associated numerical algorithm for the HD model.

Flow Equation
This study assumes that the virtual fracture includes two parallel plates that enable the fluid to pass through the aperture of the plates at relatively high velocity. The virtual space between the two parallel plates consists of the void space. Therefore, the virtual fracture can be treated as a two-dimensional (2D) porous medium [21,22,[30][31][32][33][34]. If the groundwater is assumed as an incompressible fluid, Darcy's law can be applied to formulate the steadystate flow equation: where u is the Darcy velocity (m/s), K is the hydraulic conductivity (m/s), h is the hydraulic head (m), q is the source/sink term (1/s), and notation Λ is the equidimensional model domain (m 3 ). The boundary conditions on the boundaries ∂Λ have the following types: u· n| ∂Λ u = u, on ∂Λ u (4) where h is the hydraulic head on the boundary ∂Λ h (m) and u is the Darcy velocity perpendicular to the boundary surface ∂Λ u (m/s) concerning the outer unit normal vector n (-). This study assumes that a virtual fracture has two main characteristics, including (1) The width of the fracture is uniform in an element and can be represented as fracture aperture, which is extremely smaller than the fracture size, and (2) The hydraulic conductivity of a virtual fracture is relatively higher than that of the rock matrices. Therefore, the fracture can be composed of 2D triangular meshes, whereas the host rock is composed of 3D tetrahedron meshes. The 2D triangular meshes share the nodes and faces of 3D tetrahedron meshes. Such approximation can significantly reduce the required small meshes near the interfaces between fractures and the rock matrix. A mass balance equation based on Darcy's law for calculating the flow parallel to the fracture plate could be written as follows [21,22]: where u 2 is the Darcy velocity (m/s), K eq 2 is the effective hydraulic conductivity perpendicular to the fracture plate (m/s), h 2 is the pressure head (m), ∇ 2 is the gradient in the tangent direction (-), and notation q 2 is the source/sink term (1/s).

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.

Particle Tracking Algorithm for Advective Transport
The particle tracking algorithm is a typical approach to simulate the advective transport in complex fractured rocks [40,41]. This study assumes that the particles move with the flow of fractured rocks and have no chemical reactions and diffusion and dispersion. Therefore, the approach ignores the influences of sorption, degradation, and decay on transport. The movements of a particle are determined by the seepage velocity at a specific released location. With the predefined moving time steps, the particle locations at different time steps are calculated based on solving the following advection equation: where x is the particle location in space (m), t is the time (s), and U is the seepage velocity in a certain location (m/s). When the flow fields are available, the solution to the first-order differential equation can employ classical numerical solvers such as Euler and Runge-Kutta methods. In this study, the Euler method yields the following formula: where ∆t is the specified time step (s) for the particle tracking. The calculations of the particle traces and the traveling times rely on accumulating the recorded particle locations at different time steps. Figure 2 shows an example of the developed particle tracking algorithm. In this study, the basic computational meshes are triangle and tetrahedron. The calculation of particle movement follows the seepage velocities through triangular faces on each tetrahedral element. The intersection points between the trajectory line and element faces are calculated based on the Ray-Plane intersection test [42][43][44]. Figure 2a also shows a particle trace passing through several tetrahedral elements. Note that the velocities are available at the nodes of each element. This study uses the linear interpolation algorithm to obtain the velocities at the particle locations. This developed module can record the locations on element faces, the distances in elements, and the traveling time inside elements. The total travel length and time accumulate the step-by-step distances and times in different elements along the particle trace for each released particle. Figure 2. The demonstration of the particle tracking algorithm for a particle trace and the associated element along the trace line: (a) The overall particle moving paths starting from the point "Start" and end with the point "End"; (b) The local enlargement for the selected elements (marked with grey color in (a) and the particle trace outside the selected elements; (c) The drawing of local enlargement that shows the particle trace and the associated intersection point on the element face. Figure 2b,c are the local enlargements that zoom into two adjacent elements (marked with grey color) shown in Figure 2a. Figure 2b shows the particle trace outside the two selected elements, and the entrance point is marked on one face of the right element. Again, the point can be obtained based on the Ray-Plane intersection test. The 3D velocities at the marked location are then calculated based on the velocities at nodes of the element face. The traveling path inside the right element follows the trajectory of the velocity vectors determined at the point on the element face (see Figure 2c). This trajectory can reach another face of the right element, and the intersection point on the face can be determined (marked in Figure 2c). This face is shared with the left and right elements. Based on the velocities at the intersection point, the traveling distance and time are then determined for the right element. This study has implemented a process for the particle tracking algorithm.

Workflow for the Test Cases
There are a series of steps included in the workflow of this study (see Figure 3). These main tasks are developing conceptual models, generating grids or meshes for different model concepts, simulating steady-state flow, and conducting particle tracking for the advective transport. The quantitative comparisons of head distribution and particle traces are then based on the results obtained from the DFN, ECPM, and the developed HD models. In this study, the flow and advective transport solutions of the DFN and ECPM model were obtained based on the FracMan and DarcyTools software. Notice that the DFN model (i.e., FracMan) only focuses on the fractures [25,26], while the ECPM model (i.e., DarcyTools) considers the lumped behavior of the fractures and rock matrix [27][28][29]. In test Case I, the designed fractured rock is relatively simple. However, with the more complex feature of test Case II, the FracMan model (i.e., DFN) was excluded in this case because the feature to set fracture geometries and boundary conditions is not available in the DFN model. This study aims to assess the flow and advective transport based on the developed HD model. We used the solution from FracMan to obtain the DFN results and the solution from DarcyTools to obtain the ECPM results. The conceptual models for the test cases are similar to the used DFN and ECPM models. However, the boundary conditions and release of particles might be modified to fit the fractures in the DFN model and specify the cells in the ECPM model. Figure 4 shows the conceptual models for the test cases. In test Case I, the simple benchmark case includes three orthogonal fractures in the fractured rock (see Figure 4a). The domain size is a cuboid with 3 m × 4 m × 5 m in x-, y-and z-directions. There are three fracture plates in the simulation domain, i.e., F1 to F3 (F1: x = 1.5 m; F2: y = 2 m; F3: z = 2.5 m). All remaining zones of the domain represent the rock matrix, and the hydraulic properties are assumed to be homogenous and isotropic. A constant head condition was assigned to a rectangular area on the top plane, and a similar constant head boundary condition was specified on a rectangular area on the bottom of the simulation domain. The specified hydraulic heads of H in = 4 m and H out = 1 m are on the inlet and outlet areas, respectively (see Figure 4a). All remaining areas on the boundaries were assigned no-flow boundary conditions. In this study, assigning the in-and out-stresses on the matrix zones aims to evaluate the behavior of the flow that moves in and out of the interfaces between the fractures and matrix. Note that the DFN model focuses on the fractures only. However, the ECPM model considers the lumped behavior of the fractures and rock matrix. The assigned boundary conditions are either the fracture edges (i.e., DFN) or the cell faces (i.e., ECPM). Table 1 lists the hydraulic properties for the fractures and the rock matrix. Note that the apertures of the fractures were assumed to be constant for two test cases.   With the available flow field in Case I, there are 1000 particles randomly released in the inlet area. Note that the starting locations of particles and the particle movements are different in different models. In the DFN model, which only simulates the fractures, the particles are released on the fracture edges connecting to the inlet boundary area (see the marked solid pink lines on the top plane in Figure 4a). For the ECPM model, the particles are randomly released on the cell faces on the inlet boundary area (see Figure 4a). These two types of particle release scenarios are conducted in the HD model for comparison purposes. Figure 4b shows the conceptual model of the second benchmark case. In Case II, the domain size is relatively larger than that in Case I. A cylinder object represents a deposition hole (DH) of a disposal facility in the simulation domain, and two intersected fractures cut the DH. The domain size is 8 m × 8 m × 16 m in x-, y-and z-directions, and the cylinder DH has 8 m in height and 1.75 m in diameter. There are three fracture plates in the simulation domain, i.e., F4 to F6. The sizes for F4 and F5 are 6 m × 6 m and the area for the F6 is 5 m × 14 m. The apertures for the fractures are 0.1 m. Two intersected fractures cut the DH, i.e., F4 and F5 (see Figure 4b). In this case, the specified hydraulic stress is the same as that in Case I (i.e., the hydraulic head of H in = 4 m and H out = 1 m are assigned on the inlet and outlet areas, respectively). Note that three fractures have no contact with any domain boundaries, which means that the specified boundary conditions on the inlet and outlet areas need to be on the cell faces. Therefore, the DFN model was excluded in Case II because of the special setting of the conceptual model. The condition has shown the limited features of the DFN models that focus on the fractures only. Table 1 also lists the hydraulic properties and the convergence criteria for the simulation cases. In Case II, the DH hydraulic conductivity is the same as that in the rock matrix. For most practical problems, the hydraulic conductivity of DH should be lower than 1.0 × 10 −12 m/s to avoid a high corrosion rate and a possible release of radionuclides. The DH hydraulic conductivity value of 1.0 × 10 −10 m/s in this study could be a worse scenario for the DH design. The concept was proposed by studies of the case in Sweden [45]. In this case, we had made the ECPM model capture the flow and transport influenced by the fracture geometries. Therefore, the fracture aperture was assumed as 0.1 m to reduce the total computational grid in the ECPM model.

The Conceptual Models of the Test Cases
Based on the safety assessment report proposed by Sweden and Finland, i.e., the well-experienced countries that provide the general disposal concept for spent nuclear fuel, there are three potential pathways for releasing radionuclides if the canister in the DH is destroyed by corrosion or shear movement [6,46]. One of the potential and high probability pathways is the Q1 path, which is a concept for radionuclides released through the intersection point between a fracture and a DH. The intersection point usually has the highest velocity, one of the key input parameters for the corrosion calculation. Typical approaches consider the intersection points to release particles for evaluating the Q1 pathways. These released particles are in the highest velocity grid centers or nodes at the intersections to assess the possible pathways. In Case II, the particle release followed the concept of the Q1 path. In this study, the HD model could search for the highest velocity near the intersections between DH and fractures and release a particle to track the particle movement. It is interesting to obtain the intersection position, the highest velocity, and the overall transport behavior in different models. In addition, there were 48 particles released along the intersections between fractures and the DH to assess the variations of the particle movements in the case.

The Software and Computational Meshes and Grids for the Test Cases
In Case I, three fractures were connected to the domain top and bottom surfaces so that the DFN model could generate the mesh and directly assign the boundary conditions on the edge of the fractures. Thus, the 2D triangular meshes were created for fractures in the DFN model and no other mesh or gird for the remaining parts, i.e., the rock matrix. This study used the commercial software, FracMan, to simulate the groundwater flow in the DFN model. FracMan is the commercial software developed by Golder Associates for fracture generation and flow and transport simulations. The flow simulation package in the software is based on the finite element numerical scheme. Many previous investigations employed the software to simulate flow and transport in fractured rock systems (e.g., [25,26]).
The ECPM model uses the hexahedron as the computational grid. This study used DarcyTools to simulate the groundwater flow for the ECPM model. The DarcyTools is an integrated computational package that links several computational modules for the hydrogeological analysis of nuclear waste repository in fractured rocks. The integrated package is based on a workflow that allows combining the DFN and CPM models. The workflow first generates a fracture network or reads from a FAB file, a text file for recording the fracture positions and geometry in space. One of the modules then up-scales the fracture network to become the effective hydraulic properties of grid cells. The effective hydraulic properties, including the hydraulic conductivity and diffusivity, can represent the complex features of the fractured rocks. The detailed information for the computational package can be obtained in many previous investigations (e.g., [27][28][29]).
In this study, the proposed HD model uses the 2D triangular mesh for fractures and tetrahedron mesh for the 3D rock matrix. We developed a program to read the fractured file (the FAB file) from FracMan software and the DarcyTools computational package to ensure the compatibility of generated mesh for the fractures in the test cases. For the rock matrix in the test cases, the hydraulic properties were considered homogeneous and isotropic in the entire simulation domain. In the developed HD model, the 3D tetrahedron mesh was generated associated with the 2D fracture meshes. Therefore, the fractures and rock matrix could share the same nodes on the fractures and the matrix interfaces. The coupled solutions then rely on solving the system of equations developed based on all the nodes and elements in the fractures and matrix.
In this study, Case I is relatively simple as compared to Case II. We used the DFN model (i.e., FracMan) to generate the fracture geometry and the hydraulic parameters in Case I. In the DFN model, there was no mesh generated for the matrix. The generated nodes and elements for the DFN model are 6433 and 12,624, respectively. The averaged node distance for the fractures is 0.1 m. In the HD model, we had used the averaged node distance of 0.1 m in the matrix to connect triangular fracture mesh. Figure 5 shows the generated mesh for the HD model. Figure 5a is the exterior view of the domain composed of generated fracture and rock matrix elements. Figure 5b,d present the inner views of the generated mesh for Case I. The mesh for three intersected fractures is shown in Figure 5c. Based on the average node distance for 2D triangular mesh and 3D tetrahedron mesh, the number of nodes is 49,747, while the numbers of 2D triangular and 3D tetrahedron elements are 9147 and 290,324, respectively. In Case I, a similar cell size of 0.1 m was used in the ECPM model for the up-scaled hydraulic properties and flow and transport simulations. The total cell number of the ECPM model is 131,072, based on the cell numbers of 32, 64, and 64 in x-, y-, and z-directions. Figure 6 shows the generated mesh for Case II. Figure 6a shows the inner view of the test Case II. The DH object embedded in the simulation domain needs to be specified before the mesh generation. The surface of the DH object is the inner boundary for our mesh generation. Note that the fractures and the associated 2D fracture meshes are also the internal boundaries used to constrain the 3D mesh generation. In this study, the 2D fracture geometry is obtained based on the FAB file generated from the FracMan software. Figure 6b,c further zoom in the local region in the simulation domain and present the generated meshes for the 2D fractures and 3D matrix. In this case, the average node distance for 2D triangular mesh and 3D tetrahedron mesh in the HD model is approximately 0.125 m, and the number of nodes for the fractures and matrix is 412,174. Therefore, the numbers of 2D triangular and 3D tetrahedron elements are 16,734 and 2,511,187, respectively. This study used the same node distance of 0.125 m for the cell size specified in the ECPM model (i.e., DarcyTools). The number of total cells in the ECPM model was 524,288 based on 64 × 64 × 128 in x-, y-, and z-directions. Note that the mechanism for grid generation in DarcyTools follows a 221 rule that two adjacent cells always have a cross and longitudinal size ratio of 2:1 at most [27][28][29].

Results and Discussion
This section focuses on presenting the results based on the specified conditions defined in the benchmark cases. Case I used the steady-state flow fields from DFN, ECPM, and HD models, while Case II focused on comparing flow and advective transport obtained from the ECPM and HD models. Note that a DH was assigned in Case II to evaluate the dynamics of the particles induced by the fractures. Figure 7 shows the steady-state groundwater flow fields for Case I. Figure 7a,b compare results obtained from the HD and ECPM models, while Figure 7c,d present the head distributions only on the fractures. The results have shown the consistent solutions obtained from three different models. Specifically, the developed HD model enables the simulations of flow in fractures and rock matrix, which shows the advantage of the system of equations to be solved simultaneously. The concept proposed in the HD model can accurately capture the behavior of flow interactions between fractures and matrix (Figure 7a,b) and in the fractures network (Figure 7c,d). There is a difference of four orders of magnitude in hydraulic conductivity between the fractures and rock matrix in this case. The results of the HD model agree well with those obtained from the ECPM model in the entire domain region and the DFN model in the fracture network. Figure 8 presents the flow results based on the conditions specified in Case II. This study excluded the DFN model because the three fractures have no contact with the domain surfaces. The results in Figure 8 indicate that the solutions obtained from the HD model show slightly deviated from the ECPM model, especially near the corners away from the specified head boundary conditions. In general, the overall pattern of the flow was captured by the developed HD model (Figure 8a,b). The results in Figure 8c Figure 9 shows the particle traces and the associated trace lengths in Case I for three different models. Figure 9a,b show the release locations and the particle traces for the ECPM model. The results for the DFN model are shown in Figure 9e,f, while Figure 9c,d,g,h show the release location and the particle traces for the developed HD model. For comparison purposes, the HD results in Figure 9g,h only show the particles released along the edges of the fractures. In Case I, particles were randomly released on the inlet area, i.e., the specified head on the top boundary. The fracture hydraulic conductivity is four orders of magnitudes higher than the rock matrix. Figure 9a,b show that most particles directly penetrate the F3 and move to the outlet area with relatively small trace lengths in the ECPM model. However, the developed HD model shows that the fast pathway in F3 has significantly controlled the local flow and forced the particles to move along the fracture. No particle trace is shown in the block below the F3 fracture (see Figure 9c,d). All the particles move along the F3 and then follow the flow of the F1 and F2. In the ECPM model, the behavior of particle movement shows that the characteristic of high velocity in F3 is weak because the conductivity values of the F3 are lumped with the rock matrix conductivity in the upscaled cells in the ECPM model. Therefore, the hydraulic conductivity of the cells intersected by F3 might be slightly higher than the rock matrix after the upscaling process. In the ECPM model, the cell sizes intersected by fractures should be smaller than the fracture aperture to resolve the detailed flow behavior. However, it is impractical to consider small cell sizes for large-scale and complex fractured rock systems.   Figure 9g,h, the particle traces on the virtual fractures were extracted from the HD model for comparison purposes. There are special trace circulations obtained on the F3, which was recognized as the local flow effect. The local flow effect is a unique feature in the 3D DFN model because of the interplay between fracture geometry and boundary conditions for flow in the fracture plates [47]. The local flow effect has little influence on flow in a 3D DFN but may constitute a mechanism that contributions to the long breakthrough tails and retardation observed in transport. The results obtained from the developed HD model captured the behavior of the local flow effect. However, the local flow behavior is relatively small as compared with those obtained from the DFN model. Figure 9. The particle traces and lengths obtained from the ECPM, DFN, and HD models for Case I: (a,b) The particle traces for the ECPM model; (c,d) The particle traces for the HD model; (e,f) The particle traces for the DFN model; and (g,h) The particle trace along the fractures for the HD model. Table 2 summarizes the advective transport based on particle tracking applied to different models for Case I. Note that the velocity of each particle could be calculated by the observed trace length and the travel time. In Table 2, the results show that the mean trace length and the standard deviation (STD) and the coefficient of variation (CV) of the HD model are 1.14 times and 2.91 times higher than that of the ECPM model. There is no considerable deviation on the minimum trace length, but the maximum trace length of the HD model is 1.64 times higher than that of the ECPM model. The result indicates that most particles travel longer distances in the HD model because the fractures dominate the transport when the particles reach the fractures in the simulation domain. In this study, the STD of the trace length for the HD model is higher than that of the ECPM model. In general, the results show that the ECPM model underestimates the particle trace lengths and travel times in the case.

Simulations of Advective Transport
The particle movement statistics in Table 2 also show that the mean trace length and STD of the HD model are 0.98 and 0.59 times less than those of the DFN model. The minimal trace lengths for HD and DFN models are identical. However, the maximal trace length of the HD model is less than that of the DFN model. We found that the relatively significant local flow circulations in the DFN model lead to a longer trace length because of the particles that move in the horizontal fracture F3 [47]. Table 2 also shows the statistics of the traveling time for the released particles. Similar to the results of trace length, the HD model obtained lower mean and STD values of the particle travel time than those obtained from the ECPM model. Specifically, the minimum travel time of the HD model is two orders of magnitude less than that of the ECPM model, while the observed maximum travel time of the HD model is approximately half of the time obtained from the ECPM model. The results indicate that most particles in the HD model have fast traveling paths dominated by the fracture network. The short travel time is critical in conducting the safety assessment at a specific site. In general, the DFN model obtains relatively short travel times. The mean, minimum, and maximum travel times in the DFN model are generally one to two orders of magnitude smaller than those in the ECPM. However, the travel time variation (i.e., STD and CV) for the DFN model is relatively small compared with the results obtained from the HD and the ECPM models. This result is reasonable because the DFN model has constrained the particles to transport in the fracture network. Table 2 shows that the travel time statistics for the HD model vary between the results obtained from the ECPM and DFN models. The results in Table 2 also show that the mean velocity of the HD model is 7.62 times higher than that of the ECPM model. In addition, the STD of the HD model is 9.62 times higher than that of the ECPM model. Note that the velocity calculation is based on the trace length and the travel time for a specific particle. The minimum velocity for the ECPM and HD models are in the same order. However, the maximum velocity of the HD model is considerably higher than that of the ECPM model because the particles in the HD model might move in the fractures. For the results obtained from the fracture-based DFN model, the HD model could reproduce well the overall behavior of the advective transport in fractures. The HD model enables the particles to move in fractures and the rock matrix with the feature of interactions between fractures and the rock matrix. We had discussed the computational issues for different approaches. The ECPM model requires the cell sizes to be small to resolve the detailed flow dynamics near fractures. The DFN ignores the influence of the rock matrix, which might partially decrease the velocity of advective transport. In this study, the developed HD model considers the virtual fractures embedded on the element face could be a feasible approach for problems with practical scales and complexity.
In Case II, the objective is to evaluate the influence of the DH added in the simulation domain. Figure 10 shows the released location and the particle traces for the ECPM and HD model. Table 3 lists the detailed information for the results of particle tracking simulation for a specified release point in Case II. Note that the particle release point is the highest velocity based on the flow simulations obtained from the HD and ECPM models. Previous investigations had shown that the highest velocity was usually obtained along the intersections between the DH and fractures. In Case II, we had assigned a relatively large fracture aperture of 0.1 m. The large fracture aperture could reduce the number of cells for the ECPM model to capture the flow and transport influenced by the fracture geometries. However, this large fracture aperture might create differences in the highest velocities obtained from the ECPM and HD models. The location of a cell center in the ECPM model might not be consistent with the virtual fracture assigned on the element face. Figure 10. The particle release location and the particle trace for the ECPM and HD models in Case II. (a,b) The start location and particle trace for the ECPM model; (c,d) The start location and particle traces for the HD model. In most practical problems, the highest velocity is the critical parameter for safety assessment in the near field of a disposal facility. The velocity is one of the parameters for calculating the buffer erosion and canister corrosion rate inside the deposition hole. The highest velocity location is also the point for releasing the containment in most practical problems [5]. Therefore, the value and location of the highest velocity are essential for the corresponding safety assessment of the disposal facility. The highest velocity location in the ECPM model is (2.4375 m, 3.6875 m, 7.0625 m) in Figure 10a,b, while the highest velocity location in the HD model is (3.390 m, 3.573 m, 5.000 m) in Figure 10c,d. In addition, the highest velocity value of the HD model is 47.88 times higher than that of the ECPM model. In this study, the ECPM and HD models obtain similar head distribution for Case II. However, the detected highest velocity locations in the two models are considerably different (see Table 3). The differences in the detected locations reveal two important technical issues, including (1) The different grid or mesh generations and hydraulic property fields are critical in influencing the location of the highest velocity on the intersection between a fracture and the DH, and (2) The z coordinate of the highest velocity location of the HD model is 5 m which is exactly the z coordinate of the fracture F5, while the location of ECPM model is at the center of a grid, which is 7.0625 m in this case based on the applied computational grids. In this study, the developed HD model considers the intersection as the inner boundary for mesh generation. The precise location of the highest velocity can be identified based on the solution at the node on the intersection.
This study further evaluated the uncertainty introduced by the fracture and DH intersections introduced in Case II. There were 48 particles released on the two intersections cut by fractures F4 and F5 (see Figure 10). Figure 11 shows the particle traces of the 48 released particles, and Table 4 lists the statistics based on the collected particle traces and the travel times. The results in Figure 11a show that the particle traces in the ECPM model are similar, and the trace lines are uniformly distributed on the fracture surfaces. In the ECPM model, the pattern of the particle traces could capture well the two horizontal fracture geometry (see Figure 11a). Figure 11b shows the pattern of the particle traces in the HD model. The result demonstrates that the particle traces have been constrained in the narrow areas on the two fractures. Parts of the particles released from the intersection between F4 and DH migrate directly along the fracture F6 and then migrate to the rock matrix and reach the outlet boundary. A small portion of the particles might move along the F5 and then follow the particle traces starting from the intersection between DH and the fracture F5. The results present the significant difference between the ECPM and the HD models for the local flow behavior near the DH and fractures. The evaluations of the particle trace and traveling time can vary considerably if the release locations of the particles are different.  The statistics for the advection transport are summarized in Table 4. In general, the particle travel times obtained from the HD model are relatively small, leading to a relatively high transport velocity ( Table 4). The mean velocity of the HD model could be two times higher than that in the ECPM model. Table 4 also shows that the HD model obtains relatively low STD values for the trace lengths and travel times but obtains high variation in transport velocity. Note that the HD model has considered virtual fractures in the simulation procedures. The particle release locations are precisely defined in the HD model. The fractures in the HD model act as a high-velocity pathway for groundwater flow. The particle migrations could be in fractures or rock matrices, depending on the local flow dynamics. For most practical problems, capturing the complex flow dynamics near DHs and fractures is essential for accurate evaluations of the large-scale and long-term transport processes.

Conclusions
This study has presented the concept of the HD approach for the simulation of advective transport in fractured rocks. Specifically, the developed HD model uses the 2D triangular mesh for fractures and tetrahedral mesh for the 3D rock matrix in a simulation domain and allows the system of equations to be solved simultaneously. The study employed two synthetic cases to assess the developed HD model. Solutions of flow and advective transport were compared with those obtained from the DFN model and ECPM model. In this study, the advective transport was conducted based on the particle tracking algorithm.
The simulation results show that the HD model is flexible in considering the concepts of DFN, ECPM, or both. Based on the test cases, the flow fields obtained from the HD model agree well with those obtained from the ECPM and DFN models. However, the local flow dynamics led to significant variations in the advective transport. Fractures in the simulation domain play an important role in controlling the local flow patterns. In addition, the released particle locations could also influence the particle traces and the travel times. The ECPM model might need fine grids to resolve the interaction dynamics between the fractures and rock matrix. The DFN model focuses on the fractures only and has the limited feature to account for the interaction between fractures and rock matrix. The local flow circulations are the special characteristic obtained from the DFN and the developed HD model.
The transport statistics relied on collecting the particle traces and travel times obtained from different models. The results in Case I indicated that the HD and DFN models obtained similar particle traces and travel times because the particles migrations have been restricted in the fractures. The particle released in the matrix and fractures might obtain different results from the ECPM and HD models. For the case with the particle release in the matrix (i.e., Case I), the HD model considered the detailed fracture flow and obtained a longer mean trace length than that obtained from the ECPM model. In Case II, the particles were released on the intersections between the DH and fractures. The HD model obtained a relatively small mean trace length and travel time. The identified highest velocity in the HD model is significantly greater than that in the ECPM model. The results have shown the feasibility of the HD model for accurate simulations of flow and advective transport in complex fractured rock systems.