This paper presents a numerical method for the advection-diffusion process in fractured porous media described by a kind of hybrid–dimensional fracture models, in which fractures are represented as $(n-1)$–dimensional interfaces embedded within the $n$–dimensional porous matrix. We develop a numerical approach that combines the discontinuous Galerkin (DG) method for spatial discretization with the backward Euler method for temporal discretization. Leveraging the advantage of the DG method in flexibly handling boundary conditions, the DG–backward Euler method is extended from the model with a single fracture to one with intersecting fractures. Theoretical analysis shows that the optimal order error estimate for the concentration in the discrete $H^1$–norm is obtained. Numerical experiments with a single fracture and intersecting fractures are all carried out to verify the accuracy of the theoretical results.