Abstract: The theory of fractional calculus is introduced into the Kelvin-Voigt constitutive model to describe the mechanical behavior of viscoelastic saturated soils. Applying the Laplace transform upon one-dimensional consolidation equation of saturated soils and the fractional order derivative Kelvin-Voigt constitutive equation, the analytical solutions of the effective stress and the settlement are derived in the Laplace domain. Then the semi-analytical solutions to one-dimensional consolidation problem under arbitrary loadings in physical space are obtained after implementing the Laplace numerical inverse transform using the Crump’s method. As the case of viscoelasticity, the simplified semi-analytical solutions under exponential loading in this study are the same as the available analytical solutions in literatures. It is indicated that the proposed solutions under arbitrary loading are reliable. Finally, parametric studies are conducted to analyze the effects of the related parameters on the consolidation settlement. The results show that the process of one-dimensional consolidation of viscoelastic saturated soils with fractional order derivative is related to viscosity coefficient and fractional order. The larger the fractional order is, the more quickly the consolidation settlement occurs; and the higher the viscosity coefficient is, the slower the consolidation settlement takes place. The trend of loadings is consistent with the variation pattern of soil settlement caused by the change of the load parameters, and the final settlement is identical. The present study can be of help to further understand the consolidation behavior of viscoelastic saturated soils.