In finite element software one has to solve a system of non-linear equations, which is commonly simplified to a sequence of linear system. We research the possibility to solve these systems on a GPU to improve the solve time. We are particularly interested in systems arising from geotechnical models. We compare several combinations of Krylov methods, parallel preconditioners and deflation methods and present a suitable combination. This solver is then compared with existing CPU based solvers in PLAXIS 3D. We show that compared to the current iterative solver, the iteration time can be reduced by 50% up to 85% depending on the problem. While compared to the current direct solver, the memory consumption and initialization time can be reduced significantly.