## 2 Finite Element Modeling
### Modeling Methodology  
A nonlinear elastoplastic model was adopted for the backfill materials to predict the responses of surrounding soils more accurately under various loading conditions. For the nonlinear elastic components, the hyperbolic stress-strain relationship proposed by [PERSON] and [PERSON] [15] was adopted. Since material parameters change with material-stress states, an incremental approach was used to approximate nonlinear behaviors as piecewise linear behaviors. During the application of each load increment, the materials were considered to be linearly elastic. However, different tangential elastic moduli were applied according to the changes in stresses at each incremental load. The tangential modulus (\(E_{t}\)) is defined as follows:  
\[E_{t}=E_{i}\bigg{[}1-\frac{R_{f}(1-\sin\phi)(\sigma_{1}-\sigma_{3})}{2c\cos \phi+2\sigma_{3}\sin\phi}\bigg{]}^{2} \tag{1}\]  
where \(R_{f}\) is the failure ratio, which normally ranges from 0.7 to 1.0. \(E_{i}\) is the initial modulus, which is equal to \(KP_{a}=\left(\frac{\sigma_{3}}{P_{a}}\right)^{n}\). \(P_{a}\) is the atmospheric pressure (101.3 kPa). \(K\) and \(n\) are material constants. \(c\) is the cohesion. \(\phi\) is the friction angle. \(\sigma_{1}\) and \(\sigma_{3}\) are the major and minor principal stresses, respectively. Once the stresses are obtained at each load increment, the tangential modulus can be calculated for any stress state. For the plastic components, the Drucker-Prager yield criterion with a non-associated flow rule was adopted. The Drucker-Prager yield criterion can be expressed in terms of deviatoric stress and hydrostatic stress as follows:  
\[f(s,p)=\parallel s\parallel+\sqrt{2}(3\alpha p-k)=0, \tag{2}\]where \(s\) and \(p\) are the deviatoric and hydrostatic stress, respectively, \(\parallel s\parallel\) is the norm of the deviatoric stress, and \(\alpha\) and \(k\) are material constants that can be expressed in terms of the friction angle \(\phi\) and cohesion \(c\) of the soil, where  
\[\alpha=\frac{2\sin\phi}{\sqrt{3}(3-\sin\phi)}\text{ and }k=\frac{6c\cos\phi}{ \sqrt{3}(3-\sin\phi)}. \tag{3}\]  
A non-associated flow rule is then applied by defining the plastic potential as follows:  
\[g(\sigma)=\parallel s\parallel, \tag{4}\]  
which yields the plastic strain rate \(\dot{\epsilon}^{p}\) in the following form:  
\[\left\{ \begin{array}{cc}\dot{\epsilon}^{p}=\dot{\lambda}\frac{3\alpha g}{ \omega}=\dot{\lambda}n&for\ loading\\ \dot{\epsilon}^{p}=0&otherwise\end{array} \right., \tag{5}\]  
where \(n=\frac{s}{\parallel s\parallel}\). The plastic strain rate parameter \(\dot{\lambda}\) is a positive scalar of proportionality during plastic deformation and the overdot denotes an increment. The aforementioned nonlinear elastoplastic model was implemented in ABAQUS using a user-defined material subroutine called UMAT. The buried concrete pipe was characterized as a linearly elastic plastic material. A damage model was incorporated into the plastic component, where reinforcement behaviors were modeled by considering the tension stiffening effect. The damage-related failure mechanisms of the concrete (cracking and crushing) cause degradation of the elastic stiffness, as expressed by the following degradation scalar:  
\[D^{e}=(1-d)D_{0}^{e}, \tag{6}\]  
where \(d\) is a degradation variable ranging from zero (undamaged) to unity (fully damaged) and \(D_{0}^{e}\) is the initial elastic stiffness. The yield function of the plastic damage model was proposed by [PERSON] et al. (2016) and modified by [PERSON] (2017).