\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}\kern-0.20em\uppercase{E}\kern-0.20em\uppercase{I}}}



% 
% % to change spacing in subcaption, see
% % http://www.ctan.org/tex-archive/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 elastic-plastic 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) 754-9248. fax~(530) 754-7872,
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 
elastic-plastic materials.
The establishment  of such criteria is especially important when accurate prediction
of  high-frequency  motion  is  needed and/or in the presence of markedly non-linear
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 wide-band input motions in combination with both
linear  elastic and non-linear elastic-plastic material models. The blind use of
usual  rules  of thumb is shown to be sometimes debatable, and an effort is made
to   provide   improved  discretization  criteria.  Possible  pitfalls  of  wave
simulations  are  pointed  out  by 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, elastic-plastic, 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 socio-economical 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  non-linearities,  etc.
  \citep{argyrismlejnek1991, kramer1996, zienkiewicz_etal_1999}.

\end{itemize}

When  linear  (elastic)  wave  problems  are  considered,  either time-domain or
frequency-domain  solutions  may  be  sought, whereas time-domain approaches are
usually needed in the presence of non-linearities (constitutive or geometrical).
In  this  respect,  it  is  worth  remarking  that  most  of  the interest in
seismology, applied geophysics and earthquake
engineering  is nowadays on non-linear 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                              man-made                              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
earthquake-soil-structure   interaction  can  be  only  performed  by  means  of
high-fidelity  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 up-to-date contribution to  the matter, also accounting for
the      increased     importance    assumed  in  recent  years  by  non-linear,
elastic-plastic 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 time-frequency 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  non-linearities 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 non-linearities. 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 cross-section can be regarded as symmetry planes and the soil deposit undergoes a ``double  plane-strain''  deformation, with both  horizontal direct
strains  prevented by symmetry \citep{prevost1989, borja1999}. As a consequence,
 unknown  variables  only depend on time and the vertical spatial coordinate
(the problem is geometrically one-dimensional), though the stress state is still
multi-axial  \citep{diprisco_pisano2012}.  The  initial-boundary  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, 8-27 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 elastic-plastic 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-, 8-27, and 27-node 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   well-known   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  time-steps   $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 second-order 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 well-known 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 linear-elastic to advanced elastic-plastic 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
elastic-plastic   von   Mises  model  with  linear  kinematic  hardening  
\citep{lemaitre1990, Jeremic2004d} and
(iii) the bounding surface elastic-plastic 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 time-step.


%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Elastic-plastic: von Mises kinematic hardening (VMKH) model} 
\label{VMKHmodel}

The  relationship among discretization, accuracy and material non-linearity will
be  first  explored  through  an elastic-plastic von Mises constitutive law with
linear  kinematic  hardening,  of the same kind described in \cite{lemaitre1990, Jeremic2004d}.

