In this paper, a predictive model for simulating temporal behaviors of clayey sand during seepage-induced erosion has been developed by coupling discrete element method (DEM) with computational fluid dynamics (CFD). In this model, the particle–fluid coupling simulation is solved by a “fixed coarse-grid” scheme in 3D particle flow code (PFC3D), and the suffusion of clay matrix in the initiation of erosion is converted to a degradation process of bonding strength between particles according to a degradation law. The law is derived from the well-known shear stress threshold law dealing with soil internal erosion, which is based on two erosion parameters—the critical shear stress and the erosion coefficient. Then the degradation law is implemented in the CFD–DEM model via developing customized code using the Python language. The ability of the model to predict the interfacial erosion of soils is confirmed by two numerical tests. The results are seen to match the empirical criteria, such as revealing a clearly defined critical tangential shear stress, beyond which erosion occurs, and a positive correlation between the rate of erosion and the pressure gradient. It is believed that the numerical model is able to reproduce the time-dependent evolution process of seepage-induced erosion in clayey sand.