## 3 Nonlinear Soil Behavior  
It is well established in geotechnical engineering that soil response is nonlinear beyond a certain level of deformation. Stress-strain relationships for the levels of shear deformation produced by large earthquakes are nonlinear and hysteretic, as has been confirmed by numerous results of vibratory and cyclic loading tests on soil samples.  
The performance of cyclic nonlinear models can be illustrated by a very simple example in which the shape of the backbone curve is described by \(\tau=f(\gamma)\) and depends on the parameters describing the nonlinear behavior of soils [18; 19].  
Among the hyperbolic type models, the most famous and commonly used model is the one initially proposed by [PERSON] (KZ) [20], lately redefined by [PERSON] and [PERSON] [21],  
\[\tau=f(\gamma)=\frac{G\gamma}{1+\frac{G}{\tau_{max}}\mid\gamma\mid}, \tag{21}\]  
where \(\tau\) is the shear stress, \(\gamma\) is the shear strain, \(G\) is the undisturbed shear modulus (taken at the origin), and \(\tau_{max}\) is the shear strength (the maximum stress that material can support in the initial state, \(\tau_{max}=G\gamma_{r}\)). Then, (21) results,  
\[\tau=\frac{G\gamma}{1+\frac{\left\lvert\gamma\right\rvert}{\gamma_{r}}}. \tag{22}\]  
In the numerical tests here addressed, the above equation is rewritten as  
\[\tau=G\gamma-c\frac{\tau\mid\gamma\mid}{\gamma_{r}}, \tag{23}\]  
where \(\gamma_{r}\) is the reference deformation, and \(c\) can be considered as a curve fitting parameter (let us note that, when \(c=0\), the nonlinear case reduces to the linear one).  
The other parameter related with this nonlinearity is the damping characteristics of the soil (represented by the damping ratio \(\zeta\)), which is calculated by the ratio of the energy dissipated \(W_{D}\) and the one associated to the elastic deformation \(W_{E}\), from the hysteretic loop, i.e.,  
\[\zeta=\frac{W_{D}}{4\pi W_{E}}. \tag{24}\]  
[PERSON] proposed, in Reference [22], the following expression of the damping ratio:  
\[\zeta=\frac{2}{\pi}\left(\frac{2\int_{0}^{\gamma}f(\gamma)d\gamma}{f(\gamma) \gamma}-1\right), \tag{25}\]  
that using (22) leads to  
\[\zeta=\frac{4}{\pi}\left(1+\frac{1}{\frac{\gamma_{E}}{\gamma_{r}}}\right) \left[1-\frac{ln\left(1+\frac{\gamma_{A}}{\gamma_{r}}\right)}{\frac{\gamma_{ E}}{\gamma_{r}}}\right]-\frac{2}{\pi}, \tag{26}\]  
with \(\gamma_{A}=max(\gamma)\). Now, considering the linearization introduced in the previous section, it results  
\[\mathbf{M}\frac{d^{2}\mathbf{U}(t)}{dt^{2}}+\mathbf{C}(\mathbf{U})\frac{d \mathbf{U}(t)}{dt}+\mathbf{K}\mathbf{U}(t)=c\mathbf{R}(t)+\mathbf{F}(t), \tag{27}\]  
where \(\mathbf{R}(t)\) account the nonlinear contribution of the constitutive Equation (23), i.e.,  
\[\tau^{nl}=-\frac{\tau\mid\gamma\mid}{\gamma_{r}}, \tag{28}\]  
with \(\tau\) and \(\gamma\) evaluated at the previous iteration of the nonlinear solver. Equation (27) can be expressed in the form of Equation (20) to compute efficiently the parametric solution. When considering a proportional damping, the damping ratio (26) will naturally appear explicitly in the left-hand side of the equation, which will represent a sort of parametric transfer function. In the proposed approach, the iterative procedure is initialized from the small-strain shear modulus and damping. Then, damping is updated at each iteration step, in the left-hand side, from Equation (26). This method fits perfectly with any damping model proposed in the literature. As soon as such parametric solutions are available, the non-linear problem can be solved in real time because no new calculation is needed; the non-linear solver only needs to particularize online the parametric solution calculated offline.