The coupled hydro-mechanical state in geomaterials undergoing plasticity phenomena is here evaluated by means of the subloading surface model . The most important feature of this theory is the abolition of the distinction between the elastic and plastic domain, as it happens in conventional elastoplastic models. This means that plastic deformations are generated whenever there is a change in stress and a smoother elasto-plastic transition is produced. The subloading surface takes the role of a loading surface which always passes through the current stress point σ and keeps a shape similar to that of the normal yield surface and a similar orientation with respect to the origin of stress space. Additionally, the model allows for giving a smooth response in a smooth monotonic loading process and the stress is automatically drawn back to the normal-yield surface even if it goes out from that surface, leading to a more stable and robust calculation even for large loading steps. The plasticity algorithm has been implemented within the FE PLASCON3D research code, coupling hydro-(thermo)-mechanical fields within a saturated porous medium (locally partially saturated) subjected to external loads. Applications to soils allow e.g. for assessing subsidence evolution at regional scale.