## 1 Introduction  
The computation and prediction of subsidence, due to the exploitation of an underground deposit of hydrocarbons, is generally a very complex task and involves advanced numerical simulations [1, 2]. The problem physics includes the evolution of the fluid pressure in the deposit and the correct evaluation of the time-varying stress field in the compacting soil/rock. The corresponding fluid and solid problems can be addressed in a coupled [3] or uncoupled approach [4] and other variants, including for example sequential coupling [5]. In the latter case, for a domain treated as a solid continuum, nonlinear constitutive models accounting for creep become necessary. The several stress-strain time-dependent models for soils that have been presented in the literature mainly fall within the following two main categories:(i) models based on the concept of overstress and (ii) models based on nonstationary yield theory, where the classical plasticity yield limit is generalized to a viscoplastic yield locus that depends on a time-dependent function, see e.g., [6, 7] and the reviews in [8, 9]. By following the overstress approach, the Vermeer-Neher (V-N) model [10, 11], which addresses materials with a high degree of compressibility, such as soft soils, and generalizes odometer test results to fully three-dimensional conditions accounting also for (secondary) creep through an elastic-viscoplastic model, has encountered a significant popularity and is still actively used in the oil and gas industry for subsidence modeling, to predict the deformation of the ground surface induced by hydrocarbon withdrawal from underground reservoirs. The reason for its success in engineering practice [12, 13] rests on its relative simplicity and the limited number of required parameters that can be obtained by a few standard experimental tests. As any other elastic-viscoplastic model, the V-N model requires a time integration scheme to be implemented in a finite element code. Both explicit [12; 13] and (semi-)implicit schemes [10; 14] were used in the past for the V-N model time integration at the local level. For the semi-implicit time integration approaches [10; 14], there are unsolved issues: namely, the deviatoric strain increment is not accounted for [14], so that the integration is not actually fully implicit, and in general, the convergence is at best sublinear for the time integration at the global level. Very recently [15], a fully implicit backward Euler formulation for the local time integration of the V-N model was proposed, in combination with a preconditioned conjugate gradient (PCG) technique for the global solver. In this formulation, the tangent stiffness matrix, however, is artificially symmetrized (this is a requirement for the PCG method), and its entries are computed numerically, without providing their analytical expression. Despite the popularity of the V-N model, its intense use in practical engineering applications, and the fact that its implementation is still the object of new contributions in the scientific literature, a robust and completely consistent implicit backward Euler integration of the model, together with an analytical expression of its consistent tangent matrix, is still missing. The reason for this deficiency is probably due to the technical difficulties involved in the formulation of the Newton-Raphson local iterative scheme, necessary to obtain the updated stress and viscoplastic strain increments, and most of all, in the laborious derivation of the analytical expression of the consistent tangent matrix for the global Newton-Raphson scheme. A consistent tangent matrix formulation and Newton-Raphson scheme for the implicit integration of an elastoplastic finite element model in large strains was proposed for the first time by [PERSON] [16]. [PERSON] and [PERSON] [17] clarified the difference between the tangent matrix for the rate problem and the consistent tangent matrix for the backward Euler iterative scheme. [PERSON] et al.