AbstractIn the present work, a mixed least‐squares finite element method (LSFEM) is used, which is applied on the partial differential equations arising from the Theory of Porous Media (TPM), see for example [1, 2, 3]. Since the LSFEM is not limited to the LBB condition, the method has some theoretical advantages over the well‐known (mixed) Galerkin method, see [4]. In particular, the LSFEM leads to positive definite and symmetric system matrices, even for differential equations with non‐self‐adjoint operators. In detail, we study an incompressible binary model with solid and liquid phases. The modeling of saturated porous structures serves as the basis for the presented approach. The resulting finite element is a four‐field formulation in terms of solid displacements, fluid pressure, mixture stresses, and a new variable associated with the pressure gradient. Vector‐valued Raviart‐Thomas and conventional Lagrangian functions are used to achieve a conforming discretization of the unknowns in the spaces H(div) and H1. The applicability of the LSFEM approach is demonstrated with numerical examples for fluid‐saturated porous structures, which consider incompressible, geometrically nonlinear elastic material behavior.