A model consisting of three-dimensional (3D) layered saturated ground and absorbing boundary applied at the bottom was established in this paper. By introducing the Fourier and orthogonal transformations

Biot’s equations for the saturated soil are generated into equations governing Love and Rayleigh modes. By dividing the ground into several thin layer elements along the vertical direction and introducing linear interpolation function for the elemental displacement fields

the equations governing Love and Rayleigh modes for a single thin layer element are derived and then assembled to obtain global modal equations. After solving the homogenous global modal equations

the eigenvalues and corresponding eigenvectors of the saturated ground are obtained for Love and Rayleigh modes

separately. Eventually

the Green function for the 3D saturated ground is obtained for point loads by mode superposition method. Numerical results are presented to portray the evolution processes of eigenvalues of Love and Rayleigh modes with respect to the increasing excitation frequency and soil permeability. The Lamb problem of a homogeneous and layered saturated half-space is solved and compared with the numerical integration solutions given by Philippacopoulos. Through the comparisons

the obtained Green function is proven to be valid and accurate. Meanwhile

it does not require numerical integration with respect to Hankel functions and can take the soil stratification into consideration in a natural way

so as to improve the computation efficiency.