\documentclass[12pt]{article}
%\usepackage{subfigure} % was: subfignarms.sty
%\usepackage{psfig}
\usepackage{graphicx}
%\usepackage{timesmt} % use times typeface
%\usepackage{chicnarm} % < MOD
%\usepackage{chicago} % < MOD
\usepackage[numbers,round,colon,authoryear]{natbib}
\usepackage{array,arydshln}
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \usepackage{fancyheadings}
%
\usepackage{url}
\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{url}
\usepackage{multirow}
%\usepackage{fullpage}
\usepackage{lscape}
%\usepackage{subfig}
\usepackage{subfigure}
%\usepackage{slashbox}
%
% Tweak the graphics placement.
\renewcommand{\textfraction}{.1}
\renewcommand{\bottomfraction}{.6}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%12pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%12pt% this is for 12pt
\setlength{\textwidth}{16.4cm}
\setlength{\textheight}{22.8cm}
\setlength{\hoffset}{1.3cm}
\setlength{\voffset}{1.6cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%11pt%11pt
%11pt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%11pt % this is for 11pt
%11pt\setlength{\textwidth}{16.5cm}
%11pt\setlength{\textheight}{22.7cm}
%11pt\setlength{\hoffset}{2.0cm}
%11pt\setlength{\voffset}{2.3cm}
%11pt%11pt%11pt
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \lhead{}
% \rhead{}
% \chead{}
% \lfoot{}
% \rfoot{}
% %\cfoot{\pagenumber}
% \cfoot{}
% \addtolength{\headheight}{14pt}
%
% \setlength{\headrulewidth}{0mm}
% \setlength{\headrulewidth}{0mm}
% %\setlength{\footrulewidth}{0.1mm}
% %\setlength{\footrulewidth}{0.1mm}
%
% \pagestyle{fancy}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FEI LOGO definition
\newcommand{\FEI}
{{\sf \uppercase{F}\kern0.20em\uppercase{E}\kern0.20em\uppercase{I}}}
%
% % to change spacing in subcaption, see
% % http://www.ctan.org/texarchive/obsolete/macros/latex/contrib/subfigure/subfigure.pdf
% \addtolength{\subfigcapskip}{4mm}
% \addtolength{\subfigbottomskip}{4mm}
%
% change to sans serif font for all of the document
\renewcommand\familydefault{cmss}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\today
{\number\day.\space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}
% 
%
% TIME OF DAY
%
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% razmak izmedju redova
%
%baselinestretch
%
\renewcommand{\baselinestretch}{1.2}
\small\normalsize % trick from Kopka and Daly p47.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% All this stuff is from modifying the article.cls for Balkema
%%% specifications.
%\title{...}
%\author{...}
%\abstract{...}
\begin{document}
%\maketitle
\begin{center}
{\Large \bf Discretization effects in the finite element
simulation\\ of seismic shear waves in elastic and elasticplastic media}
\end{center}
%\vspace*{5mm}
%
\begin{center}
Kohei Watanabe${}^{1}$, Federico Pisan{\`o}${}^{2}$, Boris Jeremi{\'c}${}^{3,4}$
\end{center}
%\end{flushleft}
%\begin{flushleft}
\begin{flushleft}
\vspace*{5mm}
{\small ${}^{1}$~Shimizu Corporation, Tokyo, Japan} \\
{\small ${}^{2}$~Delft University of Technology, Delft, The Netherlands} \\
{\small ${}^{3}$~University of California, Davis, California, U.S.A.} \\
{\small ${}^{4}$~Lawrence Berkeley National Laboratory, Berkeley, California, U.S.A. } \\
{\small Corresponding author: B. Jeremi{\'c}, phone~(530) 7549248. fax~(530) 7547872,
Email:~\texttt{jeremic@ucdavis.edu}}
\end{flushleft}
\date{}% No date.
%\begin{flushleft}
%{\small \it Version: \today{}, \hhmm{} }
%\end{flushleft}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{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
unsatisfactory results are obtained. In particular, detailed analysis of
discretization criteria is developed for wave propagation through both elastic and
elasticplastic materials.
The establishment of such criteria is especially important when accurate prediction
of highfrequency motion is needed and/or in the presence of markedly nonlinear
material models.
Current discretization rules for wave problems in solids are critically
assessed as a \textit{conditio sine qua non} for improving
verification/validation procedures in applied seismology and earthquake
engineering. For this purpose, propagation of shear waves through a
%boundary condition controlled
1D stack of 3D finite elements have been
performed, including the use of wideband input motions in combination with both
linear elastic and nonlinear elasticplastic 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 showing the dependence of discretization requirements may be
requirements on time duration, spatial location, material model and the specific output variable
considered. \\
%
\textbf{keywords}: wave propagation, seismic, discretization, elastic, elasticplastic, verification
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{Introduction}
The study of wave motion is of utmost importance in many applied sciences, and
supports the understanding of transient phenomena in most natural and anthropic
dynamic systems. In this field, the case of mechanical waves in solid media
plays a major role because of its connection to a number of hazardous events,
such as blasts, earthquakes, structural vibrations, etc.
\citep{Semblat_Pecker_2009}. In this respect, seismic waves propagating through
the earth crust deserve the highest consideration, especially in the light of
their destructive potential and socioeconomical impact.
In the last decades, mathematicians, geophysicists and engineers devoted massive
research efforts to enhance the capability of predicting/mitigating the effects
of earthquakes on natural environments and human lives. Depending on the
analysis approach adopted, feasible solutions of seismic wave problems can be
classified as:
\begin{itemize}
\item [] \textit{Analytical}, i.e. exact mathematical solutions of
dynamic boundary value problems. In most cases, these account for
idealized/simplified geometries, boundary conditions, etc. \citep{kolsky1963,
graff1975, kausel2000, lai2005, kausel2006};
\item [] \textit{Numerical}, i.e. obtained through approximate
computational methods and also possible in the presence of inhomogeneous
properties, geometrical irregularities, material nonlinearities, etc.
\citep{argyrismlejnek1991, kramer1996, zienkiewicz_etal_1999}.
\end{itemize}
When linear (elastic) wave problems are considered, either timedomain or
frequencydomain solutions may be sought, whereas timedomain approaches are
usually needed in the presence of nonlinearities (constitutive or geometrical).
In this respect, it is worth remarking that most of the interest in
seismology, applied geophysics and earthquake
engineering is nowadays on nonlinear wave phenomena, since they govern (i) the
occurrence of instabilities (e.g. soil liquefaction and strain localization
\citep{ishihara1996,zienkiewicz_etal_1999,diprisco2012} and related
catastrophic events in earth crust materials (landslides, avalanches, debris
flows, soil sinking, etc.); (ii) the energetic interaction between geomaterials
and manmade structures
\citep{Wolf_1985,gazetas1998,Chopra_2000,jeremic2009,Semblat_Pecker_2009,randolph2011,dipriscopisano2011}.
It is thus apparent that reliable numerical simulations of seismic motion and
earthquakesoilstructure interaction can be only performed by means of
highfidelity computational tools, capable of coping with the remarkable
complexities in the aforementioned problems. The accuracy of numerical
predictions is affected  at least  by the following four factors:
\begin{enumerate}
\item selection of a suitable numerical algorithm for solving the governing
equations of motions;
\item constitutive modeling, i.e. the mathematical description of the material
behavior;
\item computer implementation of the numerical algorithm and constitutive
equations;
\item setup of the computational model, that is, a numerical, discrete
representation of physical reality.
\end{enumerate}
%
Assessing the above four points is the main core of a
thorough verification and validation process \citep{oberkampf2004,
babuska2004,roy2011}: 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, i.e. on the setup of
seismic wave propagation models via the Finite Element (FE)
approach\footnote{The extensive literature on the Finite Difference Method is
not mentioned here}.
%
The main issue is on the selection of proper time step and element size. In the context of applied
seismology and earthquake engineering, the problem seems to have been solved
quite long ago in the form of a few ``rules of thumb'' for space/time
discretization \citep{lysmer1969, kuhlemeyer1973}. Since then, not many works on
the subject have been published to the authors' knowledge
\citep{smith1975,bayliss1985,Bao98,basabe2007}, so that most new papers on seismic
wave problems tend to take the aforementioned rules for granted. Furthermore,
the relationship between discretization and accuracy in wave problems has been
mainly investigated by performing theoretical analyses for numerical attenuation
and dispersion, which is doable only for linear problems.
In the light the above premises, the authors aim an uptodate contribution to the matter, also accounting for
the increased importance assumed in recent years by nonlinear,
elasticplastic wave problems.
For the sake of clarity, the key points of the present work are hereafter
summarized:
\begin{itemize}
\item[] only 1D shear wave propagation tests are considered, in order to
benefit from the easier (partly analytically supported) interpretation of the
numerical outcomes;
\item[] discretization effects are illustrated in both the time and
frequency domains, \textbf{and then quantified via modern misfit criteria formulated in the full timefrequency domain};
\item[] since discretization effects depend in general on the numerical algorithm
adopted, a widespread finite element approximation scheme has been here adopted;
\item[] the role of constitutive nonlinearities is discussed;
\item[] the whole study has been conceived as a numerical ``falsification
test'' for the ``rules of thumb'' previously mentioned \citep{lysmer1969,
kuhlemeyer1973}.
\end{itemize}
%
The ultimate goal of this work is to reopen the
debate on the accuracy of wave simulations from a verification/validation
perspective, and in the presence of constitutive nonlinearities. The results
reported provide renovated critical insight into, and review of, traditional
discretization rules for practical simulation purposes.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\section{FE modeling of 1D seismic wave propagation}
\label{SWPM}
1D shear wave problems originate from those ideal situations where wave propagation is nearly vertical, with no lateral
geometrical/material inhomogeneities. In these conditions, all vertical crosssection can be regarded as symmetry planes and the soil deposit undergoes a ``double planestrain'' deformation, with both horizontal direct
strains prevented by symmetry \citep{prevost1989, borja1999}. As a consequence,
unknown variables only depend on time and the vertical spatial coordinate
(the problem is geometrically onedimensional), though the stress state is still
multiaxial \citep{diprisco_pisano2012}. The initialboundary value problem
under consideration is sketched in Figure \ref{fig:sketchIBVP}.
\begin{figure}[!htb]
\begin{center}
\includegraphics[width=9cm]{figures/IBVP.pdf}
\end{center}
\caption{One dimensional (1D) shear wave propagation through a soil layer}
\label{fig:sketchIBVP}
\end{figure}
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 all modeling and simulation.
The Real ESSI Simulator is a software, hardware and documentation system
developed specifically for high fidelity, realistic modeling and simulation of
earthquake soil structure interaction (ESSI).
%
The Real ESSI program
features a number of simple and advanced modeling features.
%
For example, on the finite element side, available are
solids elements (8, 20, 27, 827 node, dry and saturated bricks), structural elements
(trusses, beams, shells), contact elements (frictional slip and gap, dry and saturated), isolator and dissipator elements.
%
On the material modeling side, available are elastic (isotropic, anisotropic,
linear and nonlinear) and elasticplastic models
(isotropic, anisotropic hardening).
%
The seismic input can be applied using the Domain Reduction Method \citep{Bielak2001,Yoshimura2001}.
%
Both sequential and parallel versions (the latter uses the Plastic Domain Decomposition
(PDD) method \citep{Jeremic2008a} are available. \textbf{Recent applications of Real ESSI to seismic problems are documented, for instance, in \cite{jeremic2008, jeremic2009,cheng2009,taiebat2010,
jeremic2013smirt,jeremic2013,tasiopoulou1a2015,tasiopoulou2015,jeremic2015smirt}.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Space discretization and time marching}
\label{NS}
The essence of the finite element (FE) approximation scheme lies in the way in which the governing field
equations are discretized in space and time. A number of options are
available \citep{zienkiewicztaylor1}, however, the most widespread approach for Solid
Dynamics has been here followed in order to draw helpful conclusions for as many users as possible.
The Real ESSI program is based on a standard displacement FE formulation, where
displacement components are taken as unknown variables in the numerical
approximation \citep{zienkiewicztaylor1}. As for space discretization, the 1D FE
model has been built by using a stack of properly constrained 3D brick elements
 as was previously done, for instance, by \cite{borja1999}. Real ESSI program enables the
use of 8, 20, 827, and 27node elements, so that a number of options are given in terms of
spatial interpolation degree.
% the number of Gauss quadrature points is according
%standard rules for ensuring accurate spatial integration
%\citep{zienkiewicztaylor1, bathe2006}.
The wellknown Newmark method has been adopted for time marching
\citep{newmark1959}. The main feature
of the integration algorithm relates to the approximate series expansion for
displacement and velocity components, $u$ and $\dot{u}$ respectively:
%%%%
\begin{eqnarray}
{}^{n+1} u &=& {}^n u + \Delta t \; {}^n \dot{u} + \Delta^2 t \; \left[\left(\frac{1}{2}  \beta\right) \; {}^n \ddot{u} + \beta \; {}^{n+1} \ddot{u}\right]
\label{upU_Newmark_EQ1} \\
{}^{n+1} \dot{u} &=& {}^n \dot{u} + \Delta t \; \left[\left(1  \gamma\right) \; {}^n \ddot{u} + \gamma \; {}^{n+1} \ddot{u}\right]
\label{upU_Newmark_EQ2}
\end{eqnarray}
%%%%
between two subsequent timesteps $n$ and $n+1$. Importantly, the expansion
uses two parameters, $\beta$ and $\gamma$, governing the accuracy and stability
properties of the algorithm. While the reader is addressed to \cite{hughes2012}
for an exhaustive mathematical analysis, it is worth reminding that the
algorithm is unconditionally stable as long as:
%%%%
\begin{equation}
\gamma \ge \dfrac{1}{2}, \;\;\; \beta = \dfrac{1}{4}\left(\gamma+\dfrac{1}{2}\right)^2
\label{upU_Newmark_ab}
\end{equation}
%%%
$\gamma = 1/2$ is required for secondorder accuracy, whereas any $\gamma$ larger
than $1/2$ introduces numerical attenuation (damping). In this study, the
pair $\gamma=0.5$ and $\beta=0.25$ (no algorithmic/numerical
dissipation) have been exclusively considered.
%The family of Newmark
%time integration methods includes very wellknown members, such as:
%\begin{itemize}
% \item[] average acceleration method for $\beta = 1/4$ and $\gamma = 1/2$;
% \item[] linear acceleration method for $\beta=1/6$ and $\gamma = 1/2$
%\item[] (explicit) central difference method for $\beta = 0$ and $\gamma =
%1/2$.
%\end{itemize}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Material modeling}
\label{MM}
The Real ESSI program provides a number of material modeling options, ranging from
the simple linearelastic to advanced elasticplastic constitutive relationships
for cyclically loaded soils \citep{zienkiewicz_etal_1999, diprisco2012}.
Hereafter, the material models adopted for wave propagation analyses are briefly
described, namely (i) the standard linear elastic material model, (ii) the
elasticplastic von Mises model with linear kinematic hardening
\citep{lemaitre1990, Jeremic2004d} and
(iii) the bounding surface elasticplastic model by \cite{pisano2014}.
%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Linear elastic model}
Discretization issues will be first investigated for fully linear elastic problems. While the main concepts in
linear elasticity relevant to wave dynamics can be found in \cite{graff1975}, it
is worth reminding the relationship between the shear wave velocity $V_s$
and the two elastic constitutive parameters (Young's modulus $E$ and Poisson's ratio
$\nu$):
%
\begin{equation}
\label{eq_vs}
V_s=\sqrt{\dfrac{E}{2\left(1+\nu\right)\rho}}=\sqrt{\dfrac{G}{\rho}}
\end{equation}
%
where $\rho$ is the soil mass density and $G=E/(2 \left( 1+\nu \right))$ the
elastic shear modulus. As will be stressed in the following, $V_s$ plays a major
role in the selection of proper element size and timestep.
%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Elasticplastic: von Mises kinematic hardening (VMKH) model}
\label{VMKHmodel}
The relationship among discretization, accuracy and material nonlinearity will
be first explored through an elasticplastic von Mises constitutive law with
linear kinematic hardening, of the same kind described in \cite{lemaitre1990, Jeremic2004d}.
The VMKH model is very wellknown in literature and widely employed for
cyclically loaded metals; the application to soil dynamics is limited to
undrained loading conditions, i.e. when pore water drainage is prevented and
soils can be regarded as cohesive materials in the framework of total stress analysis
\citep{zienkiewicz_etal_1999, nova2012}. Although the assumption of linear
hardening is not the most accurate for soils\footnote{Nonlinear hardening models should
rather be used  see e.g. \cite{borjaamies1994,borja1999}}, it has been here introduced for
numerical convenience. In fact, owing to linear hardening, the postyielding
stiffness is constant, not straindependent: this implies an unrealistic
unbounded strength, but allows to identify the elasticplastic shear stiffness
with no ambiguity. Only four constitutive parameters need to be set:
\begin{itemize}
\item[] two elastic parameters  $E$ and $\nu$  for the preyielding
elastic response;
\item[] one yielding parameter  $k$  proportional to the initial size of
the cylindrical yield locus in the stress space;
\item[] one hardening parameter  $h$  governing the postyielding
(elasticplastic) stiffness.
\end{itemize}
The results from wave propagation analyses will clearly show the kind of
stressstrain response arising from the VMKH formulation.
%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Elasticplastic: Pisan{\`o} bounding surface (PBS) model }
\label{sec:PBS}
The more sophisticated constitutive relationship recently proposed by
\cite{pisano2014} has been also used. At variance with the aforementioned VMKH
formulation, the Pisan{\`o} bounding surface (PBS) model can quite accurately reproduce important aspects of monotonic/cyclic soil behavior \citep{kramer1996,Wood_2004,nova2012}, such as:
%
\begin{enumerate}
\item development of inelastic strains at the very onset of loading. This is
properly reproduced by exploiting the concept of ``vanishing yield
locus'';
\item frictional shear strength, i.e. depending on the effective confining
pressure;
\item nonlinear hardening, implying a continuous transition from smallstrain
to failure stiffness;
\item coupling between shear and volumetric strains;
\item average stiffness degradation and damping during cyclic shear loading.
\end{enumerate}
%
A remarkable quality of the PBS constitutive formulation is the low number of
input parameters required (only seven), which makes the model particularly
suitable for practical use:
\begin{itemize}
\item[] two elastic parameters  $E$ and $\nu$  to characterize the
material behavior at vanishing strains;
\item[] one shear strength parameter  $M$  directly related to the
material frictional angle;
\item[] two parameters  $k_d$ and $\xi$  governing the development of
plastic volumetric strains during shearing;
\item[] two hardening parameters  $h$ and $m$  to be identified on the
basis of stiffness degradation and damping cyclic curves.
\end{itemize}
For the sake of brevity, interested readers are addressed to \cite{pisano2014}
for details about formulation, performance and calibration. Evidence of the
simulated PBS stressstrain response under cyclic loading will be given in the
following sections.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Initial/boundary conditions and input motion}
\label{IM}
All the FE results hereafter presented have been obtained under the following
initial/boundary conditions (Figure \ref{fig:sketchIBVP}):
\begin{enumerate}
\item the system is initially still (nil initial velocities and accelerations)
with an atrest geostatic stress state \citep{nova2012};
\item a $x$displacement time history is imposed at the bottom boundary to
reproduce rigid bedrock conditions;
\item the top boundary is unloaded (free surface);
\item the aforementioned ``double planestrain'' conditions has been achieved by
(i) preventing the $y$displacement throughout the whole model, and (ii)
imposing master/slave connections to nodes at the same level $z$ (tied nodes).
\end{enumerate}
%
The last item in the list regards the fulfilment of 1D symmetry requirements and
guarantees the brick stack to undergo a ``shear beam''like deformation
\citep{kramer1996}.
As far as the input displacement is concerned, the Ormsby wavelet
\citep{ryan1994} fits the authors' intent:
\begin{equation}
\label{eq:ormsby}
\begin{split}
u\left(t\right)
=
A\left[\dfrac{\left(\pi f_4\right)^2}{\pi f_4\pi f_3}\text{sinc}^2\left(\pi f_4 t\right)

\dfrac{\left(\pi f_3\right)^2}{\pi f_4\pi f_3}\text{sinc}^2\left(\pi f_3 t\right)\right]

\\
\left[\dfrac{\left(\pi f_2\right)^2}{\pi f_2\pi f_1}\text{sinc}^2\left(\pi f_2 t\right)

\dfrac{\left(\pi f_1\right)^2}{\pi f_2\pi f_1}\text{sinc}^2\left(\pi f_1 t\right)\right]
\end{split}
\end{equation}
%
where $t$ denotes the physical time and $A$ the signal amplitude,
$\text{sinc}\left(x\right) =\left(\sin x\right) /x$ is the cardinal sine
function, $f_i$ ($i=1,2,3,4$) stand for the lowcut, lowpass, highcut and
highpass frequencies, respectively. The meaning of the $f_i$ frequencies can be inferred from Figure \ref{fig:ormsby_b}, illustrating the amplitude Fourier
spectrum of function \eqref{eq:ormsby}. In particular, the suitability of the
Ormsby wavelet has a twofold motivation:
\begin{enumerate}
\item function \eqref{eq:ormsby} has a number of sign reversals and will induce
several loading/unloading cycles into the soil undergoing wave motion (Figure
\ref{fig:ormsby_a});
\item the peculiar flat branch in the amplitude Fourier spectrum (Figure
\ref{fig:ormsby_b}) is a very convenient feature for accuracy assessments in the
frequency domain (see next section).
\end{enumerate}
\begin{figure}[!htb]
\centering
\subfigure[Time history]
{
\includegraphics[width=0.7\textwidth]{figures/Ormsby_20Hz_1mm_timehistory.pdf}
\label{fig:ormsby_a}
}
\
\subfigure[Amplitude of Fourier spectrum]
{
\includegraphics[width=0.7\textwidth]{figures/Ormsby_20Hz_1mm_fft.pdf}
\label{fig:ormsby_b}
}
\caption{Ormsby wavelet ($f_1$=0.1 Hz, $f_2$=1 Hz, $f_3$=18 Hz, $f_4$=20 Hz)}
\label{fig:ormsby}
\end{figure}
The above features of the Ormsby wavelet will allow to explore discretization
effects over frequency ranges of choice. In fact, 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.
\subsection{Misfit criteria}
\label{sec:misfit}
\textbf{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.}
\textbf{A significant improvement in this area was introduced by \cite{kristekova2006}, who compared seismograms on the basis of the timefrequency representation (TFR) obtained through continuous wavelet
transformation \citep{holschneider1995}. The TFR of signal misfit allows to extract the time
evolution of the spectral content, and thus to define the following local timefrequency envelope difference:}
\begin{equation}
\label{eq:E}
\Delta E\left(t,f\right)=\vert W\left(t,f\right)\vert\vert W_{REF}\left(t,f\right)\vert
\end{equation}
\textbf{and timefrequency phase difference:}
\begin{equation}
\label{eq:P}
\Delta P\left(t,f\right)=
\vert W_{REF}\left(t,f\right)\vert
\dfrac{\text{arg}\left[W\left(t,f\right)\right]\text{arg}\left[W_{REF}\left(t,f\right)\right]}{\pi}
\end{equation}
\textbf{where $W\left(t,f\right)$ and $W_{REF}\left(t,f\right)$ are the TFR (wavelet transform) of the signal ``under evaluation'' and the reference seismogram, respectively.
As explained by \cite{kristekova2006}, it is also possible to obtain purely time or frequencydependent misfit measures by projecting $\Delta E$ and $\Delta P$ onto one of the two domains. In particular, the following singlevalues measures for envelope misfit (EM) }
\begin{equation}
\label{eq:EM}
EM=\sqrt{\dfrac{\sum_{f}\sum_t\vert\Delta E\left(t,f\right)\vert^2}{\sum_f\sum_t\vert W_{REF}\left(t,f\right)\vert^2}}
\end{equation}
\textbf{and phase misfit (PM)}
\begin{equation}
\label{eq:FM}
PM=\sqrt{\dfrac{\sum_f\sum_t\vert\Delta P\left(t,f\right)\vert^2}{\sum_f\sum_t\vert W_{REF}\left(t,f\right)\vert^2}}
\end{equation}
\textbf{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 carrying more information than a single amplitude value at given time.}
%
\textbf{While the theoretical background for the above misfit criteria is widely described by \cite{kristekova2006,kristekova2009}, opensource routines for misfit analysis are available at \url{www.nuquake.eu/Computer
Codes/} (TFMISFITS 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 \citep{perez2007,moczo2007,benjemaa2007,kaser2008,fichtner2008}.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\section{Linear elastic wave simulations}
\label{SWPS}
In this section, the influence of discretization on accuracy is first discussed
for fully linear elastic problems. In this context, the insight coming from the
wellknown analytical solution can be exploited to critically reconsider
standard rules for space/time discretization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Standard rules for space/time discretization}
\label{sec:rules}
The selection of appropriate grid spacing\footnote{Henceforth, $\Delta x$ will
always denote the vertical node spacing, coinciding with the element thickness
in the case of 8node bricks.} and timestep size is usually based on very simple rules. As for space discretization, \cite{lysmer1969} stated that ``\textit{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 \citep{jeremic2009}:
\begin{equation}
\label{eq:gridspace}
\Delta x\leq\dfrac{\lambda_{min}}{10}=\dfrac{V_s}{10f_{max}}
\end{equation}
On the other side, the timestep size also needs to be limited
to ensure accuracy and stability \citep{argyrismlejnek1991}. In principle, the smallest fundamental period of the system should be represented with about ten timesteps
 same as for space discretization. However, $\Delta 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 $\Delta t$ values):
\begin{equation}
\label{eq:step}
\Delta t\leq\dfrac{\Delta x}{V_s}
\end{equation}
%
Condition \eqref{eq:step} ensures algorithmic stability in many explicit schemes
for hyperbolic differential problems \citep{quarteroni2008}, though it is often
regarded as an accuracy criterion for implicit (unconditionally stable) time
marching as well (see section \ref{NS}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{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 $\rho=2000$ kg/m\textsuperscript{3}, $V_s=1000$ m/s
and $\nu=0.3$ (corresponding to $G=2$ GPa). No Rayleigh damping has been introduced.
As for the input motion, two different Ormsby wavelets have been employed,
corresponding with the following input parameters in Equation \eqref{eq:ormsby}:
\begin{itemize}
\item[] input 1: $f_1$=0.1 Hz, $f_2$=1 Hz, $f_3$=18 Hz, $f_4$=20 Hz (plotted
in Figure \ref{fig:ormsby});
\item[] input 2: $f_1$=0.1 Hz, $f_2$=1 Hz, $f_3$=45 Hz, $f_4$=50 Hz;
\item[] the amplitude parameter $A$ has been always set to produce at the
bottom layer a maximum displacement of 1 mm.
\end{itemize}
As previously mentioned (section \ref{IM}), both inputs 1 and 2 have been used
to explore the interplay between discretization effects and the width of the frequency range.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Discussion of numerical results}
The influence of grid spacing and timestep 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 consideration is also paid to accelerations,
postcalculated through secondorder central differentiation.
Table \ref{tab:elastic_analyses} provides a list of the comparative simulations
performed for fully linear problems. Each case is denoted by: (i) maximum frequency $f_{max}$ in the input wavelet ($f_4$ in \eqref{eq:ormsby}); (ii)
standard grid spacing $\Delta x_{std}$ and (iii) timestep size $\Delta t_{std}$ from discretization rules \eqref{eq:gridspace}\eqref{eq:step}; (iv)
$\Delta x$ and (v) $\Delta t$ actually used; (vi) type of brick elements adopted.
\begin{table}[h!]
\centering
\begin{tabular}{ccccccc}
\hline\hline
case \# & $f_{max}$ [Hz] & $\Delta x_{std}$ [m] & $\Delta t_{std}$ [s] & $\Delta x$ [m] & $\Delta t$ [s] & brick type\\
\hline
EL1 & 20& 5 & 0.005 & 2, 5, 10 & 0.005 & 8node\\
EL2 & 20& 5 & 0.005 & 2, 5, 10 & 0.002 & 8node\\
EL3 & 50& 2& 0.002 & 0.8, 2, 4 & 0.002 & 8node\\
EL4 & 50& 2& 0.002 & 0.8, 2, 4 & 0.001 & 8node\\
EL5 & 20& 5 & 0.005 & 2, 5, 10 & 0.002 & 27node\\
\hdashline
EL6 & 20 & 5 & 0.005 & 5 & 0.002, 0.005, 0.01 & 8node \\
EL7 & 20 & 5 & 0.005 & 2 & 0.001, 0.002, 0.005 & 8node \\
EL8 & 50 & 2 & 0.002 & 2 & 0.001, 0.002, 0.005 & 8node\\
EL9 & 50 & 2 & 0.002 & 0.8 & 0.0005, 0.001, 0.002 & 8node\\
EL10 & 20 & 5 & 0.005 & 5 & 0.002, 0.005, 0.01 & 27node\\
\hline\hline
\end{tabular}
\caption{List of elastic simulations}
\label{tab:elastic_analyses}
\end{table}
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).\textbf{ Additional quantitative insight is also gained through the EM and PM misfit criteria introduced in section \ref{sec:misfit}.} Unless differently stated,
numerical outputs at the top of the soil layer are considered.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Influence of grid spacing}
\textbf{Grid spacing effects at the top of the FE model are illustrated in Figures
\ref{fig:EL1disp}\ref{fig:EL5disp} for the cases EL1EL5 (Table
\ref{tab:elastic_analyses}) in terms of:
(ab) displacement time histories; (c) Fourier amplitudes and (d) phase differences at the surface; (e) EM and PM misfits (for each numerical solution, misfits are calculated with respect to the exact analytical solution). Starting from Figure \ref{fig:EL2disp}, displacement time histories are not compared with the input motion (as done in Figure \ref{fig:elastic_8node_20Hz_case4_dis1}) for the sake of brevity, whereas only a reduced time window around the output motion is displayed for a clearer visualisation (e.g. as in Figure \ref{fig:elastic_8node_20Hz_case4_dis2})}
%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (0.04.0 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_3a.pdf}%
\label{fig:elastic_8node_20Hz_case4_dis1}%
\vspace{10\baselineskip}%
}%
\
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_3b.pdf}%
\label{fig:elastic_8node_20Hz_case4_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_3c.pdf}%
\label{fig:elastic_8node_20Hz_case4_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_3d.pdf}%
\label{fig:elastic_8node_20Hz_case4_fphs}%
}%
\
\subfigure[EM/PM misfits (ref. solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig3.pdf}%
\label{fig:elastic_8node_20Hz_case4_mis}%
}%
\caption{Influence of grid spacing, displacement plot, case EL1 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=$ 2, 5, 10 m, $\Delta t=0.005$ s, 8node brick)}
\label{fig:EL1disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_4b.pdf}%
\label{fig:elastic_8node_20Hz_case3_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_4c.pdf}%
\label{fig:elastic_8node_20Hz_case3_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_4d.pdf}%
\label{fig:elastic_8node_20Hz_case3_fphs}%
}%
\
\subfigure[EM/PM misfits (ref. solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig4.pdf}%
\label{fig:elastic_8node_20Hz_case3_mis}%
}%
\caption{Influence of grid spacing, displacement plot, case EL2 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=$ 2, 5, 10 m, $\Delta t=0.002$ s, 8node brick)}
\label{fig:EL2disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.22.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_5b_r.pdf}%
\label{fig:elastic_8node_50Hz_case3_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_5c_r.pdf}%
\label{fig:elastic_8node_50Hz_case3_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_5d_r.pdf}%
\label{fig:elastic_8node_50Hz_case3_fphs}%
}%
\
\subfigure[EM/PM misfits (ref. solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig5r.pdf}%
\label{fig:elastic_8node_50Hz_case3_mis}%
}%
\caption{Influence of grid spacing, displacement plot, case EL4 ($f_{max}=50$ Hz, $\Delta x_{std}=2$ m, $\Delta t_{std}=0.002$ s, $\Delta x=$ 0.8, 2, 4 m, $\Delta t=0.001$ s, 8node brick)}
\label{fig:EL4disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_6b.pdf}%
\label{fig:elastic_27node_20Hz_case3_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_6c.pdf}%
\label{fig:elastic_27node_20Hz_case3_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_6d.pdf}%
\label{fig:elastic_27node_20Hz_case3_fphs}%
}%
\
\subfigure[EM/PM misfits (ref. solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig6.pdf}%
\label{fig:elastic_27node_20Hz_case3_mis}%
}%
\caption{Influence of grid spacing, displacement plot, case EL5 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=$ 2, 5, 10 m, $\Delta t=0.002$ s, 27node brick)}
\label{fig:EL5disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\begin{figure}[!htb]
\centering
\subfigure[Acceleration time history (2.2 3.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_7a.pdf}%
\label{fig:elastic_8node_acc_case1_acc4}%
}%
\subfigure[EM/PM misfits (ref. solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig7a.pdf}%
\label{fig:elastic_8node_acc_case1_acc4_mis}%
}%
\
\subfigure[Acceleration time history (2.2 3.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_7b.pdf}%
\label{fig:elastic_8node_acc_case1_acc3}%
}%
\
\subfigure[EM/PM misfits (ref. solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig7b.pdf}%
\label{fig:elastic_8node_acc_case1_acc3_mis}%
}%
\caption{Influence of grid spacing, acceleration plot, cases (ab) EL1 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=$ 2, 5, 10 m, $\Delta t=0.005$ s, 8node brick) and (cd) EL2 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=$ 2, 5, 10 m, $\Delta t=0.002$ s, 8node brick)}
\label{fig:EL1EL2acc}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Acceleration time history (2.2 2.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_8a_r.pdf}%
\label{fig:elastic_8node_acc_case2_acc4}%
} %
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig8ar.pdf}%
\label{fig:elastic_8node_acc_case2_acc4_mis}%
}%
\
\subfigure[Acceleration time history (2.2 2.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_8b_r.pdf}%
\label{fig:elastic_8node_acc_case2_acc3}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig8br.pdf}%
\label{fig:elastic_8node_acc_case2_mis}%
}%
\caption{Influence of grid spacing, acceleration plot, cases (ab) EL3 ($f_{max}=50$ Hz, $\Delta x_{std}=2$ m, $\Delta t_{std}=0.002$ s, $\Delta x=$ 0.8, 2, 4 m, $\Delta t=0.002$ s, 8node brick) and (cd) EL4 ($f_{max}=50$ Hz, $\Delta x_{std}=2$ m, $\Delta t_{std}=0.002$ s, $\Delta x=$ 0.8, 2, 4 m, $\Delta t=0.001$ s, 8node brick)}
\label{fig:EL3EL4acc}
\end{figure}
Figures \ref{fig:EL1disp}\ref{fig:EL5disp} suggest the following observations (some of which expected):
\begin{itemize}
\item[] \textbf{even though $\Delta x_{std}$ is set on the basis of the maximum frequency $f_{max}$, its suitability as grid spacing 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 Figures \ref{fig:EL1disp}(c) and \ref{fig:EL2disp}\ref{fig:EL5disp}(b)). Grid spacing affects output Fourier spectra both in amplitude and phase;}
\item[] \textbf{in all cases, envelope and phase misfits, EM and PM, are quantitatively very similar (Figures \ref{fig:EL1disp}(e) and \ref{fig:EL2disp}\ref{fig:EL5disp}(d));}
\item[] \textbf{reducing $\Delta x$ below $\Delta x_{std}$ is beneficial only if $\Delta t$ is also lower than $\Delta t_{std}$. This is apparent in Figure \ref{fig:EL1disp}(e), where an increase in EM and PM is observed as $\Delta x$ gets lower than $\Delta x_{std}$. Conversely, monotonic EM/PM trends are shown in Figures \ref{fig:EL2disp}\ref{fig:EL4disp}(d);}
\item[] \textbf{at given grid spacing $\Delta x$, reducing the timestep improves the numerical solution mostly in terms of Fourier phase, not amplitude (compares Figures \ref{fig:EL1disp}(cd) and \ref{fig:EL2disp}(bc)). It may be generally stated that, when $\Delta x$ is not appropriate, reducing the timestep size does not produce substantial improvements;}
\item[] \textbf{based on these initial examples, a grid spacing $\Delta x$ in the order of $V_s/20f_{max}=\Delta x_{std}/2$ ensures high accuracy (EM and PM $<$10\% ) in combination with $\Delta t=\Delta
x/2V_s=\Delta t_{std}/2$. These enhanced discretization rules hold for loworder FEs (8node brick elements) but are not affected by the frequency bandwidth of the input signal. In the latter respect, Figures \ref{fig:EL2disp}\ref{fig:EL4disp}(d) show quantitatively similar EMPM trends for $f_{max}$ equal to 20 Hz and 50 Hz. Also, minimum misfits are attained in the EL2 case, where a smaller $\Delta t/\Delta t_{std}$ ratio has been purposely set.}
\end{itemize}
\textbf{The above conclusions apply to 8node brick elements, while Figure \ref{fig:EL5disp} shows that ``ten nodes per
wavelength'' are still suitable when higherorder elements (here 27node bricks\footnote{\textbf{For a given number of nodes per wavelength, the size $\Delta x$ of 27node elements along the propagation direction is double than for 8node bricks.}}) are employed. However, this lighter
requirement for grid spacing seems to perform well in combination with $\Delta t\leq\Delta x/2V_s$, and results in EM and PM lower than 10\% even for $\Delta x/\Delta x_{std}=2$ (5 nodes per wavelength).}
\textbf{It is also important to evaluate grid spacing effects on acceleration components, as they will affect the
inertial forces transmitted to manmade 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 Figures \ref{fig:EL1EL2acc}
and \ref{fig:EL3EL4acc}, grid spacing plays qualitatively as in Figures \ref{fig:EL1disp}\ref{fig:EL4disp}, though the EM/PM trends  similar in shape  are shifted upwards. This means that, in the presence of loworder elements, more severe discretization requirements should be fulfilled if very accurate accelerations need to be computed.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
%\newpage
\subsubsection{Influence of timestep size}
For given grid spacings, the influence of $\Delta t$ has been studied by varying the timestep size with respect to the limit emerging from Equation \eqref{eq:step}, i.e. $\Delta t_{std}={\Delta x}/{V_s}$. Time discretization effects are illustrated in Figures \ref{fig:EL6disp}\ref{fig:EL8EL9acc} and suggest the following comments:
\begin{itemize}
\item[] as observed in the previous subsection, $\Delta t$ mainly affects the
Fourier phase, \textbf{with comparable EM and PM values in all cases}. Phase differences with respect to the
exact solution decrease as $\Delta t$ is reduced  see for instance in Figures \ref{fig:EL6disp}\ref{fig:EL10disp}(c);
\item[] \textbf{in combination with $\Delta x=V_s/20f_{max}=\Delta
x_{std}/2$, $\Delta t=\Delta
t_{std}$ may still result in some highfrequency phasedifference with the respect to the analytical solution, (Figures \ref{fig:EL6disp}\ref{fig:EL10disp}(c)). As found by investigating grid spacing effects, $\Delta t={\Delta x}/{2V_s}=\Delta t_{std}/2$ yields sufficient accuracy (EMPM lower than 10\%) to most practical purposes (see Figures \ref{fig:EL6disp}\ref{fig:EL10disp}(d))};
\item[] \textbf{when 27node bricks are used, the use of $\Delta x=\Delta x_{std}$ and $\Delta t\leq\Delta t_{std}/2$ is still an appropriate option, giving rise to EM and PM lower than 5\% (Figures \ref{fig:EL10disp}). Even in this case, discretization errors are still governed by phase differences, while excellent performance in terms of Fourier amplitude is observed};
\item[] \textbf{ Figures \ref{fig:EL6EL7acc} and \ref{fig:EL8EL9acc} show that the above inferences apply qualitative to acceleration time histories as well. However, EM and PM value are quite high (significantly larger than 10\%) while $\Delta t\geq\Delta t_{std}$, regardless of the grid spacing ratio. Accuracy is quickly gained when $\Delta t$ is reduced and $\Delta x<\Delta x_{std}/2$.}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_9b.pdf}%
\label{fig:elastic_8node_20Hz_case2_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_9c.pdf}%
\label{fig:elastic_8node_20Hz_case2_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_9d.pdf}%
\label{fig:elastic_8node_20Hz_case2_fphs}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig9.pdf}%
\label{fig:elastic_8node_20Hz_case2_mis}%
}%
\caption{Influence of timestep size, displacement plot, case EL6 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=
5$ m, $\Delta t=$ 0.002, 0.005, 0.010 s, 8node brick)}
\label{fig:EL6disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_10b.pdf}%
\label{fig:elastic_8node_20Hz_case1_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_10c.pdf}%
\label{fig:elastic_8node_20Hz_case1_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_10d.pdf}%
\label{fig:elastic_8node_20Hz_case1_fphs}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig10.pdf}%
\label{fig:elastic_8node_20Hz_case1_mis}%
}%
\caption{Influence of grid spacing, displacement plot, case EL7 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=2$ m, $\Delta t=$ 0.001, 0.002, 0.005 s, 8node brick)}
\label{fig:EL7disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.22.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_11b_r.pdf}%
\label{fig:elastic_8node_50Hz_case1_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_11c_r.pdf}%
\label{fig:elastic_8node_50Hz_case1_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_11d_r.pdf}%
\label{fig:elastic_8node_50Hz_case1_fphs}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig11r.pdf}%
\label{fig:elastic_8node_50Hz_case1_mis}%
}%
\caption{Influence of timestep size, displacement plot, case EL9 ($f_{max}=50$ Hz, $\Delta x_{std}=2$ m, $\Delta t_{std}=0.002$ s, $\Delta x=0.8$ m, $\Delta t=$ 0.0005, 0.001, 0.002 s, 8node brick)}
\label{fig:EL9disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_12b.pdf}%
\label{fig:elastic_27node_20Hz_case2_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_12c.pdf}%
\label{fig:elastic_27node_20Hz_case2_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_12d.pdf}%
\label{fig:elastic_27node_20Hz_case2_fphs}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig12.pdf}%
\label{fig:elastic_27node_20Hz_case2_mis}%
}%
\caption{Influence of timestep size, displacement plot, case EL10 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=5$ m, $\Delta t=$ 0.002, 0.005, 0.010 s, 27node brick)}
\label{fig:EL10disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\
\subfigure[Acceleration time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_13a.pdf}%
\label{fig:elastic_8node_acc_case1_acc2}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig13a.pdf}%
\label{fig:elastic_8node_acc_case1_acc2_mis}%
}%
\label{fig:elastic_8node_acc_case1}
\subfigure[Acceleration time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_13b.pdf}%
\label{fig:elastic_8node_acc_case1_acc1}%
}%
\
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig13b.pdf}%
\label{fig:elastic_8node_acc_case1_acc1_mis}%
}%
\caption{Influence of timestep size, acceleration plot, cases EL6 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=
5$ m, $\Delta t=$ 0.002, 0.005, 0.010 s, 8node brick) and EL7 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=2$ m, $\Delta t=$ 0.001, 0.002, 0.005 s, 8node brick)}
\label{fig:EL6EL7acc}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Acceleration time history (2.22.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_14a_r.pdf}%
\label{fig:elastic_8node_acc_case2_acc2}%
}%
\label{fig:elastic_8node_acc_case2}
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig14a.pdf}%
\label{fig:elastic_8node_acc_case2_acc2_mis}%
}%
\
\subfigure[Acceleration time history (2.22.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_14b_r.pdf}%
\label{fig:elastic_8node_acc_case2_acc1}%
}%
\subfigure[EM/PM misfits (ref: solution: analytical)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig14br.pdf}%
\label{fig:elastic_8node_acc_case2_acc1_mis}%
}%
\caption{Influence of timestep size, acceleration plot, cases EL8 ($f_{max}=50$ Hz, $\Delta x_{std}=2$ m, $\Delta t_{std}=0.002$ s, $\Delta x=
2$ m, $\Delta t=$ 0.001, 0.002, 0.005 s, 8node brick) and EL9 ($f_{max}=50$ Hz, $\Delta x_{std}=2$ m, $\Delta t_{std}=0.002$ s, $\Delta x=0.8$ m, $\Delta t=$ 0.0005, 0.001, 0.002 s, 8node brick)}
\label{fig:EL8EL9acc}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}
\begin{figure}[h!]
\centering
\subfigure[Displacement time history (0.013.0 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_15a.pdf}%
\label{fig:elastic_8node_20Hz_dispersion_dis1}%
}%
\
\subfigure[Displacement history (different time windows)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_15b.pdf}%
\label{fig:elastic_8node_20Hz_dispersion_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_15c.pdf}%
\label{fig:elastic_8node_20Hz_dispersion_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_15d.pdf}%
\label{fig:elastic_8node_20Hz_dispersion_fphs}%
}%
\caption{Time evolution of wave dispersion, displacement plot, case EL7 ($f_{max}=20$ Hz, $\Delta x_{std}=5$ m, $\Delta t_{std}=0.005$ s, $\Delta x=2$ m, $\Delta t=$ 0.001 s, 8node brick)}
\label{fig:EL7disp_dispersion}
\end{figure}
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 neverending wave motion is established. It is thus interesting to check how discretization errors propagate in time at the free surface, as is shown in Figure
\ref{fig:EL7disp_dispersion}. Subsequent wave arrivals are compared in the time (Figure
\ref{fig:EL7disp_dispersion}(a)(b)) and frequency (Figure
\ref{fig:EL7disp_dispersion}(c)(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 highfrequency phase difference is detected in Figure
\ref{fig:EL7disp_dispersion}(d), with negligible variation in Fourier amplitude (Figure
\ref{fig:EL7disp_dispersion}(c)). The effect of cumulative wave dispersion is apparent in Figure \ref{fig:EL7disp_dispersion}(a), and implies that selecting suitable $\Delta x$ and $\Delta t$ becomes increasingly delicate for large FE models and long durations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
%\clearpage
\section{Nonlinear elasticplastic wave simulations}
This section concerns discretization effects in presence of material nonlinearity. As most commonly done in Geomechanics
\citep{zienkiewicz_etal_1999}, the nonlinear cyclic response of geomaterials can be described in the framework of elastoplasticity, and here the VMKH and PM models described in Section \ref{MM} have been adopted.
Prior to presenting numerical results, some preliminary remarks should be made:
\begin{itemize}
\item[] the nonlinear 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 $\Delta x$$\Delta t$ refinement;
\item[] 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 hereafter on several expected pitfalls, all relevant to the global verification process \citep{oberkampf2004, babuska2004, roy2011};
\item[] the accuracy of nonlinear computations is highly affected by the input amplitude. This governs the amount of nonlinearity mobilized by wave motion and, as a consequence, the accuracy of numerical solutions at varying discretization.
\end{itemize}
In nonlinear (elasticplastic) problems, discretization is not only responsible for the numerical representation of waves
(dissipation, dispersion, stability), but also governs the accuracy of
constitutive integration \citep{simo1998,borja2013}. For instance, changes in timestep 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 elasticplastic response can be established for any stressstrain combination, (lowest) elasticplastic (shear)
stiffness may be used to develop suitable discretization via Equation~\ref{eq_vs}. Apparently, this approach assumes that the stressstrain 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 stressstrain 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 \citep{Desai84, Chen}. Although implicit algorithms may possess better accuracy/stability properties, explicit integration is often preferred for advanced constitutive formulations and cyclic loading \citep{jeremic2008}.
There is also wide consensus on the poor performance in elasticplastic computations of timestep sizes derived through elastic parameters and Equation \eqref{eq:step}, especially in combination with explicit stresspoint algorithms. For this reason, the following time marching rule may be regarded as an upper bound for nonlinear problems (instead of \eqref{eq:step}):
\begin{equation}
\label{eq:stepNL}
\Delta t\leq\dfrac{\Delta x}{10V_s}
\end{equation}
In the following, rules \eqref{eq:gridspace} and \eqref{eq:stepNL} will be assumed as starting discretization criteria and critically assessed. For shorter discussion, only input 1 ($f_{max}=$ 20 Hz) and 8node brick elements are employed for nonlinear simulations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{VMKH model}
%
\subsubsection{Model parameters and parametric analysis}
A heterogeneous 1 km thick soil deposit has been considered, formed by a 200 m thick VMKH sublayer 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 Section \ref{VMKHmodel}) have been set (same elastic parameters for both the VMKH and the elastic sublayers), with no algorithmic nor Rayleigh damping introduced in numerical computations.:
\begin{itemize}
\item[] mass density and elastic properties: $\rho=$ 2000 kg/m\textsuperscript{3},
$E=$ 5.2 GPa and $\nu=$ 0.3, whence the elastic shear wave velocity $V_s=$ 1000 m/s results (same elastic parameters employed for both the elastic and the VMKH sublayers);
\item[] yielding parameter (radius of the von Mises cylinder): $k=$ 10.4 kPa;
\item[] different $h$ values (hardening parameter) have been set: $h=$ 0.5$E$, 0.05$E$, 0.01$E$.
\end{itemize}
In the analysis of VMKH cases, the influence of the hardening parameter ($h$) and the input amplitude ($A$) is also considered, as they affect the material elasticplastic stiffness and the amount of plasticity mobilized. The VMKH simulation programme is reported in Table \ref{tab:VMKH_analyses}, where $\Delta t_{std}$ has been determined through Equation~\eqref{eq:stepNL} (i.e. $\Delta t_{std} = {\Delta x}/{10V_s}$).
\begin{table}[h!]
\centering
\begin{tabular}{ccccccc}
\hline
case \# & $\Delta x_{std}$ [m] & $\Delta t_{std}$ [s] & $\Delta x$ [m] & $\Delta t$ [s] & $h$ & $A$ [mm]\\
\hline
VMKH1 & 5 & 0.0005 & 1, 5 & 0.0001 & 0.5$E$ & 0.1 \\
VMKH2 & 5 & 0.0005 & 1, 5 & 0.0001 & 0.05$E$ & 0.1 \\
\hdashline
VMKH3 & 5 & 0.0005 &5 & 0.0002, 0.0005, 0.001 & 0.5$E$ & 0.1\\
VMKH4 & 5 & 0.0005&5 & 0.0002, 0.0005, 0.001 & 0.05$E$ & 0.1\\
VMKH5 & 5 & 0.0005 &5 & 0.0002, 0.0005, 0.001 & 0.01$E$ & 0.1\\
\hdashline\hdashline
VMKH6 & 5 & 0.0005 & 1, 5& 0.0001 & 0.5$E$ & 1\\
VMKH7 & 5 & 0.0005 &1, 5 & 0.0001 & 0.05$E$ & 1\\
\hdashline
VMKH8 & 5 & 0.0005 & 5 & 0.0002, 0.0005, 0.001 & 0.5$E$ & 1\\
VMKH9 & 5 & 0.0005 & 5 & 0.0002, 0.0005, 0.001 & 0.05$E$ & 1\\
VMKH10 & 5 & 0.0005 & 5 & 0.0002, 0.0005, 0.001 & 0.01$E$ & 1\\
\hline\hline
\end{tabular}
\caption{List of VMKH simulations}
\label{tab:VMKH_analyses}
\end{table}
%\clearpage
%\newpage
\subsubsection{Influence of grid spacing and timestep size}
%
The results in Figures~\ref{fig:VMKH1disp} and \ref{fig:VMKH2disp} exemplify the role played by space discretization in elasticplastic simulations. These results have been obtained by employing a timestep smaller than $\Delta t_{std}$ (cases VMKH12 in Table \ref{tab:VMKH_analyses}), 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:
%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (0.04.0 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_16a.pdf}%
\label{fig:vmlkh8node_dx_case3_dis1}%
}%
\
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_16b.pdf}%
\label{fig:vmlkh8node_dx_case3_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_16c.pdf}%
\label{fig:vmlkh8node_dx_case3_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_16d.pdf}%
\label{vmlkh8node_dx_case3_fphs.png}%
}%
\caption{Influence of grid spacing, displacement plot, case VMKH1 ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=$ 1, 5 m, $\Delta t=0.0001$ s, $h=0.5E$, $A=0.1$ mm)}
\label{fig:VMKH1disp}
\end{figure}
%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (0.04.0 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_17a.pdf}%
\label{fig:vmlkh8node_dx_case4_dis1}%
}%
\
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_17b.pdf}%
\label{fig:vmlkh8node_dx_case4_dis2}%
}%
\
\subfigure[Amplitude of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_17c.pdf}%
\label{fig:vmlkh8node_dx_case4_famp}%
}%
\
\subfigure[Phase difference of displacement Fourier spectrum]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_17d.pdf}%
\label{vmlkh8node_dx_case4_fphs.png}%
}%
\caption{Influence of grid spacing, displacement plot, case VMKH2 ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=$ 1 m, 5 m, $\Delta t=0.0001$ s, $h=0.05E$, $A=0.1$ mm)}
\label{fig:VMKH2disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{itemize}
\item[] \textbf{propagation through a dissipative elasticplastic material alters significantly the appearance of the input signal. All plots display significant wave attenuation/distortion, while final unrecoverable displacements are produced by soil plastifications (Figures \ref{fig:VMKH1disp}\ref{fig:VMKH2disp}(a)). Steady irreversible deformations are associated with prominent static components (at nil frequency) in the Fourier amplitude spectrum (Figures \ref{fig:VMKH1disp}\ref{fig:VMKH2disp}(c)), not present in the input Ormsby wavelet (Figure \ref{fig:ormsby_b});}
\item[] \textbf{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 $\Delta x$ values have been used in this subsection for illustrative purposes, whereas EM/PM plots have been deemed not necessary;}
\item[] \textbf{the influence of $\Delta x$ seems slightly magnified when lower $h$ values, and thus lower elasticplastic stiffness, are set (see Figure \ref{fig:VMKH2disp}). It is indeed not surprising that wave propagation in softer media may be more strictly affected by space discretization, as in linear problems. However, it should be noted that $\Delta x$ mainly influences the final irreversible displacement (Figure \ref{fig:VMKH2disp}(b)(c)), which leads to presume substantial interplay of grid effects and constitutive time integration;}
\item[] \textbf{since the effects of $\Delta x$ reduction are quite small in both time and frequency domains (for a given $\Delta t$), there is no strong motivation to suggest $\Delta x=V_s/20f_{max}$. $\Delta x=V_s/10f_{max}=\Delta x_{std}$ should be actually appropriate in common practical situations, as long as no soil failure mechanisms are triggered  as e.g. in seismic slope stability problems \citep{diprisco_pisano2012}. The occurrence of soil failure may introduce additional discretization requirements for an accurate kinematic representation of collapse.}
\end{itemize}
\begin{figure}[h!]
\centering
\subfigure[$\Delta x$=5 m]
{%
\includegraphics[width=0.48\textwidth]{figures/vmlkh8node_dxss_case2_ss01_cut.png}%
\label{fig:vmlkh8node_dxss_case2_ss01}%
}%
\
\subfigure[$\Delta x$=1 m]
{%
\includegraphics[width=0.48\textwidth]{figures/vmlkh8node_dxss_case2_ss02_cut.png}%
\label{fig:vmlkh8node_dxss_case2_ss02}%
}%
\caption{Influence of grid spacing, shear stressstrain response at the bottom of the VMKH sublayer, cases VMKH1 ($h=0.5E$) and VMKH2 ($h=0.05E$), ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta t=$ 0.0001 s, $A=0.1$ mm)}
\label{fig:VMKH1_2loop}
\end{figure}
In addition, Figure \ref{fig:VMKH1_2loop} illustrates the shear stressstrain VMKH response at the deepest
integration (Gauss) point of the VMKH sublayer. The material response is bilinear (elastic and elasticplastic), with the elastic stiffness
recovered upon stress reversal until new yielding occurs \citep{lemaitre1990}. As mentioned above, the observable (small) difference in stressstrain response at different $\Delta 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.
%%%%%%%%%%
%%%%%%%%%%
\begin{figure}[h!]
\centering
\subfigure[Displacement time history, $h=0.5E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_19a.pdf}%
\label{fig:vmlkh8node_dt_case3_dis2}%
}%
\
\subfigure[Displacement time history, $h=0.05E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_19b.pdf}%
\label{fig:vmlkh8node_dt_case4_dis2}%
}%
\
\subfigure[Displacement time history, $h=0.01E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_19c.pdf}%
\label{fig:vmlkh8node_dt_case6_dis2}%
}%
\
\subfigure[EM/PM misfits (ref. solution: $\Delta t=$ 0.0001 s)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig19.pdf}%
\label{fig:vmlkh8node_dt_case346_dis2_mis}%
}%
\caption{Influence of timestep size, displacement plot, cases VMKH3 ($h=0.5E$), VMKH4 ($h=0.05E$) and VMKH5 ($h=0.01E$) ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=5$ m, $\Delta t=$ 0.0002, 0.0005, 0.001 s, $A=0.1$ mm)}
\label{fig:VMKH3_4_5disp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textbf{The influence of the timestep size is illustrated for cases VMKH35 (Table \ref{tab:VMKH_analyses}) in Figures \ref{fig:VMKH3_4_5disp}\ref{fig:VMKH3_4_5loop}, encompassing three $h$ values ($0.5E$, $0.05E$ and $0.01E$) and also including EM/PM plots (Figure \ref{fig:VMKH3_4_5disp}(d)). In the lack of analytical solutions, misfits have been determined on the basis of a ``sufficiently accurate'' reference solution, here obtained numerically by setting $\Delta t=\Delta t_{std}/5=$ 0.0001 s. For a relatively small input amplitude ($A=$ 0.1 mm), convergence seems overall quite fast, and even $\Delta t=\Delta t_{std}$ results in both EM and PM values lower or close to 10\% (in combination with $\Delta x=\Delta x_{std}$). This inference is further corroborated by the shear stressstrain response at the bottom of the VMKH sublayer (Figure \ref{fig:VMKH3_4_5loop}), exhibiting little sensitiveness to the timestep size.
Some additional comments stem from the EM/PM plots in Figure \ref{fig:VMKH3_4_5disp}(d):}
\begin{itemize}
\item[]\textbf{ at variance with the previous elastic cases, envelope (EM) and phase (PM) misfits are quantitatively quite different (EM$>$PM);}
\item[] \textbf{EM/PM trends do not depend monotonically on the hardening parameter $h$. For $\Delta t=$ 0.0002 s, the EM/PM values at $h=0.05E$ are indeed larger than the those obtained for $h=0.5E$ and $h=0.01E$}.
\end{itemize}
\textbf{Both findings are likely to relate to the influence of time discretization on the residual displacement, which is larger than on other response features. In fact, variations in the accumulated displacement mainly affect the envelope of the output signal, not its phase attributes. The counterintuitive nonmonotonic 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. }
%%%%%%%%%%
%%%%%%%%%%
\begin{figure}[h!]
\centering
\subfigure[$\Delta t$=0.0002 s]
{
\includegraphics[width=0.9\textwidth]{figures/vmlkh8node_dt_ss_case2_ss01.png}
\label{fig:vmlkh8node_dt_ss_case2_ss01}
}
\
\subfigure[$\Delta t$=0.0005 s]
{
\includegraphics[width=0.9\textwidth]{figures/vmlkh8node_dt_ss_case2_ss02.png}
\label{fig:vmlkh8node_dt_ss_case2_ss02}
}
\
\subfigure[$\Delta t$=0.001 s]
{
\includegraphics[width=0.9\textwidth]{figures/vmlkh8node_dt_ss_case2_ss03.png}
\label{fig:vmlkh8node_dt_ss_case2_ss03}
}
\caption{Influence of timestep size, shear stressstrain response at the bottom of the VMKH sublayer, cases VMKH3 ($h=0.5E$), VMKH4 ($h=0.05E$) and VMKH5 ($h=0.01E$), ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=5$ m, $A=0.1$ mm)}
\label{fig:VMKH3_4_5loop}
\end{figure}
%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Influence of input motion amplitude}
%
In nonlinear 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 nonlinearity mobilized and the resulting local stiffness, in turn affecting the requirements for accurate constitutive integration.
In Figure \ref{fig:VMKH6_7disp}, the parametric study in Figures \ref{fig:VMKH1disp}\ref{fig:VMKH2disp} is replicated for a higher input amplitude ($A=$ 1 mm) and the same two different $h$ values (cases VMKH67 in Table \ref{tab:VMKH_analyses}).\textbf{ 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 effects and constitutive integration still apply to this case.}
The discussion on the influence of $\Delta t$ at higher input amplitude refers to Figures \ref{fig:VMKH8_9_10disp}\ref{fig:VMKH8_9_10loop}, illustrating the results obtained for $\Delta x=\Delta x_{std}$ and $h$ equal to $0.5E$, $0.05E$ and $0.01E$ (cases VMKH810 in Table \ref{tab:VMKH_analyses}); \textbf{again, EM/PM plots comes from the numerical reference solution corresponding with $\Delta t=\Delta t_{std}/5=$ 0.0001 s.}
The comparison
of Figures \ref{fig:VMKH6_7disp} and \ref{fig:VMKH8_9_10disp} suggests that, even with a much larger input amplitude, $\Delta x=\Delta x_{std}$ is still an appropriate grid spacing for elasticplastic problems, as long as $\Delta t$ is substantially reduced to comply with (explicit) constitutive integration requirements. \textbf{This inference is supported by the following observations:}
\begin{itemize}
\item[] \textbf{$\Delta t$ affects not only the residual component of displacement time histories (as in Figure \ref{fig:VMKH6_7disp}), but also their maximum/minimum transient values  i.e. the numerical representation of plastic dissipation. This is clearly visible in Figure \ref{fig:VMKH8_9_10disp}(a);}
\item[] \textbf{EM/PM values are in general higher at larger input amplitude (Figure \ref{fig:VMKH8_9_10disp}(d)), and experience a slower decrease as $\Delta t$ is reduced (still depending on the specific $h$ value);}
\item[] \textbf{the shear stressstrain loops in Figure \ref{fig:VMKH8_9_10loop} show how inaccurate the simulated constitutive response can be when $\Delta t$ is too large (e.g. $\Delta=$ 0.001 s) and substantial plastic degradation of material stiffness takes place (see the case $h=0.01E$).}
\end{itemize}
\textbf{This set of results suggests that $\Delta t$ should be at least in the order of $\Delta x/20V_s$ for acceptable constitutive integration and overall accuracy in elasticplastic wave simulations. However, such a conclusion seems quite heuristic, as it may be altered by the use of different material models (see next section) and stresspoint algorithms.}
%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history, $h=0.5E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_21a.pdf}%
\label{fig:vmlkh8node_dx_case1_dis2}%
}%
\
\subfigure[Displacement time history, $h=0.05E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_21b.pdf}%
\label{fig:vmlkh8node_dx_case2_dis2}%
}%
\
\caption{Influence of grid spacing, displacement plot, cases VMKH6 ($h=0.5E$) and VMKH7 ($h=0.05E$) ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=$1, 5 m, $\Delta t=0.0001{\rm s}$, $A=1$ mm)}
\label{fig:VMKH6_7disp}
\end{figure}
%%%%%%%%%%
%\clearpage
%%%%%%%%%%
\begin{figure}[!htb]
\centering
\
\subfigure[Displacement time history, $h=0.5E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_22a.pdf}%
\label{fig:vmlkh8node_dt_case1_dis2}%
}%
\
\subfigure[Displacement time history, $h=0.05E$ (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_22b.pdf}%
\label{fig:vmlkh8node_dt_case2_dis2}%
}%
\
\subfigure[Displacement time history, $h=0.01E$ (2.23.8 s) ]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_22c.pdf}%
\label{fig:vmlkh8node_dt_case5_dis2}%
}%
\
\subfigure[EM/PM misfits (ref. solution: $\Delta t=$ 0.0001 s)]
{%
\includegraphics[width=0.5\textwidth]{figures/empm_fig22.pdf}%
\label{fig:vmlkh8node_dt_case125_dis2_mis}%
}%
\caption{Influence of timestep size, displacement plot, cases VMKH8 ($h=0.5E$), VMKH9 ($h=0.05E$) and VMKH10 ($h=0.01E$) ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=5$ m, $\Delta t=$ 0.0002, 0.0005, 0.001 s, $A=1$ mm)}
\label{fig:VMKH8_9_10disp}
\end{figure}
%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[$\Delta t$=0.0002 s]
{
\includegraphics[width=0.9\textwidth]{figures/vmlkh8node_dt_ss_case1_ss01.png}
\label{fig:vmlkh8node_dt_ss_case1_ss01}
}
\
\subfigure[$\Delta t$=0.0005 s]
{
\includegraphics[width=0.9\textwidth]{figures/vmlkh8node_dt_ss_case1_ss02.png}
\label{fig:vmlkh8node_dt_ss_case1_ss02}
}
\
\subfigure[$\Delta t$=0.001 s]
{
\includegraphics[width=0.9\textwidth]{figures/vmlkh8node_dt_ss_case1_ss03.png}
\label{fig:vmlkh8node_dt_ss_case1_ss03}
}
\caption{Influence of timestep size, shear stressstrain response at the bottom of the VMKH sublayer, cases VMKH8 ($h=0.5E$), VMKH9 ($h=0.05E$) and VMKH10 ($h=0.01E$) ($\Delta x_{std}=5$ m, $\Delta t_{std}=0.0005$ s, $\Delta x=5$ m, $A=1$ mm)}
\label{fig:VMKH8_9_10loop}
\end{figure}
%\clearpage
%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{PBS model}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Model parameters and parametric analysis}
The influence of space/time discretization is now explored in combination with the nonlinear
PBS soil model introduced in Section \ref{sec:PBS} \citep{pisano2014}. As in real geomaterials, the PBS model features an
elasticplastic response since the very onset of loading (vanishing yield
locus), with the stiffness smoothly evolving from smallstrain 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 nonlinear PBS soil resting on a 400 m elastic sublayer. 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 \citep{pisano2014} (same elastic parameters for both the PBS and the elastic sublayers):
\begin{itemize}
\item[] $\rho=$ 2000 kg/m\textsuperscript{3}, $E=$ 1.3 GPa and $\nu=$ 0.3,
implying an elastic shear wave velocity $V_s=$ 500 m/s ;
\item[] shear strength parameter: $M=$ 1.2, corresponding with friction angle equal to 30 deg under triaxial compression;
\item[] dilatancy parameters: $k_d=$ 0.0 and $\xi=$ 0.0\footnote{Soil volume
changes under shear loading have been inhibited for the sake of simplicity. This
aspect would further affect the overall stiffness of the soil layer and require
additional parametric analyses.};
\item[] hardening parameters: $h=$ 300 and $m=$ 1.
\end{itemize}
The list of PBS simulations is reported in Table
\ref{tab:PJ_analyses}, while the next figures will also point out the good performance of the PBS model
in reproducing the cyclic soil behavior.
\begin{table}[h!]
\centering
\begin{tabular}{ccccccc}
\hline\hline
case \# & $\Delta x_{std}$ [m] & $\Delta t_{std}$ [s] & $\Delta x$ [m] & $\Delta t$ [s] & $A$ [mm]\\
\hline
PBS1 & 2.5 & 0.0005& 0.5, 2.5 & 0.0001 & 1 \\
PBS2 & 2.5 & 0.0005& 0.1, 0.5, 1 & 0.00002 & 1 \\
\hdashline
PBS3 & 2.5 & 0.0005& 2.5 & 0.0002, 0.0005, 0.001 & 1 \\
PBS4 & 2.5 & 0.0005& 2.5 & 0.00001, 0.00002, 0.0001 & 1 \\
\hline\hline
\end{tabular}
\caption{List of PBS simulations}
\label{tab:PJ_analyses}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Influence of grid spacing and timestep size }
%
Most of the issues observed in VMKH simulations appear to be magnified by the more complex PBS model. A summary of the main inferences drawn on the basis of Figures \ref{fig:PJ1disp}\ref{fig:PJ4loop} is provided here below:
\begin{itemize}
\item[] grid spacing turns out to be influential again (Figures \ref{fig:PJ1disp} and \ref{fig:PJ2disp}), 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;
\item[] \textbf{as in VMKH simulations, grid spacing mainly affects residual displacements. This is clearly shown by the EM/PM plots in Figure \ref{fig:PJ2disp}(b), where EM errors larger than 10\% arise even for a very small timestep size ($\Delta t =\Delta 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 $\Delta x$ $4\div5$ times smaller than $\Delta x_{std}=V/10f_{max}$;}
\item[] \textbf{the combination of explicit constitutive integration and high nonlinearity makes timestepping effects quite prominent, as is shown by Figures \ref{fig:PJ3disp} and \ref{fig:PJ3loop}. Further, Figures \ref{fig:PJ4displacement} leads to conclude that $\Delta t=\Delta t_{std}/50$ may be needed to obtain EM errors lower than 10\% (Figures \ref{fig:PJ4displacement}\ref{fig:PJ4loop}). Apparently,
analysts have to compromise on accuracy and computational costs in these situations;}
\item[] as expected, the shear stressstrain cycles in Figures \ref{fig:PJ1loop} and \ref{fig:PJ3loop} show that the sensitivity
to discretization builds up as increasing nonlinearity is mobilized. \textbf{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.}
\end{itemize}
%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (0.04.0 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_24a.pdf}%
\label{fig:pisano_dx_case1_dis1}%
}%
\
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_24b.pdf}%
\label{fig:pisano_dx_case1_dis2}%
}%
\caption{Influence of grid spacing, displacement plot, case PBS1 ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 0.5, 2.5 m, $\Delta t=$ 0.0001 s, $A=$ 1 mm)}
\label{fig:PJ1disp}
\end{figure}
%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[At the top of the layer]
{%
\includegraphics[width=0.48\textwidth]{figures/pisano_dx_case1_ss01.png}%
\label{fig:pisano_dx_case1_ss01}%
}%
\
\subfigure[At the bottom of the layer]
{%
\includegraphics[width=0.48\textwidth]{figures/pisano_dx_case1_ss02.png}%
\label{fig:pisano_dx_case1_ss02}%
}%
\caption{Influence of grid spacing, shear stressstrain response in the PBS sublayer, case PBS1 ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 0.5, 2.5 m, $\Delta t=$ 0.0001 s, $A=$ 1 mm)}
\label{fig:PJ1loop}
\end{figure}
%%%%%%%%%%%
%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{
\includegraphics[width=0.45\textwidth]{figures/fig_26a.pdf}
\label{fig:pisano_dx_case2_dis2}
}
\subfigure[EM/PM misfits (ref. solution: $\Delta x=$ m \textbf{KOHEI?})]
{
\includegraphics[width=0.48\textwidth]{figures/empm_fig26.pdf}
\label{fig:pisano_dx_case2_dis2_mis}
}
\caption{Influence of grid spacing, displacement plot, case PBS2 ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 0.1, 0.5, 1 m, $\Delta t=$ 0.00002 s, $A=$ 1 mm)}
\label{fig:PJ2disp}
\end{figure}
%%%%%%%%%%%
%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (0.04.0 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_27a.pdf}%
\label{fig:pisano_dt_case1_dis1}%
}%
\
\subfigure[Displacement time history (2.23.8 s)]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_27b.pdf}%
\label{fig:pisano_dt_case1_dis2}%
}%
\caption{Influence of timestep size, displacement plot, case PBS3 ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 2.5 m, $\Delta t=$ 0.0002, 0.0005, 0.001 s, $A=$ 1 mm)}
\label{fig:PJ3disp}
\end{figure}
%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[At the top of the layer]
{%
\includegraphics[width=0.9\textwidth]{figures/pisano_dt_case1_ss01.png}%
\label{fig:pisano_dt_case1_ss01}%
}%
\
\subfigure[At the bottom of the layer]
{%
\includegraphics[width=0.9\textwidth]{figures/pisano_dt_case1_ss02.png}%
\label{fig:pisano_dt_case1_ss02}%
}%
\caption{Influence of timestep size, shear stressstrain response in the PBS sublayer, case PBS3 ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 2.5 m, $\Delta t=$ 0.0002, 0.0005, 0.001 s, $A=$ 1 mm)}
\label{fig:PJ3loop}
\end{figure}
%
%
%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[Displacement time history (2.23.8 s)]
{
\includegraphics[width=0.45\textwidth]{figures/fig_29a.pdf}
\label{fig:pisano_dt_case2_dis2}
}
\subfigure[EM/PM misfits (ref. solution: $\Delta t=$ s \textbf{KOHEI?})]
{
\includegraphics[width=0.48\textwidth]{figures/empm_fig2729.pdf}
\label{fig:pisano_dt_case2_dis2_mis}
}
\caption{Influence of timestep size, displacement plot, case PBS4 ($\Delta x_{std}=$ 2.5. m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 2.5 m, $\Delta t=$ 0.00001, 0.00002,0.0001 s, $A=$ 1 mm)}
\label{fig:PJ4displacement}
\end{figure}
%
%%%%%%%%%%%
%
\begin{figure}[h!]
\centering
\includegraphics[width=0.48\textwidth]{figures/pisano_dt_case2_ss01.png}
\caption{Influence of timestep size, shear stressstrain response at the bottom of the PBS sublayer, case PBS4 ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 2.5 m, $\Delta t=$ 0.00001, 0.00002,0.0001 s, $A=$ 1 mm)}
\label{fig:PJ4loop}
\end{figure}
Since displacement components result from strains through spatial integration,
the displacement performance can be wellpredicted on condition that strains are
accurately computed all along the soil domain. For the same reason, the
discretization requirements for displacement convergence are not uniform along
the soil deposit. Figures \ref{fig:grid_depth_disp} and \ref{fig:step_depth_disp}
illustrate in the time domain the displacements simulated at different depths in the
nonlinear sublayer (the vertical $x$ axis points upward  Figure \ref{fig:sketchIBVP}) and
at different $\Delta x$ and $\Delta 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 this
implies satisfactory accuracy in a lower number of nodes and integration points.
\textbf{Conversely, the close relationship between plastic strains and residual displacements has slender influence on acceleration components. In this respect, Figures \ref{fig:grid_depth_acc} and \ref{fig:step_depth_acc} show that, as long as reasonable grid spacing is set (perhaps in the order of $\Delta x_{std}/2=V_s/20f_{max}$), the
sensitivity of acceleration components to $\Delta t$ is practically negligible (or definitely weaker than for residual displacements)}.
%\clearpage
%%%%%%%%%%%
%
\begin{figure}[h!]
\centering
\subfigure[x=500 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_31a.pdf}%
\label{fig:pisano_z2_case1_dis1}%
}%
\
\subfigure[x=480m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_31b.pdf}%
\label{fig:pisano_z2_case1_dis2}%
}%
\
\subfigure[x=460 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_31c.pdf}%
\label{fig:pisano_z2_case1_dis3}%
}%
\
\subfigure[x=440 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_31d.pdf}%
\label{fig:pisano_z2_case1_dis4}%
}%
\caption{Influence of grid spacing at different locations along the PBS layer, displacement plot, PBS2 case ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 0.1, 0.5, 1 m, $\Delta t=$ 0.00002 s, $A=$ 1 mm) }
\label{fig:grid_depth_disp}
\end{figure}
%%%%%%%%%%%
%\clearpage
%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[x=500 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_32a.pdf}%
\label{fig:pisano_z2_case2_dis1}%
}%
\
\subfigure[x=480 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_32b.pdf}%
\label{fig:pisano_z2_case2_dis2}%
}%
\
\subfigure[x=460 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_32c.pdf}%
\label{fig:pisano_z2_case2_dis3}%
}%
\
\subfigure[x=440 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_32d.pdf}%
\label{fig:pisano_z2_case2_dis4}%
}%
\caption{Influence of timestep size at different locations along the PBS layer, displacement plot, PBS4 case ($\Delta x_{std}=$ 2.5. m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 2.5 m, $\Delta t=$ 0.00001, 0.00002, 0.0001 s, $A=$ 1 mm)}
\label{fig:step_depth_disp}
\end{figure}
%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[x=500 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_33a.pdf}%
\label{fig:pisano_z2_case1_acc1}%
}%
\
\subfigure[x=480 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_33b.pdf}%
\label{fig:pisano_z2_case1_acc2}%
}%
\
\subfigure[x=460 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_33c.pdf}%
\label{fig:pisano_z2_case1_acc3}%
}%
\
\subfigure[x=440 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_33d.pdf}%
\label{fig:pisano_z2_case1_acc4}%
}%
\caption{Influence of grid spacing at different locations along the PBS layer, acceleration plot, PBS2 case ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 0.1, 0.5, 1 m, $\Delta t=$ 0.00002 s, $A=$ 1 mm) }
\label{fig:grid_depth_acc}
\end{figure}
%%%%%%%%%%%
\begin{figure}[!htb]
\centering
\subfigure[x=500 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_34a.pdf}%
\label{fig:pisano_z2_case2_acc1}%
}%
% \
\subfigure[x=480 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_34b.pdf}%
\label{fig:pisano_z2_case2_acc2}%
}%
\
\subfigure[x=460 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_34c.pdf}%
\label{fig:pisano_z2_case2_acc3}%
}%
% \
\subfigure[x=440 m]
{%
\includegraphics[width=0.48\textwidth]{figures/fig_34d.pdf}%
\label{fig:pisano_z2_case2_acc4}%
}%
\caption{Influence of timestep size at different locations along the PBS layer, acceleration plot, PBS4 case ($\Delta x_{std}=$ 2.5 m, $\Delta t_{std}=$ 0.0005 s, $\Delta x=$ 2.5 m, $\Delta t=$ 0.00001, 0.00002, 0.0001 s, $A=$ 1 mm)}
\label{fig:step_depth_acc}
\end{figure}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
\section{Concluding remarks}
Previously established criteria for space/time discretization in wave
propagation FE simulations have been reappraised and critically discussed, in order to provide more solid ground to verification procedures in Computational Dynamics. The 1D propagation of seismic shear waves (Ormsby wavelets) through
both linear and nonlinear (elasticplastic) media has been numerically simulated, with focus on capturing highfrequency motion and exploring the relationship between material response and discretization requirements. After initial linear computations, two different material models
(referred to as VMKH and PBS) have been used at increasing level of complexity.
The main conclusions inferred are hereafter
summarized:
\paragraph*{Elastic simulations} Setting grid spacing and timestep size as per standard
rules ($\Delta x_{std}=V_s/10f_{max}$ and $\Delta t_{std}=\Delta t/V_s$) has
proven not to be always appropriate, especially to reproduce highfrequency
motion components (this can be clearly visualized in the Fourier phase plane).
When linear elements (8node bricks) are used, $\Delta x\approx \Delta
x_{std}/2$ and $\Delta t\approx\Delta t_{std}/2$ seem to ensure sufficient
accuracy over the whole frequency range (both in amplitude and phase); higher
order elements (e.g. 27node bricks), will allow the use of $\Delta x=\Delta
x_{std}$ still in combination with $\Delta t\approx \Delta t_{std}/2$. Preserving accuracy in simulations including large domains and/or time intervals seems intrinsically more difficult, since attenuation/dispersion phenomena are cumulative.
\paragraph*{Elasticplastic simulations} Conclusive criteria for elasticplastic problems can be hardly established, as space/time discretization also interferes with the integration of nonlinear constitutive equations. In this respect, different outcomes may be found depending on (i) the kind and amount of nonlinearity
associated with the material model (stiffness variations during straining), (ii)
stresspoint integration algorithm (e.g. explicit or implicit), (iii) the amplitude of the input motion, affecting the amount of material plasticity mobilized during wave motion. The experience gained through the use of the PBS model (explicitly integrated in 8node brick elements) suggests that $\Delta
x_{std}=V_s/10f_{max}$ and $\Delta t_{std}=\Delta t/10V_s$ may need to be reduced by factors up to $4\div5$ and $500$, respectively, in the presence of strong input motions and severe stiffness variations. Importantly, these conclusions also depend on what output component is being considered and where within the computational domain.
\textbf{It is also worth remarking that the present study may not be regarded as conclusive, especially when it comes to nonlinear wave problems. There are in fact several aspects that will deserve in the future further consideration, such as the implications of using higherorder finite elements. The same comment also applies to geometrical effects (e.g. scattering) in 2D/3D problems, whose influence on discretization criteria for elasticplastic simulations would be per se a whole research topic.}
In the authors' hope, this study will reopen the discussion on a subject that is as relevant as often overlooked. Nowadays, modern computational platforms are increasingly supporting engineering analysis and decisionmaking, that should always be supported by robust verification procedures. As shown for nonlinear wave problems, combining knowledge in computational methods and material mechanics seems to be the way towards more reliable numerical simulations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgements}
%
The Authors wish to acknowledge the financial support from the USNRC, USDOE and Shimizu Inc. Corporation (Japan).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\baselinestretch}{1.0}
%\small\normalsize % trick from Kopka and Daly p47.
\small % trick from Kopka and Daly p47.
%\bibliographystyle{chicnarm} % < MOD
%\bibliographystyle{chicago} % < MOD
\bibliographystyle{abbrvnat} % < MOD
%\bibliographystyle{plain} % < MOD
%
\bibliography{./ref_Watanabe_Pisano_Jeremic} % These is my BiBTeX reference list.
\end{document}