In this work, we present a finite element formulation for variably-saturated porous geomaterials undergoing elastoplastic deformations. The deforming body is treated as a multiphase continuum, and the governing mass and momentum balance equations are solved in a fully-coupled manner. It is well-known, however, that mixed formulations of the type examined here may lead to unstable approximations unless the spaces chosen for the pressure and displacement interpolation satisfy stringent stability restrictions. Failure to choose a stable pair typically leads to spurious pressure oscillations and poor convergence behavior. Unfortunately, many seemingly natural combinations---including equal-order interpolation for all field variables---do not satisfy the necessary requirements. In this work, we propose a stabilized formulation, based on a minor modification of the variational equations, which allows one to circumvent these restrictions and employ equal-order mixed elements. Several numerical examples are used to demonstrate the computationally appealing features of this alternative formulation. The resulting implicit, nonlinear algebraic systems are then solved using an inexact Newton algorithm. We discuss methods for solving the linearized systems using memory-efficient iterative solvers, both on serial and parallel computing platforms. In order to deal with inherent ill-conditioning, we propose a block-structured, multilevel preconditioner that both accelerates the convergence of the Krylov solver and exhibits excellent scaling properties as the number of unknowns and number of processors increase. To demonstrate the effectiveness of these approaches, the analysis framework is applied to modeling hydrologically-driven slope failure. This analysis is motivated by a recent landslide that occurred at a steep experimental catchment (CB1) near Coos Bay, Oregon. Simulations are used to quantify the rainfall-induced slope deformation and assess the failure potential. Results of parametric studies suggest that for a steep hillside slope underlain by shallow bedrock similar to the CB1 site, failure would occur by a multiple slide block mechanism, with progressive failure surfaces forming at the bedrock interface and then propagating to the slope surface. A key observation is that significant computational resources are required to capture these complex solid/fluid interaction mechanisms at sufficient resolution, further justifying the use of the proposed approaches over conventional methods.