Discretization effects in the finite element simulation of seismic waves in elastic and elastic-plastic media

Presented here is a numerical investigation that (re-)appraises standard rules for space/time discretization in seismic wave propagation analyses. Although the issue is almost off the table of research, situations are often encountered where (established) discretization criteria are not observed and inaccurate results possibly obtained. In particular, a detailed analysis of discretization criteria is carried out for wave propagation through elastic and elastic-plastic media. The establishment of such criteria is especially important when accurate prediction of high-frequency motion is needed and/or in the presence of highly non-linear material models. Current discretization rules for wave problems in solids are critically assessed as a conditio sine qua non for improving verification/validation procedures in applied seismology and earthquake engineering. For this purpose, the propagation of shear waves through a 1D stack of 3D finite elements is considered, including the use of wide-band input motions in combination with both linear elastic and non-linear elastic-plastic material models. The blind use of usual rules of thumb is shown to be sometimes debatable, and an effort is made to provide improved discretization criteria. Possible pitfalls of wave simulations are pointed out by highlighting the dependence of discretization effects on time duration, spatial location, material model and specific output variable considered.


Introduction
The study of wave motion is of utmost importance in many applied sciences, as it supports the understanding of transient phenomena in many natural and anthropic dynamic systems. In particular, seismic waves propagating through the earth crust deserve the highest consideration, especially in light of their destructive potential and socio-economical impact.
In the last decades, mathematicians, geophysicists and engineers have devoted massive research efforts to the prediction of seismic motion, based on either analytical [21,[32][33][34]39] or numerical methods [2,53,63]. When linear elastic wave problems are considered, either time-domain or frequency-domain solutions may be sought, whereas time-domain approaches are usually needed in the presence of non-linearities (constitutive or geometrical). In this respect, it should be remarked that much interest in earthquake engineering is nowadays on non-linear wave phenomena, since they govern (i) the occurrence of natural catastrophes (e.g., landslides and debris flows) induced by soil instabilities, such as liquefaction and strain localization [18,24,63]; (ii) the interaction between geomaterials and man-made structures [13,16,20,28,53,59].
Abstract Presented here is a numerical investigation that (re-) appraises standard rules for space/time discretization in seismic wave propagation analyses. Although the issue is almost off the table of research, situations are often encountered where (established) discretization criteria are not observed and inaccurate results possibly obtained. In particular, a detailed analysis of discretization criteria is carried out for wave propagation through elastic and elastic-plastic media. The establishment of such criteria is especially important when accurate prediction of high-frequency motion is needed and/or in the presence of highly non-linear material models. Current discretization rules for wave problems in solids are critically assessed as a conditio sine qua non for improving verification/ validation procedures in applied seismology and earthquake engineering. For this purpose, the propagation of shear waves through a 1D stack of 3D finite elements is considered, including the use of wide-band input motions in combination with 1 3 It is thus apparent that reliable numerical simulations of seismic motion and earthquake-soil-structure interaction can only be performed by means of high-fidelity computational tools, capable of coping with the remarkable complexity of the aforementioned problems. The accuracy of numerical predictions is in turn affected by, at least, the following four factors: 1. selection of the numerical solution algorithm; 2. mathematical description of material behavior (constitutive model); 3. computer implementation; 4. set-up of the computational discrete model.
The assessment of the above four items is the main core of a thorough verification and validation process [3,45,51]: is the mathematical problem numerically solved to the desired degree of accuracy? Do numerical results reasonably reproduce real world phenomena?
The present work focuses on the fourth item in the list, and specifically on the selection of appropriate time-step and element size in dynamic Finite Element (FE) computations. This problem seems to have been solved quite long ago in the form of "rules of thumb" for space/time discretization [38,41], so that not many works on the subject have been published ever since [4,5,14,55]. Furthermore, the relationship between discretization and accuracy in wave simulations has been mainly investigated for linear elastic problems.
In light of the above premises, the authors aim an upto-date contribution to the matter, also accounting for the large importance assumed in recent years by non-linear, elastic-plastic wave problems. The key features of the present work are hereafter summarized: -only 1D shear wave propagation tests are performed for a more straightforward interpretation of numerical results; -discretization effects have been illustrated in both time and frequency domains, and then quantified via modern misfit criteria formulated in the full time-frequency domain; -since discretization effects depend in general on the numerical algorithm adopted, a widespread FE approximation scheme has been here adopted; -the role of constitutive non-linearity (plasticity) is discussed; -the whole study should be regarded as a numerical "falsification test" for the "rules of thumb" previously mentioned [38,41].
The ultimate goal of this work is to reopen the debate on the accuracy of wave simulations from a verification/ validation perspective, also in the presence of constitutive non-linearities. The results reported provide renovated critical insight into, and review of, traditional discretization rules for practical simulation purposes.

