Fractured media have been extensively studied due to their significant impact on hydraulic properties of subsurface systems. However, accurately modeling and simulating flow and transport in complex 3D fractured media still remains a challenge. For this purpose, this work aims to develop a new 3D modeling approach and investigate flow and mass transport in fractured media. First, a modeling approach is presented for generating 3D complex structures and high quality finite element grids. This approach allows for the presence of a range of complex geometries, including curved fracture, multiple fractures and the coexistence of fracture‐inclusions. The resulting models and grids are of high quality, enabling accurate simulation and analysis of complex systems. Furthermore, stress‐dependent aperture and roughness are integrated in the computational scheme using Barton–Bandis law. The Galerkin finite element method is used for numerical discretization, then pressure and concentration are calculated in a unified formulation. Later, numerical tests are performed to analyze the role of 3D fractures in flow and mass transport. Simulation results show that fractures have a larger impact on pressure distribution and flow rate compared to inclusions. Fracture roughness increases the aperture then the effective conductivity is improved as well. Curved fractures may reduce flow rate more than planar fractures due to the tortuosity effect. The variation of equivalent permeability tensor is analyzed with varying orientation, roughness, and number of fractures. The impact of stress and fracture roughness causes a variation of fracture permeability and a high pressure gradient on the fracture‐matrix interface, resulting in different patterns of concentration evolution. Overall, this investigation provides insights into the role of 3D fractures in subsurface flow and transport, with implications for various applications.