## 2 Mathematical Framework  
### Geomechanical Governing Equations and Model Parameterization  
The geomechanical model of a producing reservoir computes the stress and displacement fields generated by spatial and temporal pore pressure changes due to fluid withdrawal. The analysis is carried out by solving numerically the governing partial differential equations of 3-D poro-mechanics [41; 42] with the aid of the FE method. A one-way coupled approach is used, as it is usually accepted in petroleum engineering. This implies that the pore pressure change in space and time is first computed by a dynamic flow model, and then used as an external source of strength in the equilibrium equations. Let \(\Omega\subset\mathbb{R}^{3}\) denote a 3-D domain in the x-y-z reference system with boundary \(\Gamma\). Neumann and Dirichlet boundary conditions are imposed on boundaries \(\Gamma_{N}\) and \(\Gamma_{D}\), respectively, with \(\Gamma=\Gamma_{N}\cup\Gamma_{D}\) and \(\Gamma_{N}\cap\Gamma_{D}=\emptyset\). The equilibrium equations governing the consolidation of a porous medium at some instant \(t\in I=]0,T[\), with \(T\) the simulation period read:  
\[ \begin{cases}\
abla\cdot\underline{\sigma}-\alpha\
abla p=\mathbf{b}&\text{ on }\Omega\\ \mathbf{u}=\hat{\mathbf{u}}&\text{ on }\Gamma_{D}\\ \underline{\sigma}\cdot\mathbf{n}=\mathbf{f}&\text{ on }\Gamma_{N},\end{cases}  \tag{1}\]  
where \(\underline{\sigma}(\mathbf{x},t)\) is the effective stress tensor, \(\alpha(\mathbf{x})\) the Biot coefficient, \(p(\mathbf{x},t)\) the pore pressure change, \(\mathbf{b}(\mathbf{x},t)\) the external body forces, and \(\mathbf{n}(\mathbf{x})\) the outer normal to \(\Gamma_{N}\). The system of Equations (1) is solved in terms of displacements \(\mathbf{u}(\mathbf{x},t)\), subject to the prescribed displacements \(\hat{\mathbf{u}}(\mathbf{x},t)\) along \(\Gamma_{D}\) and the force per unit surface \(\mathbf{f}(\mathbf{x},t)\) along \(\Gamma_{N}\). To close the system, a non-linear constitutive relationship is introduced to link the effective stress to the strain tensor \(\underline{\epsilon}(\mathbf{x},t)\):  
\[\underline{\sigma}=\boldsymbol{D}\left(\underline{\epsilon}\right), \tag{2}\]  
where \(\underline{\epsilon}\) is related to \(\mathbf{u}\) according to the small strain hypothesis:  
\[\underline{\epsilon}=\frac{1}{2}\left(\
abla+\
abla^{T}\right)\mathbf{u}. \tag{3}\]In this study, the VN [43] visco-elasto-plastic constitutive law is used. Equation (2) is written in differential form as:  
\[\underline{\boldsymbol{\sigma}}=\textbf{C}:\left(\underline{\boldsymbol{\epsilon }}-\hat{\gamma}\frac{\partial p_{c}}{\partial\underline{\boldsymbol{\sigma}}} \right), \tag{4}\]  
where **C** is the standard rank-4 elasticity tensor, \(\hat{\gamma}\) is the plastic multiplier rate and \(p_{c}\) is the plastic potential. The elasticity tensor **C** depends on [PERSON]'s ratio \(\
u\) and a stress-dependent Young's modulus \(E\):  
\[E=\left(1-2\
u\right)\frac{3\rho_{\kappa}p}{\lambda^{*}}, \tag{5}\]  
with \(p\) the volumetric stress (first stress invariant), \(\lambda^{*}\) the modified compression index, and \(\rho_{\kappa}\) the ratio between the modified compression and swelling index.