Advanced dynamic nonlinear schemes for geotechnical earthquake engineering applications: a review of critical aspects

Nonlinear time domain numerical approaches together with elasto-plastic effective stress soil constitutive models are nowadays available to geotechnical researchers and practitioners interested in geotechnical earthquake engineering. The use of such advanced two- and three-dimensional schemes allows the analysis and design of complex geotechnical structures within a performance-based framework, considering the build-up of excess pore water pressure during the earthquake, dynamic interaction between the soil deposit and the above surface buildings and infrastructures and effects of multi-directional seismic loading. Within this context, the paper focuses on the dynamic finite element (FE) method and presents a review of the key ingredients governing its predictive capabilities. These include i) the description of the fully coupled solid–fluid interaction formulation and time integration, ii) the calibration of Rayleigh damping and the soil constitutive model parameters, iii) prescription of the boundary conditions for the generated mesh and iv) input motion selection/scaling strategies. For each of the above points, the paper summarises the current knowledge and best practice, with the aim of providing protocols for a confident application of nonlinear FE schemes to evaluate the performance of critical geotechnical infrastructures. Useful hints to promote familiarity of advanced nonlinear soil dynamic analysis among geotechnical practitioners and to indicate areas for further improvement are also highlighted.


Introduction
The design of geotechnical structures subjected to earthquake loading can be conducted using a variety of analytical and numerical approaches, which are characterised by different levels of complexity. Simplified dynamic analyses are based on the pseudo-static method, where the seismic action can be obtained from preliminary free-field ground response simulations. To predict site effects, seismic input motions are propagated through the deposit approximating the dynamic characteristics of the soil by means of equivalent linear or nonlinear models Abstract Nonlinear time domain numerical approaches together with elasto-plastic effective stress soil constitutive models are nowadays available to geotechnical researchers and practitioners interested in geotechnical earthquake engineering. The use of such advanced two-and three-dimensional schemes allows the analysis and design of complex geotechnical structures within a performance-based framework, considering the build-up of excess pore water pressure during the earthquake, dynamic interaction between the soil deposit and the above surface buildings and infrastructures and effects of multidirectional seismic loading. Within this context, the paper focuses on the dynamic finite element (FE) method and presents a review of the key ingredients governing its predictive capabilities. These include i) the description of the fully coupled solid-fluid interaction formulation and time integration, ii) the (e.g. Elia 2015). In the equivalent linear visco-elastic method (Schnabel et al. 1972;Idriss et al. 1973), the shear modulus (G) and damping ratio (D) are iteratively updated in order to approximate the soil nonlinear behaviour through a sequence of linear analyses. For each iteration, the exact solution proposed by Roesset (1977) is used within a total stress approach. Nonlinear time-domain schemes are also widely adopted in the literature for one-dimensional ground response analyses. Some of them employ the hyperbolic model, in its original or extended pressure dependent formulation (e.g. Hashash et al. 2011), but still adopting a total stress approach. In others, nonlinearity is simulated by incremental plasticity models, which allow to predict the permanent deformation of a soil deposit, the generation of excess pore pressure, the decay of the soil stiffness and the corresponding hysteretic dissipation induced by the earthquake action (e.g. Elgamal et al. 2006).
In addition to the classical two-dimensional (2D) and three-dimensional Finite Element Method (FEM) (e.g. Zienkiewicz et al. 1977;Potts and Zdravković 1999) and Finite Difference Method (FDM) (e.g. Smith 1985;Lemos 2012), numerical approaches also include the Boundary Element Method (e.g. Banerjee and Butterfield 1981;Venturini 1983;Shinokawa and Mitsui 1993) and the Discrete Element Method (e.g. Cundall and Strack 1979;O'Sullivan 2017). Through a 2D or 3D numerical discretisation of the boundary value problem (BVP), the analysis can account for i) the interaction during earthquakes between buried complex morphologies and surface or subsurface buildings and infrastructures (i.e. dynamic soil-structure interaction) and ii) the effects of the multi-directional seismic loading conditions. This come at the expense of increasing the computational cost, especially when 3D geometries are considered, and a thorough understanding of the mathematical formulation implemented in the software. In addition, the soil mechanical behaviour can be described using advanced elasto-plastic soil constitutive laws, which are able to predict within a coupled effective stress formulation the solid-fluid interaction (Biot 1941;Zienkiewicz et al. 1999). Advanced constitutive models are, indeed, increasingly available in commercial FE and FD codes (e.g. Mazzoni et al. 2006;Itasca Consulting Group Inc 2016;Plaxis 2019). Nevertheless, their calibration is not straightforward and requires special in-situ and laboratory data (e.g. geophysical tests, resonant column/torsional shear, cyclic triaxial, bender element tests), often not included in a standard geotechnical characterisation. Therefore, advanced dynamic analyses are adopted by non-expert users with difficulty, as both the model calibration procedures and the software usage protocols are unclear or poorly documented, leading to uncertain numerical results and, as such, obscuring the possible benefits associated with their use in the geotechnical design process (e.g. Kwok et al. 2007;Amorosi et al. 2010). Nonetheless, the study of wellinvestigated and monitored case histories, where an accurate characterisation of the soil dynamic properties is available, can effectively be used to inform calibration of the adopted constitutive model and validation of numerical results.
This paper focuses on the application of the FE method in the analysis and design of geotechnical BVPs subjected to seismic loading, covering a review of the key ingredients governing its predictive capabilities. In a sequential order, the work describes and discusses i) fully coupled solid-fluid interaction formulation and time integration of the dynamic equations of motion, ii) setting of the Rayleigh viscous damping, iii) calibration of the soil constitutive model parameters and the initialisation of their state variables, iv) definition of the dynamic boundary conditions along the sides of the numerical model and v) input motion selection/scaling strategies, as they all can play a crucial role in the numerical results. Useful hints for each of the above points, based on the most recent research, are provided to promote the adoption of dynamic FE schemes with more confidence.