The  VMKH  model  is  very  well-known  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{Non-linear 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 post-yielding
stiffness  is  constant,  not  strain-dependent:  this  implies  an  unrealistic
unbounded  strength,  but allows to identify the elastic-plastic shear stiffness
with no ambiguity. Only four constitutive parameters need to be set:

\begin{itemize}

  \item[--]  two  elastic  parameters  --  $E$ and $\nu$ -- for the pre-yielding
  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 post-yielding
  (elastic-plastic) stiffness.

\end{itemize}

The  results  from  wave  propagation  analyses  will  clearly  show the kind of
stress-strain response arising from the VMKH formulation.

%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Elastic-plastic: 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 non-linear hardening, implying a continuous transition from small-strain
  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  stress-strain 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 at--rest 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 plane-strain'' 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 low-cut, low-pass, high-cut and
high-pass  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 time-frequency 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 time-frequency 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 time-frequency 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 frequency-dependent misfit measures by projecting $\Delta E$ and $\Delta P$ onto one of the two domains. In particular, the following single-values 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}, open-source routines for misfit analysis are available at \url{www.nuquake.eu/Computer
Codes/} (TF-MISFITS package). Discretization effects in wave propagation simulations will be assessed in the following on the basis of EM and PM criteria, as previously done by a number of authors \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
well-known  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  8-node  bricks.}  and  time-step 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 time-step 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 time-steps
--  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 time-step size are
discussed separately for  the  sake  of clarity.
Since  the Real  ESSI program  is  based  on  a  displacement  FE  formulation, displacement
components  are  the  most  reliable  output;  however,  some consideration is also paid to accelerations,
post-calculated  through second-order 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) time-step  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}{c|cccccc}
\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 & 8-node\\
EL2 & 20& 5 & 0.005 & 2, 5, 10 & 0.002 & 8-node\\
EL3 & 50&  2& 0.002 & 0.8, 2, 4 & 0.002 & 8-node\\
EL4 & 50&  2& 0.002 & 0.8, 2, 4 & 0.001 & 8-node\\
EL5 & 20& 5 & 0.005 & 2, 5, 10 & 0.002 & 27-node\\
\hdashline
EL6 & 20  & 5 & 0.005 & 5 & 0.002, 0.005, 0.01 & 8-node \\
EL7 & 20 & 5 & 0.005 & 2 & 0.001, 0.002, 0.005 & 8-node \\
EL8 & 50  & 2 & 0.002 & 2 & 0.001, 0.002, 0.005 & 8-node\\
EL9 & 50  & 2 & 0.002 & 0.8 & 0.0005, 0.001, 0.002 & 8-node\\
EL10 & 20 &  5 & 0.005  & 5 & 0.002, 0.005, 0.01 & 27-node\\
\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   EL1--EL5  (Table
\ref{tab:elastic_analyses}) in terms of:
(a--b)  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.0-4.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.2--3.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, 8-node brick)}
    \label{fig:EL1disp}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
   \centering
     \subfigure[Displacement time history (2.2-3.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, 8-node brick)}
    \label{fig:EL2disp}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
   \centering
    \subfigure[Displacement time history (2.2-2.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, 8-node brick)}
    \label{fig:EL4disp}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
   \centering
    \subfigure[Displacement time history (2.2--3.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, 27-node 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 (a--b) 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, 8-node brick) and (c--d) 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, 8-node 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 (a--b) 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, 8-node brick) and (c--d) 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, 8-node 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 time-step improves the numerical solution mostly in terms of Fourier phase, not amplitude (compares Figures \ref{fig:EL1disp}(c--d) and \ref{fig:EL2disp}(b--c)). It may be generally stated that, when  $\Delta  x$ is not appropriate, reducing the time-step 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 low-order FEs (8-node 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 EM-PM 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 8-node brick elements, while  Figure  \ref{fig:EL5disp}  shows that ``ten nodes per
wavelength'' are still suitable when higher--order elements (here 27-node bricks\footnote{\textbf{For a given number of nodes per wavelength, the size $\Delta x$ of 27-node  elements along the propagation direction is double than for 8-node 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 low--order elements, more severe discretization requirements should be fulfilled if very accurate accelerations need to be computed.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
%\newpage
\subsubsection{Influence of time-step size}

For given grid spacings, the influence of $\Delta t$ has been studied by varying the time-step 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 high-frequency phase-difference 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 (EM-PM lower than 10\%) to most practical purposes (see Figures \ref{fig:EL6disp}-\ref{fig:EL10disp}(d))};

\item[--] \textbf{when 27-node 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.2--3.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   time-step 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, 8-node brick)}
    \label{fig:EL6disp}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
    \centering
    \subfigure[Displacement time history (2.2--3.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, 8-node brick)}
    \label{fig:EL7disp}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
   \centering
    \subfigure[Displacement time history (2.2--2.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   time-step 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, 8-node brick)}
    \label{fig:EL9disp}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
    \centering
    \subfigure[Displacement time history (2.2--3.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   time-step 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, 27-node brick)}
    \label{fig:EL10disp}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[!htb]
    \centering
     \
    \subfigure[Acceleration time history (2.2--3.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.2--3.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 time-step 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, 8-node 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, 8-node brick)}
    \label{fig:EL6EL7acc}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[!htb]
    \centering
    \subfigure[Acceleration time history (2.2--2.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.2--2.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 time-step 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, 8-node 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, 8-node brick)}
    \label{fig:EL8EL9acc}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\graphicspath{{./figures/}}
\begin{figure}[h!]
     \centering
    \subfigure[Displacement time history (0.0-13.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, 8-node 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 never-ending wave motion is established. It is thus interesting to check how discretization errors propagate in time at the free surface, as is shown in 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 high-frequency 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{Non-linear elastic-plastic wave simulations}

This   section  concerns discretization effects in   presence  of  material  non-linearity.  As  most commonly done in Geomechanics
\citep{zienkiewicz_etal_1999}, the non-linear cyclic response of geomaterials can be described in the framework of elasto-plasticity, 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  non-linear  problem  under  consideration  cannot be solved analytically.  Therefore, the quality of discretization settings may  only be assessed by evaluating the converging behavior of numerical solutions upon $\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  non-linear  computations is highly affected by the input  amplitude.  This governs the amount of non-linearity mobilized by wave  motion  and,  as a consequence, the accuracy of numerical solutions at varying discretization.
\end{itemize}

In  non-linear (elastic-plastic) problems, discretization  is not only responsible for the numerical representation of waves
 (dissipation,  dispersion,  stability), but also governs the accuracy of
 constitutive integration \citep{simo1998,borja2013}. For instance, changes  in  time-step size will affect the strain size driving the constitutive integration algorithm and, in turn, the final simulation results. 
%
This    dependence  of  the constitutive  response  (material model and constitutive
integration  algorithm) on  the dynamic step size precludes direct development of automatic
criteria for discretization. However, as tangent elastic-plastic response can be established  for any stress-strain combination, (lowest) elastic-plastic (shear)
stiffness   may   be   used   to   develop  suitable  discretization via Equation~\ref{eq_vs}.  Apparently,  this  approach  assumes  that  the stress-strain response   is   already  known,  as  is  not  the  case  when discretization is being set. This means that an iterative approach is in principle needed, whereby one will first
design discretization based on an estimate of the strain level, run the dynamic simulation, and record  the  actual  stress-strain  response. After few iterations, a stable discretization will be usually achieved.
%
%


In  this  study, VMKH and PM constitutive equations have been integrated via the standard forward Euler, explicit algorithm   \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 elastic-plastic computations of time-step sizes derived through elastic parameters and Equation \eqref{eq:step}, especially  in combination with explicit  stress-point algorithms.  For this reason, the following time marching rule may be regarded as an upper bound for non-linear 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 8-node brick elements are employed for non-linear 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 sub-layer resting on an elastic stratum (remaining 800 m).
At  the surface, a thin layer (5 m) of elastic material has been added to prevent numerical problems with very strong motions and the so-called whip effect.
The following constitutive
parameters  (see  Section  \ref{VMKHmodel})  have  been  set (same  elastic parameters for both the VMKH and the elastic sub-layers), with no algorithmic nor Rayleigh damping introduced in numerical computations.:
\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 sub-layers);
\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 elastic-plastic 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}{c|cccccc}
\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 time-step size}
%

The results in Figures~\ref{fig:VMKH1disp} and \ref{fig:VMKH2disp} exemplify the role played by space discretization in elastic-plastic simulations. These results have been obtained by employing a time-step smaller than $\Delta t_{std}$ (cases VMKH1--2 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.0--4.0 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_16a.pdf}%
        \label{fig:vmlkh-8node_dx_case3_dis1}%
    }%
    \
    \subfigure[Displacement time history (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_16b.pdf}%
        \label{fig:vmlkh-8node_dx_case3_dis2}%
    }%
    \
    \subfigure[Amplitude of displacement Fourier spectrum]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_16c.pdf}%
        \label{fig:vmlkh-8node_dx_case3_famp}%
    }%
     \
    \subfigure[Phase difference of displacement Fourier spectrum]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_16d.pdf}%
        \label{vmlkh-8node_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.0--4.0 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_17a.pdf}%
        \label{fig:vmlkh-8node_dx_case4_dis1}%
    }%
    \
    \subfigure[Displacement time history (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_17b.pdf}%
        \label{fig:vmlkh-8node_dx_case4_dis2}%
    }%
    \
    \subfigure[Amplitude of displacement Fourier spectrum]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_17c.pdf}%
        \label{fig:vmlkh-8node_dx_case4_famp}%
    }%
     \
    \subfigure[Phase difference of displacement Fourier spectrum]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_17d.pdf}%
        \label{vmlkh-8node_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 elastic-plastic 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 elastic-plastic 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/vmlkh-8node_dx-ss_case2_ss01_cut.png}%
        \label{fig:vmlkh-8node_dx-ss_case2_ss01}%
    }%
    \
    \subfigure[$\Delta x$=1 m]
    {%
        \includegraphics[width=0.48\textwidth]{figures/vmlkh-8node_dx-ss_case2_ss02_cut.png}%
        \label{fig:vmlkh-8node_dx-ss_case2_ss02}%
    }%
    \caption{Influence of grid spacing, shear stress-strain response at the bottom of the VMKH sub-layer, 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 stress-strain VMKH  response  at  the  deepest
integration  (Gauss) point of the VMKH sub-layer. The material response is bilinear (elastic and elastic-plastic), with the elastic stiffness 
recovered upon stress   reversal   until  new  yielding occurs \citep{lemaitre1990}.  As mentioned above, the  observable (small) difference  in stress-strain 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.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_19a.pdf}%
        \label{fig:vmlkh-8node_dt_case3_dis2}%
    }%
  \
    \subfigure[Displacement time history, $h=0.05E$  (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_19b.pdf}%
        \label{fig:vmlkh-8node_dt_case4_dis2}%
    }%
     \
    \subfigure[Displacement time history, $h=0.01E$  (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_19c.pdf}%
        \label{fig:vmlkh-8node_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:vmlkh-8node_dt_case346_dis2_mis}%
    }%
    \caption{Influence of time-step 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 time-step size is illustrated for cases VMKH3--5 (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 stress-strain response at the bottom of the VMKH sub-layer (Figure \ref{fig:VMKH3_4_5loop}), exhibiting little sensitiveness to the time-step 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 non-monotonic relationship between $h$ and displacement EM/PM  values has not been detected in the corresponding acceleration EM/PM plots (not reported here), due to the obvious lack of residual accelerations. }

%%%%%%%%%%
%%%%%%%%%%
\begin{figure}[h!]
    \centering
    \subfigure[$\Delta t$=0.0002 s]
    {
        \includegraphics[width=0.9\textwidth]{figures/vmlkh-8node_dt_ss_case2_ss01.png}
        \label{fig:vmlkh-8node_dt_ss_case2_ss01}
    }
    \
    \subfigure[$\Delta t$=0.0005 s]
    {
        \includegraphics[width=0.9\textwidth]{figures/vmlkh-8node_dt_ss_case2_ss02.png}
        \label{fig:vmlkh-8node_dt_ss_case2_ss02}
    }
    \
    \subfigure[$\Delta t$=0.001 s]
    {
        \includegraphics[width=0.9\textwidth]{figures/vmlkh-8node_dt_ss_case2_ss03.png}
        \label{fig:vmlkh-8node_dt_ss_case2_ss03}
    }
    \caption{Influence of time-step size, shear stress-strain response at the bottom of the VMKH sub-layer, 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    non-linear    problems,   it  is  hard to  draw  general conclusions  on the interaction between space/time discretization  and input amplitude. The latter governs the
amount  of soil non-linearity mobilized and the resulting local stiffness, in turn affecting the requirements for accurate constitutive integration.

In 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 VMKH6-7 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 VMKH8-10 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 elastic-plastic 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 stress-strain 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 elastic-plastic 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 stress-point algorithms.} 



%%%%%%%%%%
\begin{figure}[!htb]
     \centering
        \subfigure[Displacement time history, $h=0.5E$ (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_21a.pdf}%
        \label{fig:vmlkh-8node_dx_case1_dis2}%
    }%
 \
    \subfigure[Displacement time history, $h=0.05E$ (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_21b.pdf}%
        \label{fig:vmlkh-8node_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.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_22a.pdf}%
        \label{fig:vmlkh-8node_dt_case1_dis2}%
    }%
     \
    \subfigure[Displacement time history, $h=0.05E$ (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_22b.pdf}%
        \label{fig:vmlkh-8node_dt_case2_dis2}%
    }%
 \
    \subfigure[Displacement time history, $h=0.01E$ (2.2--3.8 s) ]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_22c.pdf}%
        \label{fig:vmlkh-8node_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:vmlkh-8node_dt_case125_dis2_mis}%
    }%
    \caption{Influence of time-step 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/vmlkh-8node_dt_ss_case1_ss01.png}
        \label{fig:vmlkh-8node_dt_ss_case1_ss01}
    }
    \
    \subfigure[$\Delta t$=0.0005 s]
    {
        \includegraphics[width=0.9\textwidth]{figures/vmlkh-8node_dt_ss_case1_ss02.png}
        \label{fig:vmlkh-8node_dt_ss_case1_ss02}
    }
    \
    \subfigure[$\Delta t$=0.001 s]
    {
        \includegraphics[width=0.9\textwidth]{figures/vmlkh-8node_dt_ss_case1_ss03.png}
        \label{fig:vmlkh-8node_dt_ss_case1_ss03}
    }
    \caption{Influence of time-step size, shear stress-strain response at the bottom of the VMKH sub-layer, 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 non-linear
PBS  soil  model introduced in Section \ref{sec:PBS} \citep{pisano2014}.  As in real geomaterials, the  PBS model   features  an
elastic-plastic  response  since  the  very  onset  of  loading (vanishing yield
locus),  with the stiffness smoothly evolving from small-strain elastic behavior to failure (nil stiffness).

The  results presented hereafter concern a $500$ m thick soil layer, whose upper 100 m
are  made  of  a  non-linear  PBS  soil  resting  on  a  400 m  elastic sub-layer. As done for the VMKH simulations, a thin layer (2.5 m) of elastic
material  has  been added to prevent numerical problems with very strong motions and the whip effect at the ground surface.
Input 1 with $A=1$ mm has been exclusively considered, along with the following set of PBS  parameters \citep{pisano2014} (same elastic parameters for both the PBS and the elastic sub-layers):

\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}{c|cccccc}
\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 time-step 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 time-step 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 non-linearity, 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 non-linearity makes time-stepping 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 stress-strain cycles in Figures \ref{fig:PJ1loop} and \ref{fig:PJ3loop} show that the sensitivity
to  discretization  builds up  as increasing non-linearity 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.0--4.0 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_24a.pdf}%
        \label{fig:pisano_dx_case1_dis1}%
    }%
    \
    \subfigure[Displacement time history (2.2--3.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 stress-strain response in the PBS sub-layer, 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.2--3.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.0--4.0 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_27a.pdf}%
        \label{fig:pisano_dt_case1_dis1}%
    }%
    \
    \subfigure[Displacement time history (2.2--3.8 s)]
    {%
        \includegraphics[width=0.48\textwidth]{figures/fig_27b.pdf}%
        \label{fig:pisano_dt_case1_dis2}%
    }%
    \caption{Influence of time-step 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 time-step size, shear stress-strain response in the PBS sub-layer, 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.2--3.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 time-step 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 time-step size, shear stress-strain response at the bottom of the PBS sub-layer, 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 well-predicted 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
non-linear  sub-layer  (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 time-step 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 time-step 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  re-appraised 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  non-linear  (elastic-plastic)  media  has  been  numerically simulated,  with  focus  on  capturing  high-frequency  motion and exploring the relationship   between material   response   and   discretization 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 time-step 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 high-frequency
motion  components  (this can be clearly visualized in the Fourier phase plane).
When   linear  elements  (8-node  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.  27-node  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*{Elastic-plastic simulations}     Conclusive     criteria    for elastic-plastic problems can be hardly established, as space/time discretization also interferes with  the  integration  of  non-linear  constitutive equations. In this respect, different  outcomes  may  be  found depending on (i) the kind and amount of non-linearity
associated with the material model (stiffness variations during straining), (ii)
stress-point 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    8-node    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 non-linear wave problems. There are  in fact several aspects that  will deserve in the future further consideration, such as the implications of using higher-order finite elements. The same comment also applies to geometrical effects (e.g. scattering) in 2D/3D problems, whose influence on discretization criteria for elastic-plastic 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 decision-making, that should always be supported by robust verification procedures.  As   shown  for  non-linear  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 US-NRC, US-DOE 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}

