This paper studies the time-dependent behavior of a pile drilled in layered saturated viscoelastic rock-soil mass due to a vertical load. By virtue of the finite element method (FEM), the pile is discretized into three-node axial bar elements. On the basis of the consolidation solution of layered transversely isotropic saturated rock-soil mass, the fractional Poyting-Thomson model and fractional Merchant model are used to simulate the rheological properties of rock and soil, respectively. The viscoelastic solution of layered saturated rock-soil mass under annular linear loads is derived according to the elastic-viscoelastic correspondence principle. Taking the above solutions as the kernel functions of the boundary element method (BEM), and combining with the stiffness matrix equation of a pile, a coupling formula of pile-soil-rock interaction is further established by the FEM-BEM coupling method. A MATLAB code is developed for numerical calculation, and the correctness of the present theory and calculation method are verified by comparing with the existing solutions, field tests and ABAQUS analyses. Finally, the key factors influencing the time-dependent behavior of piles are in-depth discussed.