Fully coupled solid-fluid interaction formulation
Considering the soil as a porous material with a single fluid filling the pores, the mechanical behaviour can be described by the effective stress carried by the solid skeleton, defined as (Biot 1941): where σ is the total stress vector, p the pore pressure and m = [1 1 1 0 0 0]. In Eq. (1), α is a factor which becomes close to unity when the bulk modulus of the grains is much higher than that of the whole material (i.e. α ≈ 1 for soils). Geotechnical calculation models in special cases can be restricted to either fully drained or undrained conditions; however, the real soil behaviour is time dependent, with pore pressures related to the soil permeability, rate of loading and hydraulic boundary conditions. Hence, the numerical analysis of time dependent problems (such as those related to soil dynamics) requires simultaneous calculation of the soil deformation and groundwater flow by coupling the seepage equations with the equilibrium and constitutive equations. The fully coupled dynamic formulation for solid-fluid interactions stems from the general 3D Biot's theory of consolidation [1941], combining force equilibrium and flow continuity equations. Assuming that the relative acceleration of the fluid phase with respect to the solid skeleton and its convective part are negligible, the resulting system of ordinary equations of motion can be written as follows (Chan 1988;Zienkiewicz et al. 1999): where u is the solid phase displacement vector and p is the pore fluid pressure vector, M is the mass matrix, K is the stiffness matrix, C is the viscous damping matrix, Q is the coupling matrix between the motion and flow equations, H is the permeability matrix, S is the fluid compressibility matrix, f p is the force vector for the fluid phase and f s is the force vector for the solid phase. The system (2) of mixed equations of displacements and pore pressures, valid for fully saturated soils, represents the so-called "u-p formulation" for geotechnical dynamic analysis. System (2) reduces to the general consolidation theory if the inertia of the soil is neglected (i.e. ̈ = 0 ), to the seepage equation if also ̇ = = 0 and ̇ = 0 , and to the static conditions when ̈ =̇ = 0 and ̇ = 0. Zienkiewicz et al. [1999] showed that the u-p approximation is suitable for the analysis of dynamic problems when the frequency of the input is within a typical range of frequencies characterizing real earthquakes (i.e. 0 ÷ 20 Hz). If shocks (i.e. explosions) or very high frequencies are involved, the exact Biot's solution accounting for all the acceleration terms should be adopted.
(2) M̈ + Ċ + Ku-Qp = s Q̇ + Ṡ + Hp = p Frequency-dependent viscous damping is included in the equation of motion relative to the solid phase (Eq. (2)) via the Rayleigh matrix as follows (e.g. Clough and Penzien 2003): where the coefficients a 0 and a 1 can be calculated by selecting one target value of the damping ratio D and two angular frequencies, ω m and ω n , as follows: The viscous damping will be equal to or lower than D for frequencies inside the interval [ω m , ω n ] and higher outside. Typically, the two frequencies ω m and ω n are expressed in Hertz and called f m and f n , respectively.
Park and Hashash (2004) proposed a method to extend the classical two frequencies formulation of Eq. (4) to incorporate any number of frequencies. A frequency-independent viscous damping matrix for one-dimensional site response analysis is suggested by Phillips and Hashash (2009) and, more recently, by Nghiem and Chang (2019).
Depending on how the formulation is implemented in the FE software, the Rayleigh matrix defined by Eq. (3) can be calculated using the initial or the current (tangent) stiffness of the system (e.g. Mazzoni et al. 2006). Hall (2006) and Charney (2008) discussed the potential problems associated with the introduction of some viscous dissipation in dynamic analysis of structures, suggesting to limit and bound the damping forces deriving from the Rayleigh scheme. The following section describes some recent advances for the appropriate use of Rayleigh damping in geotechnical earthquake engineering applications.

