Most of the land subsidence numerical models introduced in the literature are based on limiting assumptions, such as one-dimensional land movement, uncoupled formulations, and linear elasticity; and some even does not satisfy mass conservation requirements. As a result, many of the models may fail to realistically and efficiently represent the behavior of the porous media undergoing hydraulic loading. In this study, a locally conservative numerical technique is proposed for the analysis of land subsidence in saturated porous media in both the plane strain and axisymmetric conditions. To this end, the flow equation, discretized using a dual discrete finite volume method, is coupled iteratively with the elastoplastic deformation model, solved using the finite element method. The developed model is then verified, and its benefits compared to the standard finite element method, are discussed. The application of the proposed method is further demonstrated through a quasi-real conceptualization of a multi-aquifer system, considering the effect of different seasonal pumping schedules. The results show that soil deformations in the horizontal direction may be as large as those in the vertical direction and cannot be neglected. It is also shown that the pumping schedule can have considerable effects on the land subsidence regime.