The simulation of seismic wave propagation in fluid-filled fractured porous media poses a formidable challenge, requiring the resolution of a hydromechanical coupling problem, described by Biot's equations, in heterogenous media characterized by the presence of complex networks of thin heterogeneities (fractures). Key challenges in this endeavor include the generation of meshes to accurately represent realistic fracture networks, a significant hurdle in itself. Additionally, the low regularity of solutions impacts the convergence of discretization methods and high-performance solution techniques, such as geometric multigrid. In \cite{favino2020fully}, we proposed an automatic strategy to automatically generated meshes for fractured media. The idea is to start from a uniform mesh and a fracture distribution and refine the elements which have a non-empty overlap with at least one fracture. Iterating this process, this strategy produces a mesh which does not resolve the boundary of the fractures but it allows to approximate it with a desired accuracy. However, the resulting mesh, and hence, the linear system associated to the discretization of these meshes is extremely large. Another common approach to model fractures involves using hybrid-dimensional models, where fractures are downscaled to objects of a lower dimension and separate problems are solved in the background and fracture domains with suitable coupling conditions~\cite{flwmich,BaVe06}. The goal of this talk is twofold. On the one hand, we study application of the an algebraic multigrid method for the solution of the flow problem in fractured porous media. We show that the convergence rate of the method slightly depends on the mesh refinement and the material properties in the fractures. On the other hand, we propose a novel formulation for modeling fractures, which considers fractures as cuts in the domain and models their effects using Wentzell interface conditions~\cite{Ci}. The resulting model is elegant and straightforward, consisting in a primal formulation where the only unknown is the pressure in the background domain. The model can naturally handle fractures that cross each other. We analyze the discretization error for different orders of finite element and material properties, and we compare the arising conditions with standard analytical models used for fractured poroelastic media~\cite{Li88}.