Calibration of Rayleigh damping
The Rayleigh damping is defined by selecting the coefficients a 0 and a 1 , which depend on both D and the adopted frequency interval [ω m , ω n ] according to Eq. (4). Woodward and Griffiths (1996) have clearly shown how the amount of viscous damping introduced in the FE dynamic analysis of an earth dam is difficult to quantify a priori but, at the same time, it can significantly control the results of the simulations in terms of permanent displacements induced by the seismic action. Different possible calibration procedures have been proposed in the literature to identify the two controlling frequencies f m and f n corresponding to ω m and ω n . A well-established rule (Hudson et al. 1994) suggests to select f m as the first natural frequency of the deposit f 1 , while f n is assumed equal to n times f m , where n is the closest odd integer larger than the ratio f p /f 1 between the predominant frequency of the input earthquake motion (f p ) and the fundamental frequency of the soil deposit (f 1 ). This choice was based on the elastic analysis of the vibration modes of a shear beam. Kwok et al. (2007) proposed to select, as a first approximation, the first mode of the site and five times this frequency for f m and f n , respectively. Hence, their suggestion to use n equal to 5, as shown in Fig. 1. More recently, Amorosi et al. (2010) discussed a new calibration of the Rayleigh coefficients based on the linear equivalent amplification function and the frequency content of the input motion. In their work, it is suggested to select the first target frequency (f m ) equal to the first site frequency significantly excited by the seismic motion, and to set the second frequency (f n ) equal to the value where the amplification function gets lower than one. Adopting this rule, the overestimation of the peak ground acceleration with respect to the corresponding visco-elastic solution can be considerably reduced.
Even if the frequency interval [ω m , ω n ] can be appropriately defined, the real amount of viscous dissipation introduced in numerical simulations is not really know a priori. It could, nevertheless, be evaluated a posteriori from the examination of the free vibrations of the system after the end of the earthquake, as shown by Amorosi and Elia (2008) and Elia et al. (2011) for the dynamic analysis of an earth dam in Southern Italy. In this case study, a kinematic hardening constitutive model (Kavvadas and Amorosi 2000) was used to simulate the cyclic behaviour of the dam soils. To compensate for the model underestimation of hysteretic dissipation in the small-strain range with respect to the experimental data, a 2% Rayleigh damping was added at the beginning of the FE dynamic simulation. This was calibrated using the initial stiffness of the system and the procedure proposed by Kwok et al. (2007). Then, the evolution of the horizontal displacement recorded at the dam crest in the postseismic stage of the analysis, shown in Fig. 2, was investigated. It should be noted that the damping ratio was evaluated through the logarithmic decrement method (Clough and Penzien 2003), the same method adopted for the determination of damping during resonant column tests, considering the reduction in amplitude of the residual oscillations of the system composed by the dam and its foundation layer. As the amplitude of the observed free vibrations was very small, the computed value of the damping ratio could not be associated with the material (hysteretic) dissipation provided by the constitutive model, but it was representative of the viscous dissipation, introduced through the Rayleigh matrix, acting at the first natural frequency of the system. The damping calculated a posteriori (i.e. 1.7%- Fig. 2) was found to be consistent with the target value of 2% defined a priori, thus showing how it can be possible to control with considerable confidence the amount of viscous dissipation introduced in the FE dynamic simulations. The algebraic counterparts of Eq.
(2) are obtained by applying a time integration scheme. Assuming that the values of displacements, pore pressures and their time derivatives have been obtained at time t n , the integration consists of updating the same quantities at the next time step t n+1 . The integration of the dynamic equations of motion can be performed using time stepping schemes characterised by different accuracy, stability, algorithmic damping and run-time (Bathe and Wilson 1976). A review of time integration schemes for dynamic geotechnical analysis is presented by Kontoe (2006). Adopting, for example, the Generalised Newmark scheme (Katona and Zienkiewicz 1985), the solid phase can be solved as follows: Similarly for the fluid phase: where the coefficients: are typically chosen for unconditional stability of the recurrence scheme (Zienkiewicz et al. 1999). The substitution of the above approximations into Eq. (2) leads to a system of coupled nonlinear equations which are solved iteratively by the FE code using the Newton-Raphson procedure. If the coefficients of Eq. (7) are all set equal to 0.5, the time stepping scheme reduces to the trapezoidal scheme, the algorithm remains implicit, but it does not provide numerical (or algorithmic) damping during the integration of the governing equations. In this case, numerical oscillations may occur during the analysis if no other sources of dissipation (e.g. viscous or hysteretic) are introduced in the simulation (Zienkiewicz et al. 1999). Therefore, coefficient values larger than 0.5 are typically adopted, leading to a numerical dissipation which is more pronounced at high frequencies (Amorosi et al. 2010).

