Multiphase geophysical mass flows, such as debris, mud or slush flows, pose a significant risk to people and infrastructure in alpine regions. The interaction between solid and liquid constituents is essential for understanding these flows, as it governs key phenomena like liquefaction of saturated sediments under rapid undrained loading, pore-pressure diffusion and consolidation. Yet, current operational models do not directly account for a solid-liquid coupling and are often based on a depth averaged approach, lacking a full three-dimensional resolution. In recent years, the multiphase extension of the Material Point Method (MPM) has proven to be a suitable candidate for modeling large-deformation multiphase dynamics. However, most contributions fall into two categories: (1) contributions from the graphics community, where visually appealing results often lack a physical foundation, and (2) advanced modeling efforts that primarily focus on simple verification tests. Additionally, many studies employ a single-point approach which cannot capture free liquid during phase separation, while the majority also assume fully saturated conditions, disregarding suction effects and the dependence of permeability and effective compressibility on soil saturation. In this study, we develop a three-dimensional two-phase two-point semi-implicit unsaturated MPM. This model strikes a balance between robustness against pressure oscillations and computational efficiency that enables slope-scale modeling of multiphase mass movement events. We assess the stability of the scheme and evaluate its performance in complex scenarios, extending beyond the standard verification tests such as the well-known one-dimensional consolidation. We validate the ability of the model to reproduce Reynolds dilatancy using the Modified Cam Clay model for soil plasticity and investigate the role of liquefaction and saturation degree in debris flow dynamics and bed erosion.