FE modeling of 1D seismic wave propagation
1D shear wave problems originate from the ideal situation in which wave propagation is nearly vertical, with no lateral geometrical/material inhomogeneities. In these conditions, all vertical cross-section can be regarded as symmetry planes and the soil deposit undergoes a "double plane-strain" deformation, with both horizontal direct strains prevented by symmetry [10,49]. As a consequence, all variables only depend on time and vertical elevation (the problem is geometrically one-dimensional), whereas the stress state is still multi-axial [17]. The initial-boundary value problem under consideration is sketched in Fig. 1. Like in general 3D problems, the numerical analysis of 1D seismic wave propagation requires a suitable computational platform for (i) space/time discretization, (ii) material modeling and (iii) simulation under given initial/ boundary conditions. The Real ESSI Simulator has been used here for these purposes.

Space discretization and time marching
The Real ESSI program is based on a standard displacement FE formulation, where displacement components are taken as unknown variables in the numerical approximation [62]. As for space discretization, the 1D FE model has been built using a stack of properly constrained 3D brick elements-as was previously done, for instance, by [10]. Real ESSI program enables the use of 8-, 20-and 27-node elements, so that several options are given in terms of spatial interpolation order. The well-known Newmark method has been adopted for time marching [43]. The main feature of the integration algorithm relates to the approximate series expansion for displacement and velocity components, u and u respectively: between two subsequent time-steps n and n + 1. Importantly, the expansion uses two parameters, β and γ, governing the accuracy and stability properties of the algorithm. It is worth reminding that the algorithm is unconditionally stable as long as [23]: γ = 1/2 is required for second-order accuracy, whereas any γ value larger than 1 / 2 introduces numerical attenuation (damping). In this study, the pair γ = 1/2 and β = 1/4 (no algorithmic dissipation) is exclusively considered.

Material modeling
The Real ESSI program provides a number of material modeling options, ranging from simple linear-elastic to advanced elastic-plastic constitutive relationships for cyclically loaded soils [18,63]. Hereafter, the material models adopted for wave propagation analyses are briefly described, namely (i) the standard linear elastic material model, (ii) the elastic-plastic von Mises model with linear kinematic hardening [26,40] and (iii) the bounding surface elastic-plastic model by [48].

Linear elastic model
Discretization issues will be first addressed with reference to linear elastic problems. While relevant concepts in elastodynamics can be found in [21], it is only worth reminding here the relationship between the shear wave velocity V s and the two elastic parameters (Young's modulus E and Poisson's ratio ν): where ρ is the soil mass density and G = E/[2(1 + ν)] the elastic shear modulus.

Elastic-plastic: von Mises kinematic hardening (VMKH) model
The relationship among discretization, accuracy and material non-linearity will be first explored through the elasticplastic von Mises kinematic hardening (VMKH) model, of the same kind described in [26,40]. The VMKH model is very well-known in literature and widely employed for cyclically loaded metals, while the application to soil dynamics is limited to undrained loading conditions in combination with total stress analysis [44,63]. Although the assumption of linear hardening is not the most accurate for soils 1 , it has been here introduced for numerical convenience. In fact, owing to linear hardening, the post-yielding stiffness is constant, not strain-dependent: this implies an unrealistic unbounded strength, but allows to identify the elastic-plastic shear stiffness with no ambiguity. Only four constitutive parameters need to be set: -two elastic parameters-E and ν; -one yielding parameter-k-proportional to the initial size of the cylindrical yield locus in the stress space; -one hardening parameter-h-governing the post-yielding (elastic-plastic) stiffness.

