The analytical layer-elements for a single poroelastic soil layer and the underlying half-space are established using an algebraic manipulation and Hankel transform. According to the boundary conditions and adjacent continuity conditions of general stresses and displacements, a global matrix equation in the transform domain for multilayered saturated soil media is assembled and solved. Solutions in the frequency domain can be further obtained with an inverse Hankel transform. Numerical examples are used to examine accuracy of the present method and demonstrate effects of soil parameters and load conditions on dynamic responses of the multilayered poroelastic saturated soils. The analytical layer-elements for a single poroelastic soil layer and the underlying half-space are established using an algebraic manipulation and Hankel transform. According to the boundary conditions and adjacent continuity conditions of general stresses and displacements, a global matrix equation in the transform domain for multilayered saturated soil media is assembled and solved. Solutions in the frequency domain can be further obtained with an inverse Hankel transform. Numerical examples are used to examine accuracy of the present method and demonstrate effects of soil parameters and load conditions on dynamic responses of the multilayered poroelastic saturated soils.