In this work, prestressed elasticity problem as a model of the so-called glacial isostatic adjustment (GIA) process is studied. The model problem is described by a set of partial differential equations (PDE) and discretized with a mixed finite element (FE) formulation. In the presence of prestress the so-constructed system of equations is non-symmetric and indefinite. Moreover, the resulting system of equations is of the saddle point form. We focus on a robust and efficient block lower-triangular preconditioning method, where the lower diagonal block is and approximation of the so-called Schur complement. The Schur complement is approximated by the so-called element-wise Schur complement. The element-wise Schur complement is constructed by assembling exact local Schur complements on the cell elements and distributing the resulting local matrices to the global preconditioner matrix. We analyse the properties of the element-wise Schur complement for the symmetric indefinite system matrix and provide proof of its quality. We show that the spectral radius of the element-wise Schur complement is bounded by the exact Schur complement and that the quality of the approximation is not affected by the domain shape. The diagonal blocks of the lower-triangular preconditioner are combined with inner iterative schemes accelerated by (numerically) optimal and robust algebraic multigrid (AMG) preconditioner. We observe that on distributed memory systems, the top pivot block of the preconditioner is not scaling satisfactorily. The implementation of the methods is further studied using a general profiling tool, designed for clusters. For nonsymmetric matrices we use the theory of Generalized Locally Toeplitz (GLT) matrices and show the spectral behavior of the element-wise Schur complement, compared to the exact Schur complement. Moreover, we use the properties of the GLT matrices to construct a more efficient AMG preconditioner. Numerical experiments show that the so-constructed methods are robust and optimal.