The clogging-unclogging process in gap-graded soils is a result of the migration of seepage-driven fines, which subsequently induces measurable changes in the local hydraulic gradients. This process can be temporally observed in the variations of Darcy's hydraulic conductivity (K). The current study proposes an integrated statistical Monte Carlo approach combining the discrete element method and 2D computational fluid dynamics simulations to estimate the flow-dependent constriction size distribution (CSD) for a gap-graded soil. The computational inferences were supported with experimental results using an internally stable soil, which was subjected to one-dimensional flow stimulating desired hydraulic loadings: a hydraulic gradient lower than the critical gradient applied as a multi-staged loading pattern. The 35th percentile size of the flow-dependent CSD (Dc35) for both internally stable and unstable gap-graded soils becomes approximately equal to Dc35 at steady-state. However, a greater variation of larger constrictions persists for the unstable soils. This pilot study has shown the applicability of the proposed method to estimate flow-dependent CSD for a wide range of experimentally observed K values.