The local flow field and seepage induced drag obtained from Pore Network Models (PNM) is compared to Immersed Boundary Method (IBM) simulations, for a range of linear graded and bimodal samples. PNM were generated using a weighted Delaunay Tessellation (DT), along with the Modified Delaunay Tessellation (MDT) which considers the merging of tetrahedral Delaunay cells. Two local conductivity models are compared in simulating fluid flow in the PNM. The local pressure field was very accurately captured, while the local flux (flow rate) exhibited more scatter and sensitivity to the choice of the local conductance model. PNM based on the MDT clearly provided a better correlation with the IBM. There was close similarity in the network shortest paths, indicating that the PNM captures dominant flow channels. Comparison of streamline profiles demonstrated that local pressure drops coincided with the pore constrictions. A rigorous validation was undertaken for the drag force calculated from the PNM by comparing with analytical solutions for ordered array of spheres. This method was subsequently applied to all samples, and the calculated force was compared with the IBM data. Linear graded samples were able to calculate the force with reasonable accuracy, while the bimodal samples exhibited slightly more scatter.