Abstract. Rocks break if shear stresses exceed their strength. It is therefore important for typical geoscientific applications to take shear failure mechanism and the subsequent development of mode-II shear bands or faults into account. Many existing codes incorporate non-associated Drucker-Prager or Mohr-Coulomb plasticity models to simulate this behavior. Yet, when effective mean stress becomes extensional, for example when fluid pressure becomes large, the dominant failure mode changes to a mode-I (opening) mode, which initiates plastic volumetric deformation. It is rather difficult to represent both failure modes in numerical models in a self-consistent manner, while also accounting for the nonlinear visco-elastic host rock rheology, which varies from being nearly incompressible in the mantle to being compressible in surface-near regions. Here, we present a simple plasticity model that is designed to overcome these difficulties. We employ a combination of a linearized Drucker-Prager shear failure envelope with a circular tensile cap function in way that ensures continuity and smoothness of both yield surface and flow potential in the entire stress space. A Perzyna-type viscoplastic regularization ensures that the resulting localization zones are mesh-insensitive. To deal with the near incompressibility condition, a mixed two-field finite element formulation is employed. The local nonlinear iterations at the integration-point level are used to determine the stress increments. The global Newton-Raphson iterations are applied to solve the discretized momentum and continuity residual equations. The presented plasticity model is implemented in an open-source 2D unstructured finite element code GeoTech2D. The results of several typical test cases that range from crustal scale deformation to the propagation of fluid-induced tensile failure zones demonstrate rapid convergence. The robustness of the solution scheme is enhanced by the adaptive time stepping algorithm.