# Numerical Simulation of the Influence of the Baihetan Reservoir Impoundment on Regional Seismicity
## 2 Pore-Viscoelastic Method
Next, we transform Equation (5) to simplify it:  
\[\sigma_{ij}=2\big{(}g(t-\xi)\epsilon_{ij}(\xi)\big{)}\big{|}_{- \infty}^{t}+2\int_{-\infty}^{t}\frac{\partial g(t-\xi)}{\partial(t-\xi)} \epsilon_{ij}(\xi)d\xi+\big{(}K-\frac{2}{3}g(t-\xi)\epsilon_{kk}(\xi)\big{)} \big{|}_{-\infty}^{t}\delta_{ij} \tag{6}\] \[-\frac{2}{3}\int_{-\infty}^{t}\frac{\partial g(t-\xi)}{\partial( t-\xi)}\epsilon_{kk}(\xi)d\xi\delta_{ij}+\Big{[}-\frac{3(v_{u}-v)}{B(1+v)(1+v_{u})}p \Big{]}\delta_{ij}\] \[=2g(0)\epsilon_{ij}(t)+\big{[}K-\frac{2}{3}g(0)\big{]}\epsilon_{ kk}(t)\delta_{ij}+\Big{[}-\frac{3(v_{u}-v)}{B(1+v)(1+v_{u})}p\Big{]}\delta_{ij}+r_{ij}\]  
where  
\[r_{ij}=2\int_{-\infty}^{t}\frac{\partial g(t-\xi)}{\partial(t-\xi)}\epsilon_{ ij}(\xi)d\xi-\frac{2}{3}\int_{-\infty}^{t}\frac{\partial g(t-\xi)}{\partial(t-\xi)} \epsilon_{kk}(\xi)d\xi\delta_{ij}=2\int_{-\infty}^{t}\Phi(t)\epsilon_{ij}(\xi )d\xi-\frac{2}{3}\int_{-\infty}^{t}\Phi(t)\epsilon_{kk}(\xi)d\xi\delta_{ij} \tag{7}\]  
where  
\[\sigma_{ij}=2G\epsilon_{ij}+\frac{v}{1+v}(\sigma_{kk}+3p)\delta_{ij}-\frac{3(v_{ u}-v)}{B(1+v)(1+v_{u})}p\delta_{ij} \tag{1}\]0  
Then we make some changes to \(r_{ij}\), as \(g(t)=Ge^{-\frac{G}{\eta}t}\), and use the theorem  
\[ \begin{array}{c}\frac{d}{dt}\int_{a(t)}^{b(t)}f(x,t)dx=\int_{a(t)}^{b(t)} \frac{\partial f(x,t)}{\partial t}dx+f[b(t),t]\frac{db(t)}{dt}-f[a(t),t]\frac{ da(t)}{dt}\\ \frac{\partial}{\partial t}\int_{-\infty}^{t}\Phi(t)\epsilon_{ij}(\xi)d\xi= \int_{-\infty}^{t}\frac{\partial\Phi(t)}{\partial t}\epsilon_{ij}(\xi)d\xi+ \Phi(t-t)\epsilon_{ij}(t)\frac{\partial t}{\partial t}-\Phi(t+\infty)\epsilon _{ij}(-\infty)\frac{\partial(-\infty)}{\partial t}\\ =\int_{-\infty}^{t}\frac{\partial\Phi(t)}{\partial t}\epsilon_{ij}(\xi)d\xi+ \Phi(0)\epsilon_{ij}(t)\\ =-\frac{G}{\eta}\int_{-\infty}^{t}(\Phi)\epsilon_{ij}(\xi)d\xi+\Phi(0) \epsilon_{ij}(t)\end{array}  \tag{10}\]Thus, we can transform Equation (7) to  
\[ \begin{split}\frac{\partial r_{ij}}{\partial t}&=2\Phi(0) \epsilon_{ij}(t)-\frac{2}{3}\Phi(0)\epsilon_{kk}(t)\delta_{ij}-\frac{G}{\eta}r _{ij}(t)\\ &=-2\frac{G^{2}}{\eta}\epsilon_{ij}(t)-3\frac{G^{2}}{\eta}\epsilon _{kk}(t)\delta_{ij}-\frac{G}{\eta}r_{ij}(t)\end{split}  \tag{11}\]  
For quasi-static phenomena, we use \(\sigma_{ij}\) and p to express the governing equations, where the stress \(\sigma_{ij}\) must satisfy the equilibrium equation. \[\frac{\partial\sigma_{ij}}{\partial x_{j}}-f_{i}=0,\sigma_{ij}=\sigma_{ji} \tag{12}\]  
in weak form,  
\[ \begin{split}\int_{V}\big{(}\sigma_{xx}\delta\epsilon_{xx}+ \sigma_{yy}\delta\epsilon_{yy}+\sigma_{zz}\delta\epsilon_{zz}+2\sigma_{xy} \delta\epsilon_{xy}+2\sigma_{xz}\delta\epsilon_{xz}+\sigma_{yz}\delta \epsilon_{yz}\big{)}dV\\ &=\int_{V}\big{(}f_{x}\delta u+f_{y}\delta v+f_{z}\delta w\big{)} dV+\int_{\Gamma}\big{(}T_{x}\delta u+T_{y}\delta v+T_{z}\delta w\big{)}d \Gamma\end{split}  \tag{13}\]  
where \(f_{i}\) denotes the internal force such as gravity, and the \(T_{i}\) is the face force on the surface. Using a forward difference scheme in the time domain \(\frac{\partial F}{\partial t}=\frac{F-F^{x}}{dt}\), and according to Equations (2) and (11)-(13), we can get the finite element weak solution form of the equations:  
\[ \begin{split}\kappa\big{(}\frac{\partial p}{\partial x_{i}}, \frac{\partial\wp}{\partial y_{i}}\big{)}+\frac{9(v_{ix}-v)}{2C\Theta^{2}(1+v )(1+v_{ix})dt}(p,\wp)\\ =\frac{9(v_{ix}-v)}{2C\Theta^{2}(1+v)(1+v_{ix})dt}(p^{n},\wp) \\ -\frac{3(v_{ix}-v)}{2C\Theta(1+v)(1+v_{ix})dt}(\sigma_{kk}- \sigma_{ik}^{n},\wp)+\kappa\Big{(}\frac{\partial p}{\partial n},\wp\Big{)} \end{split}  \tag{14}\]  
\[