A mixed least-squares finite element method (LSFEM) is proposed in this work for the geometrically nonlinear theory of porous media (TPM) with incompressible solid skeleton and viscous fluid [1]. In particular, a fully saturated solid-fluid binary system is studied. The deformation of solid skeleton is governed by the usual Cauchy’s linear momentum equation and Darcy’s law describes the fluid flow in pores. This forms the well known first order Div- Grad system of governing equations, suitable for least-sqaures functional minimization process. The stress and deformation of the solid phase, along with the pore pressure and velocity of the liquid phase, constitute four independent fields in this formulation. The normal-flux preserving, vector valued Raviart–Thomas (RT) functions are employed to discretize the stress and velocity fields, while the displacement and pressure fields are approximated using conventional Lagrange polynomials (P). The LSFEM avoids the intricacies associated with satisfying the inf–sup (or LBB) condition required in mixed Galerkin finite element methods [2] to ensure numerical stability. Another key advantage is that the LSFEM results in symmetric positive-definite system of matrices - even for non-self-adjoint differential operators. The error minimization approach adapted in LSFEM makes the role of residual weights crucial in deciding the convergence characteristics. Therefore, in this study, the performance of higher order RT-P elements is evaluated against the canonical consolidation problem. A parametric study is carried out to describe the effect of respective weights on the convergence of analysis. The behavior of each residual with the usage of different combination of RT and P functions is also brought into context. The suitability of the developed LSFEM for the simulation of fluid saturated porous media with incompressible constituents is thus demonstrated.