Elastic-plastic: Pisanò bounding surface (PBS) model
The more sophisticated constitutive relationship recently proposed by [48] will be also used. At variance with the aforementioned VMKH formulation, the Pisanò bounding surface (PBS) model can quite accurately reproduce important aspects of monotonic/cyclic soil behavior within the effective stress framework [35,44,60], such as: 1. development of inelastic strains from the very onset of loading. This is reproduced by exploiting the concept of "vanishing yield locus"; 2. frictional shear strength, i.e., depending on the effective confining pressure; 3. non-linear hardening, implying a continuous transition from small-strain to failure stiffness; 4. coupling between deviatoric and volumetric responses; 5. stiffness degradation and damping under cyclic shear loading.
A remarkable feature of the PBS constitutive formulation is the low number of input parameters required (only seven), which makes the model particularly suitable for practical use: -two elastic parameters-E and ν-to characterize the material behavior at vanishing strains; -one shear strength parameter-M-directly related to the material frictional angle; -two parameters-k d and ξ-governing the development of plastic volumetric strains during shearing; -two hardening parameters-h and m-to be identified on the basis of stiffness degradation and damping cyclic curves.
Interested readers are addressed to [48] for details about formulation, performance and calibration of the PBS model.

Initial/boundary conditions and input motion
All the FE results hereafter presented have been obtained under the following initial and boundary conditions ( Fig. 1): 1. the system is initially at rest (nil initial velocities and accelerations); 2. a horizontal x-displacement time history is imposed at the bottom boundary to reproduce rigid bedrock conditions; 3. no loads are applied to the top boundary (free surface); 4. the aforementioned "double plane-strain" conditions has been achieved by preventing y-displacements throughout the model, as well as imposing master/slave connections to nodes at the same elevation (tied nodes).
As for the input displacement, the Ormsby wavelet [52] fits the authors' intent: where t denotes the physical time and A the signal amplitude, sinc(x) = (sin x)/x is the cardinal sine function, f i (i = 1, 2, 3, 4) stand for the low-cut, low-pass, high-cut and high-pass frequencies, respectively. The meaning of the f i frequencies can be grasped from Fig. 2b, illustrating the amplitude Fourier spectrum of function (5). In particular, the suitability of the Ormsby wavelet has a twofold motivation: 1. function (5) has a number of sign reversals and will induce several loading/unloading cycles into the soil undergoing wave motion (Fig. 2a); 2. the peculiar flat branch in the amplitude Fourier spectrum ( Fig. 2b) is convenient for frequency domain analysis (see next section).
The above features of the Ormsby wavelet will enable the analysis of discretization effects over frequency ranges of choice. Although most seismic energy relates to frequencies lower than 20 Hz, ensuring accuracy at higher frequencies may be relevant when seismic serviceability analyses are to be performed for structures, systems and components (SSCs) related to nuclear power plants and other industrial objects.

Misfit criteria
The analysis of discretization effects requires objective criteria to quantify the discrepancy (misfit) between different numerical solutions. In numerical seismology, the difference seismogram between the numerical solution and a reliable reference solution is often adopted for this purpose, although it only enables visual/qualitative observations; simple integral criteria (e.g., root mean square misfit) can provide some quantitative insight, but still with no distinction of amplitude or phase errors. A significant improvement in this area was introduced by [36], who compared seismograms on the basis of the time-frequency representation (TFR) obtained through continuous wavelet transformation [22]. The TFR of signal misfit allows to extract the time evolution of the spectral content, and thus to define the following local time-frequency envelope difference: and time-frequency phase difference: where W (t, f ) and W REF (t, f ) are the TFR (wavelet transform) of the signal "under evaluation" and the reference seismogram, respectively. As explained by [36], it is also possible to obtain purely time-or frequency-dependent misfit measures by projecting E and P onto one of the two domains. In particular, the following single-values measures for envelope misfit (EM) and phase misfit (PM) may be employed to separate amplitude and phase errors when comparing different signal couples. It should be recalled that the envelope function of an oscillating signal is the smooth curve outlining its extremes, and therefore, carries more information than single amplitude values at given time. While the theoretical background for the above misfit criteria is widely described by [36,37], open-source routines for misfit analysis are available at http://www. nuquake.eu/ComputerCodes/ (TF-MISFITS package). Discretization effects in wave propagation simulations will be assessed in the following on the basis of EM and PM criteria, as previously done by a number of authors [6,19,31,42,47].

