Strain localization is ubiquitous and controls the development of shear zones and the establishment of plate boundaries. The outermost portion of the lithosphere, the crust, is relatively cold. In these regions, deformation is governed by elasto-plasticity, which allows tectonic plates to behave both as rigid solid and, locally, as a deformable solid. The rocks that make up the crust exhibit a frictional behavior, which causes their strength to depend on pressure (Byerlee, 1978). Moreover, the dilatancy angles of rocks are generally smaller than friction angles (Zhao & Cai, 2010), therefore frictional plasticity models should be non-associated. Geodynamic models are typically employed to simulate the deformations of the Earth's upper shell, such as the opening of rifts zone (e.g., Naliboff et al., 2017), the growth of orogenic wedges (e.g., Buiter et al., 2016), or the formation of transform boundaries (e.g., Gerya, 2013) heavily rely on the simulation of frictional deformation processes. Nonetheless, the inclusion of frictional plasticity in numerical models poses problems as standard model implementations may not achieve force equilibrium and may exhibit mesh dependence. The absence of characteristic spatial or temporal scales in standard frictional constitutive laws is the primary cause for these problems. To overcome these issues, models should be augmented by including additional physics (e.g., fluid pressure diffusion, viscoplastic dissipation) (e.g., Brantut et al., 2017), which from a mathematical point of view can be considered as regularization techniques (e.g., de Borst et al., 1993).Numerous regularization strategies have been proposed, which have the effect of introducing a length or a time scale in the model. Viscoplastic regularization relies on the inclusion of a time scale W. M. Wang et al., 1997) which implicitly rather than directly introduces a length scale like in (e.g., H. Wang, 2019). The yield function does not depend on the gradients of plastic strain: it is defined locally. As a consequence, return mapping can be achieved without solving an additional partial differential equation. Thus, the total number of degrees of freedom is not increased. The implementation of viscoplastic regularization in existing codes is