Soil constitutive modelling
The stiffness matrix of the system K in Eq. (2) depends on the constitutive model adopted to describe the soil mechanical behaviour in the numerical simulations. Many nonlinear effective stress constitutive assumptions have been proposed in the past years, all stemming from the general theory of elastoplasticity. A review of advanced soil constitutive models is outside the scope of the paper. In all cases, the implementation of these models into numerical software usually implies two issues: 1) the calibration of several parameters (sometimes even more than ten) against appropriate laboratory and in-situ data; 2) the initialisation of the model state variables accounting for the complex and, at times, unknown past stress history of the material.

Calibration of the model parameters
The higher is the complexity of the soil mechanical behaviour to be modelled (depending on the physical phenomena that needs to be captured) the larger is the number of parameters introduced in the constitutive formulation. This is not necessary an issue, provided that a clear procedure for the calibration of these coefficients can be followed. Most of the parameters have a clear physical meaning, such as the elastic moduli, the critical stress ratio or the slopes of normal compression and swelling lines in the compression plane for the elasto-plastic models derived from the Critical State Soil Mechanics theory (Roscoe and Burland 1968). Others, instead, can be obtained from an optimisation procedure through the comparison with laboratory and in-situ measurements. For dynamic BVPs, cyclic triaxial and simple shear, resonant column/torsional shear (RC/TS) and bender element tests could be adopted for the characterisation of the soil response in terms of small-strain stiffness, normalized shear modulus degradation curve, excess pore water pressure and damping curves with cyclic shear strain. This can be carried out in conjunction with geophysical, cross-hole and down-hole measurements to define the small-strain stiffness profile of the soil deposit in the FE model. In order to calibrate the model parameters against RC/TS results, single element simulations of undrained cyclic simple shear (CSS) tests can be performed, imposing different shear strain amplitudes and assessing the secant shear modulus (G sec ) and the damping ratio for each amplitude after a number of cycles sufficient to reach steady-state condition (e.g. Elia et al. 2011;Seidalinov and Taiebat 2014;Elia and Rouainia 2016). Specifically, G sec is equal to the ratio between the cyclic shear stress (τ c ) and the cyclic shear strain (γ c ), while the damping ratio can be evaluated as: where W D is the energy dissipated during the cycle (i.e. the area of the hysteresis loop A loop ) and W S is the maximum strain energy stored during the loading phase (approximated to the area ½τ c γ c ). Figure 3a shows the typical stress-strain curves obtained with a kinematic hardening soil model imposing increasing cyclic shear strains during CSS single element simulations. Referring, in particular, to the CSS test carried out with a γ c equal to 1%, the secant shear modulus is represented by the slope of the line connecting the tips of the last cycle, while the areas needed for the calculation of the damping ratio according to Eq. (8) are highlighted with shadings. Once the normalised shear modulus and the hysteretic dissipation are calculated for different imposed shear strain amplitudes, the G/G 0 and D curves predicted by the model can be plotted point by point, as shown in Fig. 3b.
An example of a real calibration is presented in Figs. 4 and 5, where the comparison between the predictions of a kinematic hardening soil model (i.e. the RMW model proposed by Rouainia and Muir Wood (2000)) and the laboratory data from undrained triaxial compression tests and combined resonant column/torsional shear tests performed on natural Avezzano Clay samples is shown, respectively (Cabangon et al. 2019). It should be noted that only 2% Rayleigh damping has been added to avoid the propagation of spurious high frequencies and to compensate for the model underestimation of damping in the small-strain range (Fig. 5). At the same time, the hysteretic dissipation provided by the kinematic hardening formulation at large strains overestimates the experimental data, as typically observed when this class of constitutive models is adopted (e.g. Amorosi and Elia 2008;Elia et al. 2011;Elia and Rouainia 2016). To address this point, some adjustments in the plastic modulus law have been proposed (e.g. Seidalinov and Taiebat 2014), but very recently Elia et al. (2021) have demonstrated that the largest loops obtained from the back-analysis of ground response FE results, which can potentially provide a high source of dissipation, are statistically irrelevant, as their frequency of occurrence during the earthquake time history is very low.
Moreover, the irregular nature of the earthquake motion implies that the stress-strain cycles induced in the FE soil deposit by the seismic event are typically non-symmetric. Figure 6 shows the predictions obtained at the local level (i.e. at a Gauss point) in terms of shear strain time history and corresponding stress-strain curve during nonlinear FE dynamic simulations of the free-field response of a soil deposit using the RMW model (Guzel 2018). The maximum value of the shear strain induced by the earthquake applied at bedrock (Fig. 6a) is reported in Fig. 6b and in Fig. 6c as the distance between points A and B. At the beginning of the analysis, when the earthquake has not attained its maximum acceleration values yet, the cycles are centred on the origin (see the red part of the curves). This is followed by larger loops associated to the higher accelerations (highlighted in green) and, then, the shear strains oscillate around a final non zero value at the end of the motion (blue part of the curves). Differently from visco-elastic dynamic results, a continuous change of shear stiffness and hysteretic damping can be observed throughout the earthquake excitation and most of the loops induced by the seismic action are not centred around the origin, due to plastic strain accumulation (Fig. 6c). Therefore, the maximum value attained by the shear strain at each depth may not be representative of the actual strain level induced by the earthquake in the soil deposit, as ratcheting effects should be taken into account (Guzel 2018). From the experimental point of view, further laboratory testing is needed to better understand the asymmetrical cyclic behaviour of soils at large strains (i.e. for γ higher than 1.0%) and subjected to multi-directional loading conditions (e.g. Yang et al. 2018; Sun and Biscontin 2019). This could guide further developments of advanced soil constitutive models.