Linear elastic wave simulations
In this section, the influence of discretization on accuracy is first discussed for linear elastic problems.

Standard rules for space/time discretization
The selection of appropriate grid spacing 2 and time-step size is usually based on very simple rules. As for space discretization, [41] stated that "the accuracy of the finite element method depends on the ratio obtained by dividing the length of the side of the largest element by the minimum wavelength of elastic waves propagating in the system. For accurate results this ratio should be smaller than 1/12". Since then, it has been believed that approximately ten nodes per wavelength are appropriate in most cases, whereas fewer than ten nodes are likely to result in undesired numerical attenuation/dispersion. Accordingly, suitable maximum grid spacing is usually determined by considering the minimum relevant wavelength (or highest frequency f max ) in the input signal [28]: On the other side, the time-step size also needs to be limited to ensure accuracy and stability [2]. In principle, the smallest fundamental period of the system should be represented with about ten time-steps-same as for space discretization. However, t is often selected on the basis of a different physical argument, i.e., to avoid that a given wave front reaches two consecutive nodes at the same time (this would happen for too large t values): Condition (11) ensures algorithmic stability in many explicit schemes for hyperbolic differential problems [50], and is also often regarded as an accuracy criterion for implicit (unconditionally stable) time marching as well (see Sect. 2.1).

Model parameters
The geometrical/mechanical parameters adopted for elastic wave simulations are here reported. A uniform soil layer has been considered, having thickness H = 1 km and made of an elastic material with ρ = 2000 kg/m 3 , V s = 1000 m/s and ν = 0.3 (corresponding to G = 2 GPa). No Rayleigh damping has been introduced. 2 Henceforth, x will always denote the vertical node spacing, coinciding with the element thickness in the case of 8-node bricks.
As for the input motion, two different Ormsby wavelets have been employed, corresponding with the following input parameters in Eq. (5): Hz; -the amplitude parameter A has been always set to produce at the bottom of the layer a maximum displacement of 1 mm.
As previously mentioned (Sect. 2.3), both inputs 1 and 2 have been used to explore the interplay of discretization effects and input bandwidth.

Discussion of numerical results
The influence of grid spacing and time-step size are discussed separately for the sake of clarity. Since the Real ESSI program is based on a displacement FE formulation, displacement components are the most reliable output; however, some attention is also paid to accelerations, postcalculated through second-order central differentiation. Table 1 provides a list of the comparative simulations performed for linear problems. Each case is denoted by: (5)); (ii) grid spacing x std and (iii) time-step size t std from standard discretization rules (10)-(11); (iv) x and (v) t actually used; (vi) type of brick elements adopted.
The results being presented aim to assess the quality of standard discretization rules, as well as the improvements attainable through refined discretization. For this purpose, the numerical results are discussed in both time and frequency domains-the Fourier spectra of considered time histories are plotted in terms of (i) amplitude and (ii) phase difference with respect to the analytical solution (known at the free surface). Additional quantitative insight is also gained through the EM and PM misfit criteria introduced in Sect. 2.4. Unless differently stated, numerical outputs at the top of the soil layer are considered.  Fig. 4, displacement time histories are not compared with the input motion (as done in Fig. 3a) for the sake of brevity, whereas only a reduced time window around the output motion is displayed for clearer visualisation (e.g., as in Fig. 3b) Figs. 3, 4, 5, 6 suggest the following observations (some of which expected):

