## 1 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 cause the copper canister to corrode. Once the canisters are destroyed by corrosion, the radionuclides will release from the deposition holes through the fractures and migrate to the human environment. Therefore, the flow and transport dynamics near the disposal facility might considerably influence the evaluations of long-term transport behavior. Continuous and discrete are two common approaches to simulate the flow and advective transport in fractured rocks. The models of the approaches are either the discrete fracture network (DFN) [7; 8; 9; 10; 11] or the equivalent continuum porous medium (ECPM) [12; 13; 14], depending on the site-specific conditions and the interesting issues. The DFN model generates a series of individual fractures based on the statistical distribution of geometrical properties (i.e., the fracture orientation, size, and center position) and hydraulic parameters (i.e., transmissivity and hydraulic aperture) to represent the complex fracture system in fractured rock. The DFN model enables the simulations of detailed flow and transport dynamics within fractures and neglects the contribution of flow and transport controlled by the rock matrix [15]. The DFN model is useful for addressing the problems with small-scale domain, low fracture number or fracture intensity, and high discrepancy of hydraulic conductivity between fractures and the rock matrix. However, the interactions between the fractures and rock matrix might be important for contaminant transport. The sorption and diffusive process in the rock matrix could lead to plume migrations with long-tail behavior in fractured rocks. Many previous studies have investigated the inconsistency between the DFN model results and the site-specific observations [7; 8; 9; 10; 11]. The key advantage of the ECPM model is the efficiency in obtaining solutions for the equivalent flow and transport processes in large-scale fractured rock systems. The numerical cells in the ECPM model are considered the representative elementary volume (REV), where the lumped flow and transport properties are treated as equivalent values in the numerical cells. In the ECPM models, the hydraulic and transport properties in a REV cell are typically preprocessed or up-scaled to be a tensor that represents the anisotropic characteristic. Thus, the ECPM models have the same features as those available in most classical porous media models. Many ECPM models have been developed in the past few decades based on the coupled thermal, hydrologic, mechanical, and chemical (THMC) processes.