Initialisation of the model state variables
When dynamic analyses of geotechnical BVPs are performed, the application of the seismic loading is usually preceded by the simulation of the geostatic initial stress condition. Nevertheless, the at-rest coefficient of earth pressure predicted by an advanced elasto-plastic model is not known in advance. Therefore, consistently with the model parameters calibration described before, a simplified, though realistic simulation of the geological and past stress history of the soil deposit and geotechnical structure under investigation should be conducted prior to any dynamic analysis. This is the case, for example, presented by Elia et al. (2011) in the evaluation of the seismic performance of the Marana Capacciotti earth dam. Before the application of the seismic loading, the initial values of the internal variables of the adopted constitutive model (i.e. initial stress state and hardening variables) have been determined from a sequence of static FE analyses reproducing the geological history of the deposit, the dam construction phases and the subsequent first reservoir impounding. At the end of the static simulations, the stress state and the values of the internal variables have been checked for consistency with those assumed in the model calibration phase. The values of the overconsolidation ratio, R, along the dam axis were in good agreement with those assumed during the model calibration (Fig. 7a). The computed values of the initial shear modulus, G 0 , were also found to agree well with the bender element results (Fig. 7b), thus confirming the consistency between laboratory measurements, model state variables values assumed in the calibration stage and those obtained at the end of the static FE analyses.