Influence of grid spacing
-even though x std is set on the basis of the maximum frequency f max , its suitability is not uniform over the input spectrum. Indeed, increasing inaccuracies in the frequency domain are clearly visible as f max is approached (check for instance the Fourier amplitudes compared in Figs. 3c and 4, 5, 6b). Grid spacing affects output Fourier spectra both in amplitude and phase; -in all cases, envelope and phase misfits, EM and PM, are quantitatively very similar (Figs. 3e and 4, 5, 6d); -reducing x below x std is beneficial only if t is also lower than t std . This is apparent in Fig. 3e, where an increase in EM and PM is observed as x gets lower than x std . Conversely, monotonic EM/PM trends are shown in Figs. 4, 5d; -at given grid spacing x, reducing the time-step improves the numerical solution mostly in terms of Fourier phase, not amplitude (compares Figs. 3c-d, 4b-c). It may be generally stated that, when x is not appropriate, reducing the time-step size does not produce substantial improvements; -based on these initial examples, a grid spacing x in the order of V s /20f max = �x std /2 ensures high accuracy (EM and PM < 10 % ) in combination with �t = �x/2V s = �t std /2. These enhanced discretization rules hold for low-order FEs (8-node brick elements) but are not affected by the frequency bandwidth of the input signal. In the latter respect, Figs. 4, 5d show quantitatively similar EM-PM trends for f max equal to 20 Hz and 50 Hz. Also, minimum misfits are attained in the EL2 case (Fig. 4d), where a smaller �t/�t std ratio has been purposely set.
The above conclusions apply to 8-node brick elements, while Fig. 6 shows that "ten elements per wavelength" are still suitable when higher-order elements (here 27-node bricks 3 ) are employed. However, this lighter requirement for grid spacing seems to perform well in combination with �t ≤ �x/2V s , and results in EM and PM lower than 10 % even for �x/�x std = 2.
It is also important to evaluate grid spacing effects on acceleration components, as they will affect the inertial forces transmitted to man-made structures on the ground surface. Since acceleration time histories are dominated by high frequencies, the poorer performance of standard discretization rules at high frequencies becomes more evident. In Figs. 7 and 8, grid spacing plays qualitatively as in Figs. 3, 4, 5, although the EM/PM trends-similar in shape-are shifted upwards. This means that, in the presence of low-order elements, more severe discretization requirements should be fulfilled if very accurate accelerations are needed.

Influence of time-step size
For given grid spacings, the influence of t has been studied by varying the time-step size with respect to the limit -when 27-node bricks are used, the use of x = x std and �t ≤ �t std /2 is still an appropriate option, giving rise to EM and PM lower than 5 % (Fig. 12). Even in this case, discretization errors are still governed by phase differences, while excellent performance in terms of Fourier amplitude is observed; -Figs. 13 and 14 show that the above findings apply qualitative to acceleration time histories as well. However, EM and PM values are quite high (significantly larger than 10 %) when t ≥ t std , regardless of the grid spacing ratio. Accuracy is quickly regained when t is reduced and �x < �x std /2.
While the above conclusions have been all drawn on the basis of the first incoming wave, many reflected waves may in reality hit the ground surface because of soil layering. In the present elastic case (no energy dissipation), perfect reflections occur at the lower rigid bedrock and never-ending wave motion is established. It is thus interesting to check how discretization errors propagate in time at the free surface, as is shown in Fig. 15. Subsequent wave arrivals are compared in the time (Fig. 15a, b) and frequency (Fig. 15c, d) domains, where a gradual "accumulation" of wave dispersion can be observed. Even though satisfactory accuracy is achieved on the first arrival, an increase in high-frequency phase difference is detected in Fig. 15d, with negligible variation in Fourier amplitude (Fig. 15c). Cumulative wave dispersion implies that selecting suitable x and t becomes increasingly delicate for large FE models and long durations.

