## 4 Consistent Tangent Stiffness
Furthermore, the overall better performance of a fully consistent Newton-Raphson scheme with respect to its quasi-Newton counterpart was also demonstrated. Methodology and formal analysis: A.G., M.C., and U.P. ; A.C., F.G., S.M. provided the industrial context, the examples in Section 5.3 and carried out the most onerous calculations. All the authors contributed to the discussion, the writing, and the revision of the paper. All authors read and agreed to the published version of the manuscript. This research was funded by ENI S.p.A. The authors declare no conflict of interest. ## Appendix A Deviatoric Projection Operator \(\mathbf{n}\)

The formulation of the implicit time integration scheme in Section 3 makes use of the following key property of the deviatoric projection operator \(\mathbf{n}\):

\[\mathbf{n}^{\text{trial}}=\mathbf{n}^{n+1}=\mathbf{n} \tag{10}\]

This can be verified by considering the definition of the deviatoric stress tensor:

\[\mathbf{s}^{n+1}=\mathbf{Dev}:\sigma^{n+1}=\mathbf{Dev}:\sigma^{\text{trial}} -\mathbf{Dev}:\mathbf{D}:\left(\frac{\partial p^{\text{eq}}(p,q)}{\partial \sigma}\right)^{n+1}\Delta\gamma \tag{11}\]

where the deviatoric operator \(\mathbf{Dev}\) was defined in (42). Making use of the volumetric-deviatoric decomposition (10) of the plastic strain increment, one has:

\[\mathbf{s}^{n+1}=\mathbf{s}^{\text{trial}}-2G\;\mathbf{Dev}:\left(p_{,p}^{ \text{eq}}\mathbf{m}+p_{A}^{\text{eq}}\mathbf{n}\right)^{n+1}\Delta\gamma= \mathbf{s}^{\text{trial}}-2G\left(p_{,q}^{\text{eq}}\mathbf{n}\right)^{n+1}\Delta\gamma \tag{12}\]

Solving with respect to \(\mathbf{s}^{\text{trial}}\) and using the definition (8) of \(\mathbf{n}\), one obtains:

\[\mathbf{s}^{\text{trial}}=\left[\frac{2}{3}q^{n+1}+2G\left(p_{,q}^{\text{eq} }\right)^{n+1}\Delta\gamma\right]\mathbf{n}^{n+1} \tag{13}\]

Squaring and multiplying both sides of the equality by \(3/2\) and taking into account that \(\mathbf{n}:\mathbf{n}=3/2\) give:

\[\frac{3}{2}\mathbf{s}^{\text{trial}}:\mathbf{s}^{\text{trial}}=\frac{3}{2} \cdot\frac{3}{2}\left[\frac{2}{3}q^{n+1}+2G\left(p_{,q}^{\text{eq}}\right)^{n +1}\Delta\gamma\right]^{2} \tag{14}\]

leading to:

\[q^{\text{trial}}=q^{n+1}+3G\left(p_{,q}^{\text{eq}}\right)^{n+1}\Delta\gamma \tag{15}\]

Using the definitions of \(q^{\text{trial}}\) above and of Equations (9), (14), and (15), one can also write:

\[q^{\text{trial}} =\mathbf{n}:\sigma^{n+1}+3G\left(p_{,q}^{\text{eq}}\right)^{n+1} \Delta\gamma= \tag{16}\] Figure 11:0

Having in mind that one can also write \(q^{\text{trial}}=\mathbf{n}^{\text{trial}}:\sigma^{\text{trial}}\), from (16), one obtains:

Figure 11:1



Figure 11:2 Figure 11:3 Figure 11:4 Figure 11:5 Figure 11:6

## References

* [1] [PERSON]; [PERSON] Finite element modelling for simulating the surface subsidence above a compacting hydrocarbon reservoir. _Int. J. Numer. Anal. Methods Geomech._**1994**, _18_, 619-639. [CrossRef]
* [2] [PERSON]; [PERSON]; [PERSON] A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model. _Eur. J. Mech. A/Solids_**2002**, _21_, 1-12. [CrossRef]
* [3] [PERSON]; [PERSON] A fully coupled consolidation model of the subsidence of Venice. _Water Resour. Res._**1978**, _14_, 223-230. [CrossRef]
* [4] [PERSON]; [PERSON] Mathematical simulation of the subsidence of Venice. 1. Theory. _Water Resour. Res._**1973**, \(9\), 721-733. [CrossRef]
* [5] [PERSON]; [PERSON]; [PERSON]; [PERSON] Numerical studies on two-way coupled fluid flow and geomechanics in hydrate deposits. _SPE J._**2012**, \(2\), 485-501. [CrossRef]
* [6] [PERSON]; [PERSON] [PERSON]; [PERSON] Application of a Modified Strain Rate Dependent Constitutive Model for Well-Bore Stability Analyses in Chalk Reservoirs. In Proceedings of the U.S. Rock Mechanics/Geomechanics Symposium, All Days, New York, NY, USA, 28 June 2020; ARMA-2020-1665. * [7] [PERSON]; [PERSON] A modified strain rate dependent constitutive model for chalk and porous rock. _Int. J. Rock Mech. Min. Sci._**2020**, _134_, 104406.