Wave propagation in partially saturated porous continua is a fundamnetal subject in Civil Engineering, Petroleum Engineering, Bioengineering, Earthquake Engineering, and Geophysics, etc. For such problems, there exist different types of theories can be used to model as Biot’s theory, 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. With the material data of Massillon sandstone, 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 test results. Based on the Convolution Quadrature method, the time domain results referring to the displacement and the pore pressure are obtained and compared with the corresponding results of saturated poroelasticity. With the proposed governing equations, the fundamental solution is deduced following the Hörmander’s method. The boundary integral equations are established based on the weighted residual method. The regularized fundamental solution is then implemented with the help of the open source C++ BEM library HyENA. The latter numerical examples aim to validate the proposed method, and some applications are also presented.