Non-linear elastic-plastic wave simulations
This section concerns discretization effects in presence of material non-linearity. As most commonly done in  Geomechanics [63], the non-linear cyclic response of geomaterials can be described in the framework of elastoplasticity, and here the VMKH and PM models described in Sect. 2.2 have been adopted. Prior to presenting numerical results, some preliminary remarks should be made: -the non-linear problem under consideration cannot be solved analytically. Therefore, the quality of discretization settings may only be assessed by evaluating the converging behavior of numerical solutions upon x-t refinement; -with no analytical solution at hand, one needs engineering judgement to establish when the (unknown) exact solution is reasonably approached. In this respect, light is shed on several expected pitfalls, all relevant to the global verification process [3,45,51]; -the accuracy of non-linear computations is highly affected by the input amplitude. This governs the amount of non-linearity mobilized by wave motion and, as a consequence, the accuracy of numerical solutions at varying discretization.
In non-linear (elastic-plastic) problems, discretization is not only responsible for the numerical representation of waves (dissipation, dispersion, stability), but also governs the accuracy of constitutive integration [8,54]. For instance, changes in time-step size will affect the strain size driving the constitutive integration algorithm and, in turn, the final simulation results. This dependence of the constitutive response (material model and constitutive integration algorithm) on the dynamic step size precludes direct development of automatic criteria for discretization. However, as tangent elastic-plastic response can be established for any stress-strain combination, (lowest) elastic-plastic (shear) stiffness may be used to develop suitable discretization via Equation 4. Apparently, this approach assumes that the stress-strain response is already known, as is not the case when discretization is being set. This means that an iterative approach is in principle needed, whereby one will first design discretization based on an estimate of the strain level, run the dynamic simulation, and record the actual stress-strain response. After few iterations, a stable discretization will be usually achieved. In this study, VMKH and PM constitutive equations have been integrated via the standard forward Euler, explicit algorithm [11,15]. Although implicit algorithms may possess better accuracy/stability properties, explicit integration is often preferred for advanced constitutive formulations and cyclic loading [27]. There is also wide consensus on the poor performance in elastic-plastic computations of time-step sizes derived through elastic parameters and Equation (11), especially in combination with explicit stress-point algorithms. For this reason, the following time marching rule may be regarded as an upper bound for nonlinear problems (instead of (11)): In the following, rules (10) and (12) will be assumed as starting discretization criteria and critically assessed. For shorter discussion, only input 1 ( f max = 20 Hz) and 8-node brick elements are employed for non-linear simulations.

Model parameters and parametric analysis
A heterogeneous 1 km thick soil deposit has been considered, formed by a 200 m thick VMKH sub-layer resting on an elastic stratum (remaining 800 m). At the surface, a thin layer (5 m) of elastic material has been added to prevent numerical problems with very strong motions and the socalled whip effect. The following constitutive parameters (see Sect. 2.2.2) have been set (same elastic parameters for both the VMKH and the elastic sub-layers), with no algorithmic nor Rayleigh damping introduced in numerical computations.: In the analysis of VMKH cases, the influence of the hardening parameter (h) and the input amplitude (A) has been also considered, as they affect the material elastic-plastic stiffness and the amount of plasticity mobilized. The VMKH simulation programme is reported in Table 2, where t std has been determined through Equation (12) (i.e., �t std = �x/10V s ).

Influence of grid spacing and time-step size
The results in Figs. 16 and 17 exemplify the role played by space discretization in elastic-plastic simulations. These results have been obtained by employing a timestep smaller than t std (cases VMKH1-2 in Table 2), a low input amplitude (A = 0.1 mm corresponds with a peak ground acceleration approximately equal to 0.05g), and two different values of the hardening parameter (h = 0.5E and h = 0.05E). The following observations arise from the two figures: -propagation through a dissipative elastic-plastic material alters significantly the shape of the input signal. All plots display significant wave attenuation/distortion, while final unrecoverable displacements are produced by soil plastifications (Figs. 16, 17a). Steady irreversible deformations are associated with prominent static components (at nil frequency) in the Fourier amplitude spectrum (Figs. 16, 17c), not present in the input Ormsby wavelet (Fig. 2b); -the numerical representation of wavelengths is dominated by soil plasticity, producing more deviation from the input waveform than variations in grid spacing. For this reason, only two x values have been used in this subsection for illustrative purposes, whereas EM/PM plots have been deemed not necessary; -the influence of x seems slightly magnified when lower h values, and thus lower elastic-plastic stiffness, are set (see Fig. 17). It is indeed not surprising that wave propagation in softer media may be more affected by space discretization, as in linear problems. However, it should be noted that x mainly influences the final irreversible displacement (Fig. 17b, c), which leads to presume substantial interplay of grid effects and constitutive time integration; -since the effects of x reduction are quite small in both time and frequency domains (for a given t), there is no strong motivation to suggest �x = V s /20f max . �x = V s /10f max = �x std should be actually appropriate in common practical situations, as long as no soil failure mechanisms are triggered -as for example in seismic  slope stability problems [17]. The occurrence of soil failure may introduce additional discretization requirements for an accurate representation of the collapse mechanism.
In addition, Fig. 18 illustrates the shear stress-strain VMKH response at the deepest integration (Gauss) point of the VMKH sub-layer. The material response is bilinear (elastic and elastic-plastic), with the elastic stiffness recovered upon stress reversal until new yielding occurs [40]. As mentioned above, the observable (small) differences in stress-strain response at different x may not be straightforwardly attributed to grid spacing deficiencies, but rather to the coupled influence of discretization in space and time on the global dynamics of the system.
The influence of the time-step size is illustrated for cases VMKH3-5 (Table 2) in Figs. 19, 20, encompassing three h values (0.5E, 0.05E and 0.01E) and also including EM/PM plots (Fig. 19d). In the lack of analytical solutions, misfits have been determined on the basis of a "sufficiently accurate" reference solution, here obtained numerically by setting �t = �t std /5 = 0.0001 s. For a relatively small input amplitude (A = 0.1 mm), convergence seems overall quite fast, and even t = t std results in both EM and PM values Both findings are likely related to the influence of time discretization on the residual displacement, which is larger than on other response variables. In fact, variations in the accumulated displacement mainly affect the envelope of the output signal, not its phase attributes. However, such a non-monotonic relationship between h and displacement EM/PM values has not been detected in the corresponding acceleration EM/PM plots (not reported here), due to the obvious lack of residual accelerations.

