Abstract A three-dimensional boundary element model is presented for the computation of time-harmonic dynamic stiffness coefficients of piles and pile groups embedded in two-phase poroelastic soils. Piles are modelled as continuum elastic solids and soil as a fluid-filled poroelastic half-space governed by Biot’s theory. Piles are bonded to the surrounding medium along the contact surface where rigorous pile–soil interaction conditions are imposed. The technique has been applied to the calculation of vertical and horizontal dynamic impedance functions of piles and 2 × 2 pile groups. Selected numerical results are presented to portray the influence of excitation frequency, pile stiffness, porous medium properties and kind of pile–soil contact conditions.