The consolidation theory was developed in a three-dimensional (3D) setting by Biot, giving rise to a system of PDEs that combines the elastic equilibrium of the porous body with the continuity of ground water flow. Because of the mathematical complexity of Biot’s consolidation theory, analytical solutions are possible only for simple problems with a symmetric configuration, while for more complex realistic situations it is necessary to use numerical techniques. At present the implementation of a finite element (FE) consolidation model can be considered a quite common practice. Nevertheless, the solution to realistic problems may still be very difficult, particularly because the resulting linear system can be ill-conditioned, especially for small time steps, and typically has a large size, thus requiring significant computational cost and effort. The present paper addresses some difficulties which arise from the FE so-lution to the coupled consolidation equations and suggests efficient numerical techniques for solving large size realistic models. A criterion for the critical time step leading to ill-conditioning is developed and compared with the accuracy condition by Vermeer and Verruijt. Finally, the use of iterative methods is explored, with a special attention to reordering, preconditioning, and scaling strategies which can accelerate the solver convergence.