Influence of input motion amplitude
In non-linear problems, it is hard to draw general conclusions on the interaction between space/time discretization and input amplitude. The latter governs the amount of soil non-linearity mobilized and the resulting local stiffness, in turn affecting the requirements for accurate constitutive integration. In Fig. 21, the parametric study in Figs. 16, 17 is replicated for a higher input amplitude (A = 1 mm) and the same two different h values (cases VMKH6-7 in Table 2). The time-domain plots provided testify the effects of grid spacing on the predicted response: again, they mostly concern the final residual displacement, more pronouncedly as h decreases. The same previous uncertainties about the interplay of grid spacing and constitutive integration still apply to this case.
The discussion on the influence of t at higher input amplitude refers to Figs. 22, 23, illustrating the results obtained for x = x std and h equal to 0.5E, 0.05E and 0.01E (cases VMKH8-10 in Table 2); EM/PM plots comes from the numerical reference solution corresponding with �t = �t std /5 = 0.0001 s.
The comparison of Figs. 21 and 22 suggests that, even with a much larger input amplitude, x = x std is still an appropriate grid spacing for elastic-plastic problems, as long as t is substantially reduced to comply with (explicit) constitutive integration requirements. This inference is supported by the following observations: t affects not only the residual component of displacement time histories (as in Fig. 21), but also their maximum/minimum transient values -i.e., the numerical representation of plastic dissipation. This is clearly visible in Fig. 22a; -EM/PM values are in general higher at larger input amplitude (Fig. 22d), and experience a slower decrease as t is reduced (still depending on the specific h value); -the shear stress-strain loops in Fig. 23 show how inaccurate the simulated constitutive response can be when t is too large (e.g., t = 0.001 s) and substantial plastic degradation of material stiffness takes place (see the case h = 0.01E).
This set of results suggests that t should be at least in the order of �x/20V s for acceptable constitutive integration and overall accuracy in elastic-plastic simulations. However, this heuristic conclusion may be altered by the use of different material models (see next section) and stress-point algorithms.

