The long-term seepage characteristics of fluid flow in saturated soils are of great importance. The Riemann-Liouville (R-L) fractional derivative was adopted to modify classical Darcy’s law (hereinafter referred to as modified seepage model) to describe the evolution of soil permeability during the long-term seepage process. Data fitting of experimental results given in published literature shows that the modified seepage model could more accurately describe the nonlinear evolution of fluid velocity with time. Moreover, the anomalous permeability coefficient value obtained with the modified seepage model is found to be reasonable. The R-L fractional diffusion equation was derived by integrating the modified seepage model into the one-dimensional Biot consolidation model. The explicit (time domain)-implicit (space domain) difference method was employed to discretize the above equation, and the correctness of the algorithm was verified through two numerical examples. On this basis, the influence of the modified seepage model parameters on the one-dimensional consolidation process of saturated soils was investigated. The results show that the fractional order β reflects the degree of soil permeability decay. The higher the fractional order, the lower the soil permeability, which leads to a further decrease in the consolidation rate; additionally, the permeability coefficient k plays a more dominant role than kβ. Thus, k should not be assumed to be zero to simplify the fitting process.