Finite element modelling
The accuracy of the FE solution of the system of Eq.
(2) is greatly controlled by the adopted time and spatial discretisation (i.e. the value of Δt used in the time stepping scheme and the dimension and type of elements, respectively). In addition, the definition of appropriate boundary conditions, including the seismic input motion applied at bedrock, is required.

Model discretisation and boundary conditions
When performing numerical simulations of BVPs, infinite continuous soil deposits are reduced to finite discretised domains with the adoption of some constraints along their boundaries, which allow to artificially simulate the far-field condition. Users are, therefore, required to define the appropriate extension of the numerical model, the geometrical dimensions of the finite elements and the boundary conditions to be applied. While all these aspects are well understood in the context of static simulations, the literature concerning numerical analyses in dynamic conditions is less exhaustive. The characteristic dimension of the finite elements adopted for dynamic analyses must satisfy the condition that the spacing between two consecutive nodes, Δl node , has to be smaller than approximately one-tenth to one-eighth of the wavelength associated with the maximum frequency component, f max , of the input wave (Kuhlemeyer and Lysmer 1973): where V S,min is the lowest shear wave velocity in the soil deposit. Equation (9) has been proposed for the propagation of a single, harmonic, one-dimensional wave into an elastic material. A finer discretisation should be adopted when the mechanical behaviour is expected to be nonlinear, due to plasticity.
To avoid spurious wave reflections at the vertical boundaries of symmetric numerical models, "tiednodes" can be employed in geotechnical dynamic simulations. Their effectiveness in absorbing the energy induced by the seismic action has been demonstrated by Zienkiewicz et al. (1999). For non-symmetric 2D and 3D problems, the hypothesis of tied (9) Δl node ≤ min ∕(8 ÷ 10) = V S,min ∕(8 ÷ 10)f max displacements along the lateral boundaries needs to be abandoned. Possible alternatives are represented by viscous (Lysmer and Kuhlemeyer 1969) and freefield boundaries or by the Domain Reduction Method (DRM) (Bielak et al. 2003;Yoshimura et al. 2003;Kontoe et al. 2008Kontoe et al. , 2009. Each technique used to define the boundary conditions has its own advantages and limitations (e.g. Roesset and Ettouney 1977), but not all of them are commonly implemented in FE/FD codes (as the DRM). In 2D space, the standard viscous boundary is mechanically equivalent to a system of dashpots oriented in normal and tangential directions of the mesh boundary, such that: where v x , v y , σ xx and τ are the normal and tangential node velocities and stresses, respectively, ρ is the mass density of the soil and V P and V S are the P-and S-wave velocities, respectively. This option is generally used for problems where the dynamic source is applied inside the mesh (e.g. vibrations induced by a generator or by pile driving). When adopting viscous boundaries, questions concerning the appropriate lateral extension of the mesh arise, as they may not be able to totally absorb the incoming seismic wave, especially for low angles of incidence (Lysmer and Kuhlemeyer 1969). For instance, Amorosi et al. (2010) have shown that the use of viscous boundaries in 2D ground response analyses requires a mesh characterised by a width to height ratio between five and eight to avoid spurious reflections in the centre of the FE model. As an alternative to employing dampers, the user can set a high value of Rayleigh damping (e.g. 25 ÷ 30%) to the finite elements closer to the vertical mesh boundaries. Again, Elia et al. (2011) have demonstrated the efficacy of this solution in imposing the far-field response when the model boundaries are sufficiently distant from the centre of the mesh.
A free-field boundary consists of a one-dimensional soil column (in 2D problems) coupled to the main domain by viscous dampers. The free-field motion computed for the soil column is transferred to the numerical model by applying along the boundary the equivalent normal and tangential stresses. These are obtained from Eq. (10) using the relative particle velocities between the main grid and the column element at (10) xx = V P v x xy = = V S v y the same depth. This option is more suitable for seismic analysis and has been validated for 2D BVPs (e.g. Cundall et al. 1980). The effectiveness of viscous and free-field boundary conditions in the dynamic analysis of 3D cases needs, instead, to be investigated further.
The dynamic input can be applied directly as a prescribed displacement or velocity time history at the base of the numerical model if a rigid bedrock condition is assumed. In this case, the bottom side of the mesh will be a non-absorbing boundary, i.e. it will reflect back into the model any outgoing wave. In contrast to a rigid base hypothesis, compliant base boundaries (Joyner and Chen 1975) are nowadays available in some commercial software to simulate a deformable (i.e. dissipating) bedrock, when the relative stiffness between the bedrock and the soil deposit (i.e. the impedance) is not infinite. With this option, the velocity record of input motion is transformed into a shear stress time history and applied to a nonreflecting (i.e. viscous) boundary as follows: where ρ bedrock and V S,bedrock are the mass density and the shear wave velocity of the bedrock, respectively, and the term ½v outcrop is one half of the velocity time history recoded at the outcrop of the bedrock. In this way, half of the input is absorbed by the viscous dashpots located along the base and half is transferred to the main domain. Equation (11) does not depend on the bedrock depth or on the thickness of the bedrock implemented in the numerical model, but it is based on the assumption of an elastic behaviour of the bedrock. The compliant base hypothesis makes the accelerograms at the bottom of the numerical model an output of the simulation, rather than an input data. This allows to avoid a specific deconvolution procedure of the input motion from the outcrop to the base of the FE model. This boundary condition can successfully be employed to simulate complex geometries of deformable bedrocks, as shown, among others, by Falcone et al. (2018) for 3D ground response analyses aimed at microzonation studies.

Some insights on the input motion selection and scaling strategies
The correct definition of the design seismic actions, based on seismic hazard studies, is essential to fully implement a performance-based design approach proposed by code prescriptions (e.g. Eurocode 8). The majority of the research on the appropriate methods to select and scale ground motions for linear and nonlinear time history analyses has been devoted, in the last years, to structural dynamics problems more than to ground response (e.g. NIST 2011). In the earthquake geotechnical engineering field, the seismic input signals are typically scaled to the PGA (Peak Ground Acceleration) of the specific site, representing the maximum acceleration expected at the bedrock outcropping surface. Nevertheless, different ground motion modification methods might be used to adapt the selected acceleration time histories, in order to minimise the bias and reduce the number of simulations needed to obtain statistically stable and robust results. In fact, the input motion could be linearly scaled by using a suitable scale factor without altering the shape of the response spectrum. The spectral shape of the input motion can also be modified by adding wavelets to match a target spectrum. Whereas the adoption of such scaling/matching techniques is common in the structural engineering literature (Shome et al. 1998;Hancock et al. 2006;Haselton 2009;Galasso 2010;Michaud and Léger 2014), their application to geotechnical earthquake engineering problems is still limited (Amirzehni et al. 2015;Guzel et al. 2017;Guzel 2018;Elia et al. 2019). As an example, Elia et al. (2019) investigated the effect of different earthquake scaling/matching strategies on the nonlinear dynamic response of an anchored diaphragm wall supporting a deep excavation. Five scaling/matching methods are chosen to linearly or spectrally modify seven input motions, selected on the basis of the seismic hazard analysis of the site, specifically: PGA, Sa(T 1 ), ASCE and MSE scaling techniques and the Spectral Matching (SM) method. Two seismic intensity levels were considered. The results of the nonlinear dynamic simulations are illustrated in Fig. 8 in terms of mean and standard deviation of the maximum horizontal displacement experienced by the wall during each set of input motions. The statistical analysis indicated that, for low intensity input motions, the Sa(T 1 ) scaling and spectral matching methods are characterised by the lowest mean and standard deviation of the output, while higher variability is obtained using the MSE, ASCE and PGA scaling techniques. When higher seismic intensity levels are considered, the best performance is provided by the MSE scaling and SM techniques. Simply scaling the ground motion records at the same PGA does not return results that are statistically matched with those obtained with more advanced scaling strategies. The PGA is not a good indicator of the frequency content of the ground motion and PGA scaling should be avoided in nonlinear dynamic analyses (e.g. Amirzehni et al. 2015;Elia et al. 2019). The Sa(T 1 ) scaling technique does return good statistical results when a moderate elongation of the fundamental period of the system is induced by soil nonlinearity. The variability of the dynamic response is, instead, significantly reduced when the selection strategy seeks the compatibility with the target response spectrum over a range of periods (e.g. MSE scaling), especially for high seismic intensity levels. Although the full spectral matching strategy implies a fictitious modification of the frequency content of natural records (and, therefore, it is not always accepted in the national and international code provisions), statistically robust output can be obtained by using spectrally matched accelerograms (e.g. Amirzehni et al. 2015;Elia et al. 2019).

Conclusions
Important advances have been made during recent years in the field of geotechnical earthquake engineering and soil dynamics modelling, although some challenges remain to be addressed. The paper has shown, with particular reference to advanced numerical methods for the integration of nonlinear constitutive equations in seismic design of geotechnical structures, that: -the effect of the Rayleigh damping on FE dynamic results can be fully regulated by an appropriate selection of the controlling frequencies and by the evaluation of its real amount in the post-seismic stage of the simulations; -the user of a FE code should be able to fully control and distinguish between the different sources of dissipation which can be introduced in the simulation (i.e. viscous, numerical and hysteretic damping); -the main dissipation of the energy introduced by the seismic action should be guaranteed by the material (hysteretic) damping provided by the soil constitutive model employed in the numerical simulations, thus relegating the Rayleigh and numerical damping to a small contribution; -the adoption of effective stress based soil constitutive laws, able to capture the material cyclic/ dynamic behaviour observed in-situ and in the laboratory, should be promoted in the engineering practice by a clear calibration and initialisation strategy of their parameters and state variables; -the effectiveness of the lateral and bottom boundary conditions should always be tested in detail, especially when 3D dynamic analyses are performed, through a direct comparison with the freefield solution; -the influence of the scaling/matching techniques of the input motions on the results of nonlinear dynamic analyses of geotechnical BVPs should be investigated further in order to reduce the number of simulations needed to obtain statistically stable and robust outputs.
In conclusion, the paper has presented a review of the key ingredients governing the predictive capabilities of nonlinear FE schemes implementing sophisticated elasto-plastic soil models. It is hoped that this will provide a useful framework for their implementation by non-expert users, as an on-going transfer of research knowledge into practice represents a significant requisite to support performance-based earthquake engineering and quantitative seismic risk assessments.
Author Contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by both authors. The first draft of the manuscript was written by Gaetano Elia and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Data Availability
The datasets generated during and/or analysed during the current study are available in previously published papers by the authors.

Declarations
Competing Interests The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.