Seepage plays a crucial role in the mechanical behavior and damage modes of geotechnical materials. In this work, based on the unsteady seepage equation, a hydraulic coupling numerical simulation algorithm combining interpolation finite difference method (FDM) and discrete element method (DEM) is proposed to explore the intrinsic mechanism of the interaction between geotechnical materials and the seepage process. The method involves constructing an irregular fluid calculation grid around each particle and deriving the two-dimensional unsteady seepage governing equation and its stability conditions using interpolation and the FDM. The efficiency of the seepage calculation was investigated by numerically varying the parameters of the difference format. The method was applied to simulate the generation of gushing soil in a sinking area of a sunk shaft under hydraulic drive conditions. The results indicate that the improved FDM can effectively simulate the two-dimensional seepage of soil with high calculation efficiency. The hydraulic conductivity and time step positively correlate with the calculation efficiency of the difference format, whereas the spatial step has a negative correlation. The proposed method also accurately reflects the process of gushing soil damage. These results provide a solid theoretical basis to study the geotechnical seepage field and its associated damage mechanisms.