Many natural hazards involve unsaturated soils, such as rainfall-induced landslides, seepage-induced displacements, instability of dams and levees, as well as, embankment collapses during wetting. Analysing these phenomena with a dynamic approach that considers large deformations and displacements is interesting for effective risk assessment and management, particularly in predicting post-failure responses. Among different numerical alternatives available to simulate large deformations, the material point method (MPM) [1,2] has been successfully developed over the past decades to address multiphase geotechnical problems, taking into account the soil skeleton and pore fluids (liquid and gas). In MPM-based modelling, the complexities associated with boundary conditions of infiltration and potential seepage in unsaturated soils increase because possible large deformations imply an a priori unknown position of the boundaries where to apply the boundary conditions. This study focuses on the treatment and implementation of hydraulic boundary conditions within the open-source code Anura3D [3]. The unsaturated soil is modelled using the MPM “2 phase–1 point” formulation with suction effect [4]. This is a simplified approach in which the gas pressure is assumed constant and the degree of saturation, dependent on suction through the water retention curve, is considered in the liquid phase mass balance equation. After briefly presenting the basis of the coupled hydro-mechanical formulation of MPM, the numerical implementation of the hydraulic boundary conditions is described and validated by simulating academic cases of slope instabilities. Finally, the methodology proposed is compared with alternative approaches available in the literature.