## 4 Consistent Tangent Stiffness
The implemented procedure is summarized in Algorithms 1-4. ```
Initialization  Solve Newton-Raphson linear System (35) to obtain \(\sigma^{n+1}\), Equation (21)  Compute tangent matrix \(\left(\frac{\partial\sigma}{\partial\epsilon}\right)^{n+1}\), Equation (60)
```  
**Algorithm 1**Implicit backward difference time integration. ```
At time step \(n\), from \(t_{n}\) to \(t_{n+1}\), all quantities at time \(t_{n}\) are known. Given \(\Delta\epsilon^{n}\), compute \(\epsilon^{n+1}=\epsilon^{n}+\Delta\epsilon^{n}\). \(\mathbf{n}^{\text{trial}}\), \(\sigma^{\text{trial}}\), \(\mathbf{s}^{\text{trial}}\), \(p^{\text{trial}}\) and \(q^{\text{trial}}\) from Equations (8)\({}_{2}\) and (29). Initial (trial) values for derivatives of \(p^{\text{eq}}\) are then known from Equation (A10). ```  
**Algorithm 2**Initialization. ```
At iteration \(i\), if \(i=0\), use trial values for \(n+1\) quantities. Compute residuals \(\boldsymbol{\
abla}_{i}\) for the linearized system (35) with Equations (23)-(27). Compute the coefficient matrix \(\left.\frac{d\boldsymbol{\
abla}}{dX}\right|_{i}\) for the linearized system (35) with Equations (A11)-(A13). Solve linearized System (35) for the unknowns \(\delta\boldsymbol{X}\) in Equation (36). Update \(\boldsymbol{X}^{i+1}=\boldsymbol{X}^{i}+\delta\boldsymbol{X}\). Convergence test: if passed, \(\boldsymbol{X}^{n+1}=\boldsymbol{X}^{i+1}\), and exit loop. Else, repeat from 1 with \(\boldsymbol{X}^{i+1}=\boldsymbol{X}^{i}\) and \(i+1\longrightarrow i\). ```  
**Algorithm 3**Solve Newton-Raphson linear System (35) to obtain \(\sigma^{n+1}\). ```  At iteration \(i\), if \(i=0\), use trial values for \(n+1\) quantities. Compute residuals \(\boldsymbol{\
abla}_{i}\) for the linearized system (35) with Equations (23)-(27). Compute the coefficient matrix \(\left.\frac{d\boldsymbol{\
abla}}{dX}\right|_{i}\) for the linearized system (35) with Equations (A11)-(A13). Solve linearized System (35) for the unknowns \(\delta\boldsymbol{X}\) in Equation (36). Update \(\boldsymbol{X}^{i+1}=\boldsymbol{X}^{i}+\delta\boldsymbol{X}\). Convergence test:

## 5 Numerical Tests

The robustness, accuracy, and efficiency of the proposed model were investigated with numerical tests. The implicit version of the model described in the previous sections was first implemented in a MATLAB(tm) code and validated at the material point level by imposing a strain history and calculating the corresponding stress. Then, the material model was implemented in a user material subroutine of a commercial finite element code [30] and was applied to the simulation of subsidence at a reservoir scale, as was described in [26]. Finally, the model was tested on a large-scale geomechanical problem. In the tests, the proposed implicit approach was compared both with the explicit integration procedure proposed in [13] and with the implicit one proposed in [15]. When showing the results, the following convention for the stress components was assumed: \(\mathbf{\sigma}=[\sigma_{11}\ \sigma_{22}\ \sigma_{33}\ \sigma_{23}\ \sigma_{13}\ \sigma_{12}]\). ### Tests at the Material Point Level

The tests at the material point level were performed by imposing a strain history and evaluating the corresponding stresses. The results of the proposed implicit integration scheme were compared with the outcomes of an explicit integration. The material properties used in these tests are shown in Table 1, where OCR denotes the overconsolidation ratio, used in the definition of \(p_{0}^{\text{eq}}=\text{OCR}\cdot p_{0}\), \(p_{0}\) being the initial hydrostatic stress.