The equivalent permeability of layered fractured rocks plays an important role in hydrocarbon recovery, underground energy storage, waste disposal management, groundwater hydrology, and subsurface contaminant transport. Borehole data contain some uncertainties/sampling bias during collection and interpretation. This sampling effect may lead to an inaccurate characterization of the fractured media. Studies show that long fractures with high permeability, which are rarely seen in borehole images, dominate the flow pattern and affect the overall permeability of the fractured system. This means that samples taken at any scale smaller than the scale of interest result in imprecise permeability upscaling. To understand this sampling problem, we have established an efficient sampling method to study the existence of the representative elementary volume (REV) in naturally fractured rocks. We selected a collection of outcrop data represented by discrete fracture and matrix (DFM) models in which fracture apertures are mechanically constrained. A finite-element-finite-volume approach is utilized to characterize the flow behavior of DFM models. Multiscale random sampling is combined with flow-based upscaling to determine the equivalent permeability tensor and its anisotropy by considering the variable orientation of the stress state and fracture roughness. Our findings indicate a convergence towards a scale-invariant equivalent permeability and fracture density with increasing sample size, and the equivalent perme- ability itself has a multimodal distribution. The spatial variation of the permeability tensor and the change in the degree of anisotropy with sample size reflect the inhomogeneity of the fracture patterns.