Seepage-induced fines migration, referred to as suffusion, conceivably degrades the hydraulic performance of soil masses subjected to water flow. Suffusion susceptibility of the soil (i.e. internal instability) is empirically evaluated using the grain size distribution (GSD) with little emphasis on the permeation process of fines (i.e. mobile fines) in the matrix. Upon experiencing different hydraulic gradient histories, the mobile fines can detach (unclog) or filtrate (clog) in the pore space of the load bearing soil skeleton (i.e. primary fabric). Moreover, clogging and unclogging phenomena are instrumental to temporally altering the hydraulic responses, mainly the hydraulic conductivity (K) and critical hydraulic gradient required to initiate internal erosion of the soil. Therefore, real-time predictions of the hydraulic response of the soil can reveal the initiation and progression of suffusion. This study attempts to forecast the hydraulic response of the soil by quantifying the statistical distribution of K values under mobile fines contents (fc). To this end, the discrete element method (DEM) has been employed to estimate the void ratio distribution of the primary fabric. The statistical Monte Carlo simulations have been performed using the Kozeny-Carman equation and the void ratio distribution to estimate the statistical distribution of K values. A set of pilot-scale experiments were conducted on internally stable gap-graded soil subjected to stepwise hydraulic loading histories to allow validation for the numerical results. Based on the measured pore pressure distribution parallel to the seepage path, the possible clogging and filtration of fines were inferred, and the K variations were estimated. The change in GSD of the soil after testing was used as a measure for the mobile fines fraction. An average of 6–8% mass loss (from the fine fraction) was observed over a 45-min duration of the controlled-head flows. The experimentally observed K (1–3 × 10−4 m/s) and fc ranges fall well within the non-parametric bounds of the predicted K and fc (<10%) ranges.