Wave propagation in partially saturated porous continua is an interesting subject in Civil Engineering, Petroleum Engineering, Bioengineering, Earthquake Engineering, and Geophysics, etc. For such problems, there exist different theories, eg, an extension of Biot’s theory, the Theory of Porous Media and the Mixture theory. Based on the Mixture theory, a dynamic three–phase model for partially saturated poroelasticity is established as well as the corresponding governing equations in Laplace domain. This model is applied to a one dimensional column and the related analytical solution in Laplace domain is deduced. The three different compressional waves, the fast wave, the second, and the third slow waves are calculated and validated with the Biot–Gassmann prediction and Murphy’s experimental results. The time domain results are obtained with the convolution quadrature method. Within the limit of a saturation nearly to one the results are as well compared with the corresponding results of saturated poroelasticity. For the three-dimensional governing equations, the fundamental solutions are deduced following Hörmander’s method. The boundary integral equations are established based on the weighted residual method. After regularization, spatial discretization, and the time discretization with the convolution quadrature method the boundary element formulation in time domain for partial saturated media is obtained. The implementation is done with the help of the open source C++ BEM library HyENA. Finally, the code is validated with the analytical one-dimensional solutions of the column. Two half-space applications are as well presented.