We present our new multiscale pairwise-force smoothed particle hydrodynamics (PF-SPH) model for the characterization of flow in fractured porous media. The fully coupled multiscale PF-SPH model is able to simulate flow dynamics in a porous and permeable matrix and in adjacent fractures. Porous medium flow is governed by the volume-effective Richards equation, while the flow in fractures is governed by the Navier–Stokes equation. Flow from a fracture to the porous matrix is modeled by an efficient particle removal algorithm and a virtual water redistribution formulation to enforce mass and momentum conservation. The model is validated by (1) comparison to a finite element model (FEM) COMSOL for Richards-based flow dynamics in a partially saturated medium and (2) laboratory experiments to cover more complex cases of free-surface flow dynamics and imbibition into the porous matrix. For the laboratory experiments, Seeberger sandstone is used because of its well-known homogeneous pore space properties. The saturated hydraulic conductivity of the permeable matrix is estimated from a pore size and grain size distribution analysis. The developed PF-SPH model shows good correlation with the COMSOL model and all types of laboratory experiments. We employ the proposed model to study preferential flow dynamics for different infiltration rates. Here, flow in fracture is associated with the term “preferential flow”, providing rapid water transmission, while flow within the adjacent porous matrix enables only slow and diffuse water transmission. Depending on the infiltration rate and water inlet location, two cases can be distinguished: (1) immediate preferential/fracture flow or (2) delayed preferential flow. In the latter case, water accumulates at the surface first (ponding), then the fracture rapidly transmits water to the bottom system outlet. For the immediate fracture flow response, ponding only occurs once the fracture is fully saturated with water. In all cases, preferential flow is much more rapid than diffuse flow even under saturated porous medium conditions. Furthermore, infiltration dynamics in rough fractures adjacent to an impermeable or permeable matrix for different infiltration rates are studied as well. The simulation results show a significant lag in arrival times for small infiltration rates when a permeable porous matrix is employed, rather than an impermeable one. For higher infiltration rates, water rapidly flows through the fracture to the system outlet without any significant delay in arrival times even in the presence of the permeable matrix. The analysis of the amount of water stored in permeable fracture walls and in a fracture void space shows that for small infiltration rates, most of the injected water is retarded within the porous matrix. Flow velocity is higher for large infiltration rates, such that most of the water flows rapidly to the bottom of the fracture with very little influence of matrix imbibition process