Abstract We consider a 3D-2D-1D mixed-dimensional diffusive model in a fractured porous medium coupling the 1D model along the centerline skeleton of a tubular network, the 2D model on a network of planar fractures and the 3D model in the surrounding matrix domain. The transmission conditions are based on a potential continuity assumption at matrix fracture interfaces, and on Robin type conditions at the resolved interfaces between the tubular network and the matrix and fracture network domains. The discretization of this mixed-dimensional model is formulated in the gradient discretization framework (Droniou, J., Eymard, R. & Herbin, R. (2016) Gradient schemes: generic tools for the numerical analysis of diffusion equations. ESAIM Math. Model. Numer. Anal., 50, 749–781), which covers a large class of conforming and nonconforming schemes and provides stability and error estimates based on general coercivity, consistency and limit-conformity properties. As an example of discretization fitting this framework, the mixed-dimensional version of the vertex approximate gradient (VAG) scheme is developed. It is designed to allow nonconforming meshes at the interface between the 1D and 3D-2D domains, to provide a conservative formulation with local flux expressions and to be asymptotic preserving in the limit of high transfer coefficients. Numerical experiments are provided on analytical solutions for simplified geometries, which confirm the theoretical results. Using its equivalent finite volume formulation, the VAG discretization is extended to a drying mixed-dimensional model coupling the Richards equation in a fractured porous medium and the convection diffusion of the vapor molar fraction along the 1D domain. It is applied to simulate the drying process between an operating tunnel and a radioactive waste storage rock with explicit representation of the fractures in the excavated damaged zone.