The modeling of porous media interacting with incompressible fluids is of great importance for a wide range of problems in physics and engineering. In natural systems, such coupled behavior commonly arises in granular materials whose pore spaces are partially or fully filled with fluid. The presence of interstitial fluid can fundamentally alter the mechanical response of these materials by reducing effective stresses and thus shear strength, leading to higher mobility through liquefaction.However, existing models often struggle to capture this two-phase interaction while ensuring numerical stability and independence from the chosen permeability under large-scale, large-deformation conditions. Here we present a stabilized fractional-step Material Point Method (MPM) that enables robust and efficient simulations of large deformation hydro-mechanical problems without compromising model accuracy due to stabilization effects. We systematically analyze the performance of both non-incremental and incremental fractional-step formulations, combined with various temporal strategies for enforcing inter-phase drag. The results demonstrate that the stabilized incremental formulation is the superior choice, without the need for more complicated drag schemes. This advancement enables more accurate large-scale simulations of deformation processes that critically depend on two-phase dynamics, including geophysical flows such as avalanches of rock, debris, snow, and ice, which frequently involve significant amounts of liquid water. Ultimately, this study establishes a foundation for further developments in hydromechanical modeling, paving the way toward incorporating additional physical mechanisms, such as unsaturated porous behavior or phase transitions associated with melting ice and snow.