## 4 Consistent Tangent Stiffness
It is also worthwhile to emphasize that there are very small values for \(\delta(\sigma_{\text{approx}})_{11}\) appearing in Figure 8: they are actually on the order of magnitude of the corresponding stress component shown in Table 2, where the increments for the complete stress vector corresponding to Case 4 (with \(\mathbf{n}_{2}=[0\,0\,01\,0]^{T}\)) are also shown. Figure 7: Test 3. Stress point evolution in the \(p-q\) plane. A better correspondence between the two values, to within the round-off machine error, of the local update \(\delta(\sigma_{\text{loc}})_{i}\) and of its estimate \(\delta(\sigma_{\text{approx}})_{i}\) can be obtained only by accounting for the contribution \(\frac{\partial\sigma}{\partial t}\) evidenced in the second term of (60), namely by defining a corrected term \(\delta(\sigma_{\text{corr}})_{i}=\frac{\partial\sigma}{\partial\epsilon}\, \delta\epsilon_{i}+\frac{\partial\sigma}{\partial t}\,\delta t\). In the latter expression, the last contribution accounts for the time dependence of the stiffness matrix in the time step and cannot be inserted in standard FE codes, where typically, user customization of the constitutive law allows only providing in a user subroutine the term \(\frac{\partial\sigma}{\partial\epsilon}\), next passed to the global assembly. This very small discrepancy, of the order of 0.0005% for the considered cases, however, turned out to be appreciable only in the mentioned verification cases with tiny strain increments (order of \(10^{-8}\) in this case) that were not of interest in real applications, and this is irrelevant for practical applications. Nevertheless, it needed to be carefully considered during the validation phase. ### Subsidence Evaluation for a Real Case Study

A real gas field, located in the Adriatic Basin, was selected for the last test in order to compare not only the results obtained with the different formulations, but also the computational times on a larger scale model, representative of the size of the studies typically performed in the industry for the prediction of subsidence due to hydrocarbon exploitation. The field is located at about 60 km from the Italian coastline, where the average water depth is 60 m. The depth of the reservoir layers ranges from 900 to 1500 m asl. The gas volume originally in place is approximately 30 GSm\({}^{3}\), with an expected recovery factor of 50%, and gas is produced from 28 wells, connected to two platforms. The geomechanical model, whose details can be found in [31], was built from the corresponding reservoir model in the same way as in the synthetic case described in [26]. The geometry of the fluid-dynamic model was extended both laterally and vertically, so that the geomechanical model reached an areal extent of about \(88\times 73\) km, for a depth of 5 km. It consisted of about 5.5 \(\times\)\(10^{5}\) finite elements, for a total of around \(2\times 10^{6}\) degrees of freedom. Parameters for the

\begin{table}
\begin{tabular}{c c c c} \hline \hline  & \(\delta(\sigma_{\text{loc}})\) & \(\delta(\sigma_{\text{approx}})\) & \(\delta(\sigma_{\



Vermeer-Neher constitutive law used in the simulation are reported in Figure 3.0 Preliminary values for the parameters were obtained from the interpretation of specifically tailored laboratory tests on samples from one of the wells of the field, and subsequently calibrated to reproduce the subsidence measurements provided by a GPS station installed on a platform of the field. The simulation was protracted for 30 years after the end of production in order to evaluate the effect of pressure evolution in the mineralized region and in the connected aquifers after the end of production. Figure 9 shows the time evolution of vertical displacement at the point of maximum subsidence. Due to confidentiality reasons, values were normalized with respect to the maximum value of subsidence reached during the simulation.