Model parameters and parametric analysis
The influence of space/time discretization is now explored in combination with the non-linear PBS soil model introduced in Sect. 2.2.3 [48]. As in real geomaterials, the PBS model features an elastic-plastic response since the very onset of loading (vanishing yield locus), with the stiffness smoothly evolving from small-strain elastic behavior to failure (nil stiffness).
The results presented hereafter concern a 500 m thick soil layer, whose upper 100 m are made of a non-linear PBS soil resting on a 400 m elastic sub-layer. As done for the VMKH simulations, a thin layer (2.5 m) of elastic material has been added to prevent numerical problems with very strong motions and the whip effect at the ground surface. Input 1 with A = 1 mm has been exclusively considered, along with the following set of PBS parameters [48]  The list of PBS simulations is reported in Table 3, while the next figures will also illustrate the good performance of the PBS model in reproducing the cyclic soil behavior.   -grid spacing turns out to be influential again (Figs. 24,26), as a consequence of more severe variations (than in VMKH cases) in shear stiffness during cyclic loading. In fact, one would have to follow the stiffness reduction  curves arising from the constitutive response, and use minimum stiffness to decide on space discretization; -as in VMKH simulations, grid spacing mainly affects residual displacements. This is clearly shown by the EM/PM plots in Fig. 26b, where EM errors larger than 10% arise even when a very small time-step size is used (�t = �t std /25 = 0.00002 s); conversely, phase misfits are less affected by residual displacements and thus always quite limited. In presence of high nonlinearity, it seems safer to use x 4 ÷ 5 times smaller than �x std = V /10f max ; -the combination of explicit constitutive integration and high non-linearity makes time-stepping effects quite prominent, as is shown by Figs. 27 and 28. Further, Fig. 29 leads to conclude that �t = �t std /50 may be needed to obtain EM errors lower than 10 % (Figs. 29,30). Apparently, analysts have to compromise on accuracy and computational costs in these situations; -as expected, the shear stress-strain cycles in Figs. 25 and 28 show that the sensitivity to discretization builds up as increasing non-linearity is mobilized. This is the case for instance at the top of the PBS layer, where cycles are more dissipative than at the bottom due to lower overburden stresses and dynamic amplification.
Since displacement components result from strains through spatial integration, the displacement performance can be well-predicted on condition that strains are accurately computed all along the soil domain. For the same reason, the discretization requirements for displacement  Fig. 1) and at different x and t. These figures clearly point out that accuracy requirements may be more or less hard to satisfy depending on the specific spatial location. In 1D wave propagation problems, faster convergence is attained far from the ground surface, since it requires satisfactory accuracy in a lower number of nodes and integration points. Conversely, the close relationship between plastic strains and residual displacements has slender influence on acceleration components. In this respect, Figs. 33 and 34 show that, as long as reasonable grid spacing is set (possibly in the order of �x std /2 = V s /20f max ), the sensitivity of acceleration components to t is much weaker than for residual displacements.

Concluding remarks
Previously established criteria for space/time discretization in wave propagation FE simulations have been re-appraised and critically discussed to strengthen verification procedures in Computational Dynamics. The 1D propagation of seismic shear waves (Ormsby wavelets) through both linear and non-linear (elastic-plastic) media has been numerically simulated, with focus on capturing high-frequency motion and exploring the relationship between material response and discretization effects. After initial linear computations, two different non-linear material models (referred to as VMKH and PBS) have been used at increasing level of complexity. The main conclusions inferred are hereafter summarized: Elastic simulations Setting grid spacing (element size) and time-step size as per standard rules (�x std = V s /10f max and �t std = �t/V s ) has proven not always appropriate, especially to reproduce high-frequency motion components (this can be clearly visualized in the Fourier phase plane). When linear elements (8-node bricks) are used, �x ≈ �x std /2 and �t ≈ �t std /2 seem to ensure sufficient accuracy over the whole frequency range (both in amplitude and phase); higher-order elements (e.g., 27-node bricks) will allow the use of x = x std still in combination with �t ≈ �t std /2. Preserving accuracy in simulations with large domains and/or time durations seems intrinsically more difficult, since attenuation/dispersion phenomena are cumulative.
Elastic-plastic simulations Conclusive criteria for elastic-plastic problems can be hardly established, as space/ time discretization also interferes with the integration of non-linear constitutive equations. In this respect, different outcomes may be found depending on (i) kind of non-linearity associated with the material model (stiffness variations during straining), (ii) stress-point integration algorithm (e.g., explicit or implicit), (iii) input motion amplitude. The experience gained through the use of the PBS model (explicitly integrated in 8-node brick elements) suggests that �x std = V s /10f max and �t std = �x/10V s may need to be reduced by factors up to 4 ÷ 5 and 50, respectively, in the presence of strong input motions and severe stiffness variations. Importantly, these conclusions also depend on which output component is considered and where within the computational domain.
The present study is, however, not conclusive, especially when it comes to non-linear elastic-plastic problems. There are in fact several aspects that will deserve in the future further consideration, such as the implications of using higher-order finite elements. The same comment applies to geometrical effects (e.g., wave scattering) in 2D/3D problems, whose influence on discretization criteria for elastic-plastic simulations would be per se a whole research topic.