This paper presents an advanced formulation of the time-domain two-dimensional (2D) boundary element method (BEM) for an elastic, homogeneous unsaturated soil subjected to dynamic loadings. Unlike the usual time-domain BEM, the present formulation applies a convolution quadrature which requires only the Laplace-domain instead of the time-domain fundamental solutions. The coupled equations governing the dynamic behavior of unsaturated soils ignoring contributions of the inertia effects of the fluids (water and air) are derived based on the poromechanics theory within the framework of a suction-based mathematical model. In this formulation, the solid skeleton displacements [Formula: see text], water pressure [Formula: see text] and air pressure [Formula: see text] are presumed to be independent variables. The fundamental solutions in Laplace transformed-domain for such a dynamic [Formula: see text] theory have been obtained previously by authors. Then, the BE formulation in time is derived after regularization by partial integrations and time and spatial discretizations. Thereafter, the BE formulation is implemented in a 2D boundary element code (PORO-BEM) for the numerical solution. To verify the accuracy of this implementation, the displacement response obtained by the boundary element formulation is verified by comparison with the elastodynamics problem.