This paper presents a numerical solution to the problem of a transversely isotropic multilayered viscoelastic porous rock foundation under vertical circular loadings. A fractional Poyting-Thomson model is proposed to elaborate the rheological characteristics of the rock mass. The viscoelastic solution of the rock foundation in the transformed domain is deduced by using the elasticity-viscoelasticity correspondence principle. The actual solution is obtained by the integral inverse transform. The reliability and correctness of the fractional viscoelastic model are verified with the experimental data. Numerical examples are conducted to further verify the accuracy and precision of the presented theory. The influence of viscoelastic and transversely isotropic parameters and the rock mass stratification on the time-dependent behavior of a rock foundation is further explored.