\documentclass[12pt,number,sort&compress]{elsarticle}
% \documentclass[12pt]{elsarticle}


%\usepackage[nomarkers,figuresonly]{endfloat}
% \usepackage{endfloat}

\usepackage{lineno}


\usepackage{subfig}
\usepackage{multirow}
\usepackage{graphicx}
%\usepackage[autostyle]{csquotes}
\usepackage{csquotes}
\usepackage{leftidx}
\usepackage{amsmath}

\usepackage{xcolor}
\usepackage{textcomp}

\usepackage{listings}

\graphicspath{{./Figures/}}
 
%%%%  HYPERREF  HYPERREF  HYPERREF  HYPERREF  HYPERREF 
\usepackage[pdfauthor={Boris Jeremic}, 
            colorlinks=true,
            linkcolor=blue,
            citecolor=blue,
            urlcolor=blue]
            {hyperref}

\usepackage{enumitem}

\journal{Soil Dynamics and Earthquake Engineering}

% \bibliographystyle{elsarticle-harv}
%\bibliographystyle{elsarticle-num}
%\usepackage[numbers,round,colon,authoryear]{natbib}
% \usepackage[number,sort&compress]{natbib}

%% Water mark
%\usepackage{type1cm}
%\usepackage{eso-pic}
%\usepackage{color}

% \makeatletter
% \AddToShipoutPicture{%
%  \setlength{\@tempdimb}{.05\paperwidth}% X coord placement
%  \setlength{\@tempdimc}{.25\paperheight}% Y coord placement
%  \setlength{\unitlength}{1pt}
%  \put(\strip@pt\@tempdimb,\strip@pt\@tempdimc){
%   \makebox(-5,5){\rotatebox{90}{\textcolor[gray]{0.7}
%     {\fontsize{0.6cm}{0.5cm}\selectfont{DRAFT of Han et al, work in progress!}}}}
%  }
% }

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spacing REMOVE FOR official version
%
%
%baselinestretch
%
\renewcommand{\baselinestretch}{1.40}
\small\normalsize  % trick from Kopka and Daly p47.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\begin{frontmatter}

%\title{Wave  Potential  -  Domain Reduction Method for Earthquake Soil Structure
%Interaction Modeling and Simulation Under Inclined Incident Seismic Waves}

%\title{Modeling and Simulation of Earthquake Soil Structure
%Interaction Under Inclined Incident Seismic Waves}

\title{Modeling and Simulation of Earthquake Soil Structure
Interaction Excited by Inclined Seismic Waves}


%% or include affiliations in footnotes:
\author[mymainaddress]{Hexiang Wang}
% \ead{hexwang@ucdavis.edu}
\author[mymainaddress]{Han Yang}

%\author[mysecondaddress]{Yuan Feng}
\author[mymainaddress]{Yuan Feng}
%
%\author[myfourthaddress]{Fangbo Wang}
%
% \author[mysecondaddress]{David B McCallen}

\author[mymainaddress,mythirdaddress]{Boris Jeremi{\'c}\corref{mycorrespondingauthor}}
\ead{jeremic@ucdavis.edu}

\cortext[mycorrespondingauthor]{Corresponding author}

\address[mymainaddress]{Department   of  Civil  and  Environmental  Engineering,
University of California, Davis, CA, USA}

%\address[mysecondaddress]{TuSimple, San Diego, CA, USA}

%\address[myfourthaddress]{School of Civil Engineering, Tianjin University, Tianjin, China}
%
\address[mythirdaddress]{Environmental  and  Earth    Sciences  Area,  Lawrence 
Berkeley National Laboratory, Berkeley, CA, USA}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\begin{abstract}

Presented        is      an  application of  wave    potential formulation (WPF)
together  with    domain  reduction method (DRM) to modeling earthquake soil
structure  interaction  (ESSI)  behavior  in  horizontally  layered ground under
inclined incident seismic waves.
%
Wave  potential  formulation  is  used  to  develop  a  spatially   varying,
inclined  seismic  wave  field from incident Primary (P) and Secondary (S) waves
that propagate through layered ground.
%
Developed  seismic wave field is then used to develop effective forces for 
Domain Reduction Method that are then used for analyzing ESSI response of a soil
structure system.
%
Developed  methodology,  called  WPF-DRM,        is    verified  using  analytic
solution  for  a  free  field  response  of  layered ground subjected to inclined 
incident waves.
%

%Compared   with   conventional   substructure  method,  WPF-DRM  circumvents  the
%additional  difficulties  of  solving  foundation  wave scattering and impedance
%function.
%
Developed    WPF-DRM  methodology    is  illustrated through analysis of an ESSI
response of a deeply embedded structure, a small modular reactor (SMR) subjected
to  incident  S  wave  polarized  in  vertical  plane  (SV) with variation in 
inclinations and frequencies.
%
Presented example highlights the  influences  of  incident  wave inclination and
frequency on ESSI response of a deeply embedded structure. 
% %
% Modeling  uncertainty  in  commonly  adopted  simplification  of  1D  vertically
% propagating wave field is emphasized.
% %

\end{abstract}

\begin{keyword}
earthquake  soil  structure  interaction  \sep deeply embedded structure \sep inclined 
incident P/SV/SH waves \sep layered ground \sep small modular reactor
\end{keyword}

\end{frontmatter}

\linenumbers

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

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




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


\newpage
\setcounter{tocdepth}{4}
\tableofcontents




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Introduction}

It    has    been   recognized  that during earthquakes, inclined body waves and
surface waves have significant influence on a dynamic response of soil structure
interaction                             (SSI)                           
systems~\citep{Trifunac1982,Todorovska1990,Betti1993,Teisseyre2006,Hori2006}.
%
For    example,  incident secondary (S) waves, where soil/rock particles move in 
horizontal plane (SH), can cause torsional response of structures,
%
Similarly,  incident  primary  (P)  and secondary (S) waves, where particles move in
vertical  plane  (SV), can produce amplified rocking of structures, especially in
near-fault  regions  and  for  structures with large-plan dimensions or multiple
supports~\citep{Gicev2015}.
%
Earthquake  Soil  Structure  Interaction (ESSI) response due to inclined incident
seismic  waves  (i.e.,  P,  SH  and  SV  waves)  is  of  significant interest in
earthquake engineering.
%







The  Earthquake  Soil  Structure Interaction problem has been studied for a long
time.
%
Early  work    was   focused    on dividing an SSI problem into simpler problems
that were manageable with available methodology and tools.
%
%
Substructure   method~\citep{Wolf1985}   was   established  to  decompose  the
SSI  problem  into  three sub-problems: 

\begin{itemize}

  \item Free field seismic motion
  
  \item  Foundation    input    motion,  i.e. foundation  wave scattering and
           impedance function, and

  \item Superstructure dynamic response

\end{itemize}

%
\citet{Luco1982}  studied  dynamic  response  of SSI system under non-vertically
incident SH wave.
%
% Impedance functions for several types of foundations (rectangular,
% hemispherical, etc) and soil profiles have been developed
% \citep{Wong1985,Crouse1990, Fu2017}.
% 
SSI responses excited by in-plane wave (P, SV and Rayleigh waves) were presented
by~\citet{Todorovska1992,Todorovska1993}.
%
Effects  of site dynamic characteristics on SSI were systematically investigated
by~\citet{Liang2013,Liang2013b} for incident P, SV and SH waves.
%




Due to the limitation of substructure method and the complexity of SSI
problem, simplifications have been commonly made in many studies.
%
For  example,  underground is usually simplified as homogeneous half space or a
single homogeneous soil layer above the bedrock.
%
Rigid  foundation  with  specific  shape  is  typically  assumed,  in  order  to
calculate impedance
functions  and  wave  scattering.
%
This assumption could  lead to excessive scattering of
incident wave energy and underestimated structural response~\citep{Liang2013b}.
% %
% This  is especially true for some critical structures that are sensitive to high
% frequency, short wavelength shaking, e.g., Nuclear power plants (NPP).
% %




With  increase  in computer power, direct simulation of dynamic
SSI  using  finite  element  method  (FEM),  finite  difference method (FDM) and
boundary element method (BEM) becomes feasible.
%
\citet{Stamos1996}  studied  dynamic  response  of  infinitely  long  tunnels in
elastic  or  viscoelastic  half-space  under incident seismic waves by a special
direct BEM.
%
Translational,  torsional  and  rocking response of a building SSI system excited by
plane P, SV and SH wave using FDM was recently studied
by~\citet{Gicev2015,Gicev2016,Givcev2016b}.
%


For  direct  simulation  of  SSI,  effective input of inclined incident seismic
waves is of great importance.
%
Many  artificial  boundary types  have been developed by approximating the radiation
condition at the finite boundaries of SSI system~\citep{Lysmer1969,Wolf1985,Liu2006c,Baffet2012}.
%
%
% {\Large \bf NOTE: This was done much earlier by others
% see~\cite{Lysmer1969,Bamberger1988} and more recently by~\cite{Baffet2012} so we
% can put~\cite{Liu2006c}  however~\cite{Du2006} is in Chinese and journal will not allow
% that  for a reference, official language of the journal is English and they will
% insist on that, I know I tried once to put my name in Serbian Cyrilic...}
%
% ==== Reply by Hexiang === 
% Have added references of \cite{Lysmer1969} and \cite{Baffet2012}. The Chinese paper by r~\cite{Du2006} was removed. 
%
%
% Among   these,   the   viscous-spring   artificial   boundary   put  forward 
% by~\citet{Liu2006c} and \citet{Du2006} can efficiently simulate the radiation
% damping and the elastic resumption behavior of semi-infinite medium.
% %
%
%
Using developed viscous-spring  artificial  boundary,  various  SSI and
rock-structure interaction (RSI) systems excited by inclined
incident  plane  waves,  such as tunnels~\citep{Huang2016,Huang2017b},
powerhouse~\citep{Wang2019b} and underground large scale frame structure
(ULSFS)~\citep{Qiu2019} were analyzed.
%
In these previous studies, inclined plane waves are generally assumed to
occur in  homogeneous ground.
%
The  only  wave  reflection  and refraction is considered at the ground surface,
while  multiple  layers,  usually present in realistic geological settings, were
not considered.
%
%
It  is noted that modeling and simulation of  inclined wave propagation in
layered ground is more complicated because of multiple reflection, refraction,
reverberation and interference at both layer interfaces and ground surface.
%  % 
%  Furthermore, to the best knowledge of the authors', there is still no systematic
%  study  on  dynamic response of nuclear reactor deeply embedded in layered ground
%  subjected to inclined incident seismic waves.
%  %
%
%
Of interest is modeling and simulation of deeply embedded structures, that extend
over multiple soil layers in depth.
%
Inclined seismic wave field, propagating through a number of layers, will
interact with the embedded structure.
%
Embedement and stiffness of the structure will modify the seismic wave field.
%
This effect is usually called the kinematic interaction, and applies
for linear elastic SSI analysis, where kinematic and inertial interaction
effects can be separated, superimposed~\citep{Kramer96}.


% 
% These  spatially  varying,  inclined  ground motions can have more influences on
% deeply embedded structures than on common surface structures.
% %
% SSI of surface structures is controlled by inertia effects.
% %
% In  which  case,  simplification  of  vertically propagating wave field might be
% generally adequate.
% %
% Nevertheless,  because  of  the  large  embedement,  responses of deeply embedded
% structures  are  controlled  by  kinematic interaction with wave field, which is
% more inclined at deeper ground layers according to Snell's law.
% %








Presented  is  a methodology 
developed to investigate influence of
inclined  body  and  surface seismic wave on linear or nonlinear earthquake soil
structure interaction (ESSI) behavior of soil-structure systems.
%
Methodology    is    based    on    Wave    Potential  
Formulation (WPF) ~\citep{Thomson1950,Haskell1953} as well as Domain Reduction
Method (DRM) \citep{Bielak2001}.
% %
% The  wave  potential  formulation for inclined wave
% propagation   in   layered  media  was developed  by 
% \citep{Haskell1953}  
% %\citep{Haskell1953,Zhu2002}  
% for inclined wave
% propagation   in   layered  media  is  combined  with  domain  reduction  method
% \citep{Bielak2001} for effective input of inclined motions.
% %
%
%
Paper is organized as follows:
%
Brief presentation of Wave  Potential Formulation and Domain Reduction
Method is given in  section~\ref{sec:simulation_method}.
%  to  simulate dynamic SSI in layered ground
%subjected to inclined seismic excitations.
Combined Wave Potential Formulation and Domain Reduction Method (WPF-DRM) is then verified,
with select results presented in section~\ref{sec:numerical_study_verification}.
%
Following  that, dynamic response of a deeply embedded 
%($36$m) 
small
modular  reactor (SMR) under inclined incident SV wave at different frequencies
and inclinations is  analyzed and presented in
sections~\ref{Deeply_Embedded_Soil-Structure_Model}, 
\ref{SMR_Excited_with_Inclined_SV_Waves}
and~\ref{SMR_Excited_with_Variable_Frequenct_Inclined_SV_Waves}.
%
Findings are summarized in section~\ref{Summary}.
%
% Compared with conventional three-step substructure method, two-step WPF-DRM approach is simpler. 
%
% Only $3C$ free field motion is required to conduct SSI analysis. 
%
% Difficulties in computing wave scattering and impedance functions for arbitrary shaped, flexible foundations are avoided.  
%
% In addition, the formulation of DRM is based on the equivalence of effective earthquake faulting force, instead of principal of superposition. 
%
% It can be applied to model nonlinear SSI under $3C$ seismic excitations.
%













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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Wave Potential Formulation -- Domain Reduction Method}
\label{sec:simulation_method}

Presented WPF-DRM  methodology consists of three main steps:
%

\begin{enumerate}

\item 
Analytic development of free field ground motions for a layered half space,
excited by an incident, inclined plane wave. 
%
Development of this seismic wave field is relying on wave potential formulation
in frequency-wave number domain.
%
Time  domain spatially varying ground motions are then synthesized through
inverse Fourier transformation.
% %
% These  free field motions could not be directly imposed onto the boundary of SSI
% system.
% %

\item 
Development of the Effective Earthquake Forces, from DRM formulation, is then
performed using free field seismic motions developed in the previous step.
%
%Instead,  at the second step, dynamically consistent effective earthquake forces
%are  formulated  with  DRM and applied at a single layer of elements surrounding
%the SSI system.
%

\item 
Earthquake Soil Structure Interaction (ESSI) analysis of the soil-structure
system is then performed using effective earthquake forces that are applied to a
single layer of finite elements surrounding soil-structure system, so called
DRM layer.
%
The only waves that are radiated from the soil-structure system and exit the DRM layer are due to oscillations of the structure. 
%
These outgoing waves are absorbed by damping layers.

% The only waves radiated from the soil-structure system due to oscillations of the structure  and that exit the DRM layer are and are absorbed by damping layers.
%
%waves  are  dissipated  in the absorbing layer to keep
%spurious wave reflections within acceptable limits.
%%

\end{enumerate}




Sections~\ref{Wave_Potential_Formulation_for_Inclined_Incident_Waves_in_Layered_Media} 
and \ref{The_Domain_Reduction_Method} below provide details of Wave
Potential Formulation and Domain Reduction Method, respectively.








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Wave Potential Formulation for Inclined Incident Waves in Layered Media} 
\label{Wave_Potential_Formulation_for_Inclined_Incident_Waves_in_Layered_Media}



Considered  is  an inclined wave that propagates in the layered ground, as shown
in Figure~\ref{fig:2D_layered_ground}.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=13cm]{Figures/layering_ground_illustration.pdf}
\caption{Layered ground and free field inclined seismic motions.}
\label{fig:2D_layered_ground}
\end{figure}
%
%
There   are   $n$   layers,  with  layer  thickness  $d_m$,  density 
$\rho_{m}$, compressional wave velocity $\alpha_{m}$ and shear wave velocity
$\beta_{m}$ ($m=1,2,..,n$).
%
Focus of presented development is on inclined P and SV waves, and mode
conversion between them at layer boundaries.
%
Propagation of $SH$ wave is simpler as there is no mode
conversion, so these waves are left out of presented
considerations.
% 
%It is noted that propagation of SH waves through layered ground system is
%governed by Snell's law~\citet{Kramer96, Aki2002}.
%
%
It is noted that the  wave  potential  formulation  presented  below is general  and 
also applicable to incident SH wave~\citep{Haskell1953}.
%



Without  loss  of  generality,  incident waves is considered to be
monochromatic, single frequency, 
with angular frequency $w$ and horizontal phase velocity $c$.
%
For  incident  waves  with  arbitrary time signal and multiple frequencies, free
field motions can be Fourier synthesized from the monochromatic solutions.
%


According   to   Helmholtz   decomposition   theorem~\citep{Arfken1999},   the
displacement from wave propagation modeled using 
Equation~(\ref{eq_wave_equation}), for linear elastic
material with Lam\'{e} constants $\lambda$ and $\mu$, can be expressed with $P$ wave scalar potential $\phi$ and $S$ wave vector
potential $\boldsymbol{\Psi}$. 

% \begin{equation}
% \rho \ddot{\boldsymbol{u}} 
% = 
% (\lambda + 2\mu) \nabla \nabla \cdot \boldsymbol{u} 
% - 
% \mu \nabla \times \nabla \times \boldsymbol{u}  
% \label{eq_wave_equation}
% \end{equation}

%%% === indical notation ====
\begin{equation}
\rho u_{i, tt} = \mu u_{i, jj} + (\lambda + \mu) \mu_{j, ij}
\label{eq_wave_equation}
\end{equation}

This is shown in  Equation~\ref{eq_Helmholtz_decomposition}, where  $\phi$  is the
curl free part corresponding to volumetric deformation and $\Psi$ is divergence
free part corresponding to deviatoric deformation.
%
$e_{ijk}$ is Levi-Civita permutation symbol~\citep{Arfken1999}.   

%%% === original ====
% \begin{equation}
% \boldsymbol{u} = \nabla \phi +  \nabla \times \Psi
% \label{eq_Helmholtz_decomposition}
% \end{equation}

%%% === indical notation ====
\begin{equation}
u_i = \phi_{,i} +   \Psi_{k, j}e_{ijk}
\label{eq_Helmholtz_decomposition}
\end{equation}


Using this approach, the unknown displacements for the $m^{th}$ layer are
simplified into incident $P$  wave  potential  magnitude  $\phi^{'}_{m}$, 
reflected $P$ wave potential magnitude $\phi^{''}_{m}$, incident $SV$ wave
potential magnitude $\Psi^{'}_{m}$ and reflected $SV$ wave potential magnitude
$\Psi^{''}_{m}$, as shown in Equation~\ref{eq:potential_magnitude}.

\begin{equation}
\begin{gathered}
\phi_m 
= 
[\phi_{m}^{'}e^{ik(x-\gamma_{\alpha m}z)}
+
\phi_{m}^{''}e^{ik(x+\gamma_{\alpha m}z)}]e^{-iwt}\\
%
\Psi_m 
= 
[\Psi_{m}^{'}e^{ik(x-\gamma_{\beta m}z)}
+
\Psi_{m}^{''}e^{ik(x+\gamma_{\beta m}z)}]e^{-iwt}
\end{gathered}
\label{eq:potential_magnitude}
\end{equation}
%
In Equation~\ref{eq:potential_magnitude}, 
% he    $P$    and    $SV$    wave    potential   can   be   expressed   as   Eq.
% (\ref{eq:potential_magnitude}),  
the horizontal wave number $k$ is defined as $k = w/c$. 
%
The harmonic nature of the potential field is characterized by the time factor
$e^{-iwt}$.
%
The incident and reflected angles for $P$ and $SV$ wave are equal to $arccot\gamma_{\alpha_m}$ and $arccot\gamma_{\beta_m}$, where $\gamma_{\alpha_m}$ and $\gamma_{\beta_m}$ are given in Equation \ref{eq:angle_formulation}. 
%
%
%
% {\Large \bf NOTE: 
% OR more appropriately
%$\arccot{\gamma_{\alpha_m}}$ and $\arccot{\gamma_{\beta_m}}$ 
% $arccot{\gamma_{\alpha_m}}$ and $arccot{\gamma_{\beta_m}}$ }
%
%====== Reply from Hexiang ====== 
% changed $cot^{-1}\gamma_{\alpha_m}$ to $arccot{\gamma_{\beta_m}}$
%
%
%
%
%It will be understood and omitted hereafter.
%
% The incident and reflected angles for $P$ and $SV$ wave, $\gamma_{\alpha_m}$  and 
% $\gamma_{\beta_m}$ can be determined by Snell's law
% ~\citep{Kramer96, Aki2002} as Equation \ref{eq:angle_formulation}.

% \begin{equation}
% \begin{gathered}
% \gamma_{\alpha_m} = \sqrt{(c/\alpha_m)^2-1} \ \ when \ \ \alpha_m \leq c \\
% \gamma_{\alpha_m} = -i\sqrt{1-(c/\alpha_m)^2} \ \ when \ \ \alpha_m > c \\
% \gamma_{\beta_m} = \sqrt{(c/\beta_m)^2-1} \ \ when \ \ \beta_m \leq c \\
% \gamma_{\beta_m} = -i\sqrt{1-(c/\beta_m)^2} \ \ when \ \ \beta_m >c 
% \end{gathered}
% \label{eq:angle_formulation}
% \end{equation}

\begin{equation}
\begin{gathered}
\gamma_{\alpha_m} = \left\{ 
\begin{array}{ll}
\sqrt{(c/\alpha_m)^2-1} & \alpha_m \leq c \\
-i\sqrt{1-(c/\alpha_m)^2} & \alpha_m > c \\
\end{array}
\right.
\\
\gamma_{\beta_m} = \left\{ 
\begin{array}{ll}
\sqrt{(c/\beta_m)^2-1} & \beta_m \leq c  \\
-i\sqrt{1-(c/\beta_m)^2} & \beta_m >c  \\
\end{array}
\right.
\end{gathered}
\label{eq:angle_formulation}
\end{equation}

% 
Note  that  when compressional wave velocity $\alpha_{m}$ and/or shear wave velocity
$\beta_{m}$ are greater than the horizontal phase velocity $c$, the incidence
from $P$ or $SV$ wave is beyond the critical angle.
%
% {\Large \bf NOTE: what will this mean physically, will they be bound/constrained
% to layers?}
%
%%%%%==== Reply from Hexiang ===== 
% Professor, you are right. In the case of post critical, i.e., incidence angle beyond critical angle, the wave will be constrained to the layer. Some reflection  and refraction would not happen, for example, incident SV wave beyond crtical angle would not have normally reflected P wave, refracted P and SV wave. 
%
%
%
%
% {\Large \bf NOTE: Hexiang please make sure that the following statements are true!}
%
%
%%%%%==== Reply from Hexiang ===== 
%I have double checked this with some literature \citet{Lee1989} and book \citet{Aki2002}
%Actually this is how surface wave are generate. For example, Rayleigh wave is nothing but critial and post-critical P and SV waves travel together simultaneously along the free surface; Love wave comes from Post-Critical incidence of SH wave. For those surface waves, their magnitudes exponentially decay along the depth and therefore can not be observed in deeper ground. 
%
%
%
%
In that case, the incident and
reflected angles for $P$ and $SV$ wave, $\gamma_{\alpha_m}$  and 
$\gamma_{\beta_m}$ become  complex numbers. The plane
wave magnitude exponentially decays along the depth.
%
To   be   consistent  with  the  original  formulation  by~\citet{Haskell1953},
dilatational   wave   solution   $\Delta_{m}$  and  rotational  wave  solution
$\omega_{m}$         are         first         introduced  through 
Equation~(\ref{eq:dilatational_rotational_amplitude}).
%

\begin{equation}
\begin{gathered}
\Delta = \frac{\partial u_x}{\partial x} + \frac{\partial u_z}{\partial z}\\
\omega = \frac{1}{2}(\frac{\partial u_x}{\partial z}-\frac{\partial u_z}{\partial x})
\end{gathered}
\label{eq:dilatational_rotational_amplitude}
\end{equation}
%

% {\Large \bf NOTE: Hexiang, every variable in text has to be used with its name
% so you need to add names for each variable in text, for example
% in next sentence you have to name $\phi_{m}$ and $\Psi_{m}$, for example you will say:
%
% . . . 
% $P$ wave scalar potential for $m$-th layer $\phi_{m}$ and $S$ wave vector
% potential for $m$-th layer $\boldsymbol{\Psi}_{m}$ 
% are related to dilatational   wave   solution   $\Delta_{m}$  and  rotational 
% wave solution  $\omega_{m}$   through:  
% }
%
%%%%%==== Reply from Hexiang ===== 
% Got it. Variable name has been added before the variable symbol. 
%
%
%
%
%
%
%
where $P$ wave potential magnitude $\phi_{m}$ and $SV$ wave potential magnitude $\Psi_{m}$ of $m$-th layer are related to dilatational wave solution $\Delta_{m}$ and rotational wave solution $\omega_{m}$ through:

\begin{equation}
\begin{gathered}
\phi_{m} = -(\frac{\alpha_m}{w})^2\Delta_{m}\\
\Psi_{m} = 2(\frac{\beta_m}{w})^2\omega_{m}
\label{eq:relationship_potential_dilatation_rotation}
\end{gathered}
\end{equation} 



The   displacements   $(u_x,  u_y)$  and  interfacial  stresses  $(\sigma_{zz},
\tau_{zx})$ can be expressed using wave potential magnitudes $\phi$ and $\Psi$,
through Equations~(\ref{eq_Helmholtz_decomposition})
-
(\ref{eq:relationship_potential_dilatation_rotation}).
Similarly, the
displacement  and  stress  field  of  $m^{th}$  layer can be calculated from the
dilatational wave and rotational wave solutions $\Delta_{m}$ and $\omega_{m}$ as
%
\begin{flalign}
\begin{gathered}
u_x 
= 
\{-ik(\frac{\alpha_m}{\omega})^2[(\Delta_m^{'}
+
\Delta_m^{''})cos(k\gamma_{\alpha_m}z)
-
i(\Delta_m^{'}-\Delta_m^{''})sin(k\gamma_{\alpha_m}z)]\\
+
2ik\gamma_{\beta_m}(\frac{\beta_m}{\omega})^2[(\omega_m^{'}
-
\omega_m^{''})cos(k\gamma_{\beta_m}z)+i(\omega_m^{''} 
+
\omega_m^{'})sin(k\gamma_{\beta_m}z)]\}e^{ikx} \\  
\end{gathered}
\label{eq:stress_displacement_relation_begin}
\end{flalign}

\begin{flalign}
\begin{gathered}
u_z 
= 
\{ik\gamma_{\alpha_m}(\frac{\alpha_m}{\omega})^2[(\Delta^{'}_m-\Delta^{''}_m)
cos(k\gamma_{\alpha_m}z)-i(\Delta^{''}_m+\Delta^{'}_m)sin(k\gamma_{\alpha_m}z)]\\
+2ik(\frac{\beta_m}{\omega})^2[(\omega_m^{'}+\omega_m^{''})cos(k\gamma_{\beta_m}z)
-
i(\omega_m^{'}-\omega_{m}^{'})sin(k\gamma_{\beta_m}z)]\}e^{ikx}\\
\end{gathered}
\end{flalign}

\begin{flalign}
\begin{gathered}
\sigma_{zz} 
=
\rho_{m} \{\alpha_{m}^2(1-2\frac{\beta_{m}^2}{c^2})[(\Delta_m^{'}
+
\Delta_m^{''})cos(k\gamma_{\alpha_m}z)-i(\Delta_m^{'}-\Delta_m^{''})sin(k\gamma_{\alpha_m}z)]\\
+
4\frac{\beta_{m}^4}{c^2}\gamma_{\beta_m}[(\omega_m^{'}-\omega_m^{''})cos(k\gamma_{\beta_m}z)
-
i(\omega_m^{'}+\omega_m^{''})sin(k\gamma_{\beta_m}z)]\} e^{ikx}\\
\raisetag{2.5\normalbaselineskip}
\end{gathered}
\end{flalign}

\begin{flalign}
\begin{gathered}
\tau_{zx} 
= 
2\rho_m\beta_{m}^2 \{-\gamma_{\alpha_m}(\frac{\alpha_m}{c})^2[(\Delta_m^{'}
-
\Delta_m^{''})cos(k\gamma_{\alpha_m}z)-i(\Delta_m^{''}+\Delta_m^{'})sin(k\gamma_{\alpha_m}z)]\\
+[1-2(\frac{\beta_m}{c})^2][(\omega_m^{'}+\omega_m^{''})cos(k\gamma_{\beta_m}z)
-
i(\omega_m^{'}-\omega_m^{''})sin(k\gamma_{\beta_m}z)]\}e^{ikx}
\raisetag{1.2\normalbaselineskip}
\end{gathered}
\label{eq:stress_displacement_relation_end}
\end{flalign}


% {\Large \bf NOTE: paragraph below has to be rewritten, and all the equations have
% to be in a separate lines as regular equations with equation numbers...
%
% for example, did you mean it to be like this:
%
%
% If we use an intermediate variable  $S^{(m)}$, defined as
%
% \begin{equation}
% S^{(m)} = 
% [\dot    u_x(z_m=d_m)/c,   \dot   u_z(z_m=d_m)/c,
% \sigma_{zz}(z_m=d_m), \tau_{zx}(z_m=d_m)]^{T}
% \end{equation}
%
%
% something like this, as this paragraph below is not very clear!
%
% }
%
%
%%%%%==== Reply from Hexiang ===== 
% The paragraph below has been rewritten with equations separated out in separate lines 
%
%
% Define  the displacement and stress solution at $m^{th}$ interface as $S^{(m)}$,
% which    is    equal    to    $[\dot    u_x(z_m=d_m)/c,   \dot   u_z(z_m=d_m)/c,
% \sigma_{zz}(z_m=d_m), \tau_{zx}(z_m=d_m)]^{T}$.
%
%
%
%
%
Define the displacement and stress solutions at $m^{th}$ interface as column vector $S^{(m)}$: 

\begin{equation}
S^{(m)} = 
[\dot    u_x(z_m=d_m)/c,   \dot   u_z(z_m=d_m)/c,
\sigma_{zz}(z_m=d_m), \tau_{zx}(z_m=d_m)]^{T}
\end{equation}

Then Equations~(\ref{eq:stress_displacement_relation_begin})
-
(\ref{eq:stress_displacement_relation_end})  can  be  reduced  to  the following
matrix notations~\citep{Haskell1953}:

\begin{flalign}
S^{(m-1)} 
= 
\boldsymbol{E_m} [\Delta_{m}^{''}+\Delta_{m}^{'}, \Delta_m^{''}
-
\Delta_{m}^{'}, \omega_m^{''}-\omega_m^{'}, \omega_m^{''}+\omega_m^{'}]^{T}
\label{eq:solution_from_potential_magnitude}
\end{flalign}
\begin{flalign}
S^{(m)} 
= 
\boldsymbol{D_m} [\Delta_{m}^{''}+\Delta_{m}^{'}, 
\Delta_m^{''}-\Delta_{m}^{'}, 
\omega_m^{''}-\omega_m^{'}, 
\omega_m^{''}+\omega_m^{'}]^{T}
\end{flalign}

\noindent  where transformation matrix $\boldsymbol{E_m}$ and $\boldsymbol{D_m}$
are given as: 

% {\Large \bf NOTE: moved them from  appendix here, it might be better if they are
% locally presented}
%
%
%%%%%==== Reply from Hexiang ===== 
% OK. Moved equations from Appendix to body text. 
%
%


% \resizebox{\linewidth}{!}{
% $\boldsymbol{E_m} =
% \begin{bmatrix}
% -(\alpha_m/c)^2 & 0 & -\theta_m \gamma_{\beta_m} & 0 \\
% 0 & -(\alpha_m/c)^2\gamma_{\alpha_m} & 0 & \gamma_m \\
% -\rho_m\alpha_m^{2}(\theta_m-1) & 0 & -\rho_mc^2\theta_m^2\gamma_{\beta_m} & 0 \\
% 0 & \rho_m\alpha_m^2\theta_m\gamma_{\alpha_m} & 0 & -\rho_mc^2\theta_m(\theta_m-1)
% \end{bmatrix}$}
% \\

\begin{equation}
\resizebox{0.89\linewidth}{!}{
$\boldsymbol{E_m} =
\begin{bmatrix}
-(\alpha_m/c)^2 & 0 & -\theta_m \gamma_{\beta_m} & 0 \\
0 & -(\alpha_m/c)^2\gamma_{\alpha_m} & 0 & \gamma_m \\
-\rho_m\alpha_m^{2}(\theta_m-1) & 0 & -\rho_mc^2\theta_m^2\gamma_{\beta_m} & 0 \\
0 & \rho_m\alpha_m^2\theta_m\gamma_{\alpha_m} & 0 & -\rho_mc^2\theta_m(\theta_m-1)
\end{bmatrix}$}
\label{eq:matrix_Em}
\end{equation}
%
with $\theta_{m} = 2(\beta_m/c)^2$. 

\begin{equation}
\resizebox{0.89\linewidth}{!}{
$\boldsymbol{D_m} =
\begin{bmatrix}
-(\alpha_m/c)^2cosA_m & i(\alpha_m/c)^2sinA_m & -\theta_m \gamma_{\beta_m}cosB_m & i\theta_m\gamma_{\beta_m}sinB_m \\
i(\alpha_m/c)^2\gamma_{\alpha_m}sinA_m & -(\alpha_m/c)^2\gamma_{\alpha_m}cosA_m & -i\theta_msinB_m & \theta_mcosB_m \\
-\rho_m\alpha_m^{2}(\theta_m-1)cosA_m & i\rho_m\alpha_m^2(\gamma_m-1)sinA_m & -\rho_mc^2\theta_m^2\gamma_{\beta_m}cosB_m & i\rho_mc^2\theta_m^2\gamma_{\beta_m}sinB_m \\
-i\rho_m\alpha_m^{2}\theta_m\gamma_{\alpha_m}sinA_m & \rho_m\alpha_m^2\theta_m\gamma_{\alpha_m}cosA_m & i\rho_mc^2\theta_m(\theta_m-1)sinB_m & -\rho_mc^2\theta_m(\theta_m-1)cosB_m
\end{bmatrix}$}
\label{eq:matrix_Dm}
\end{equation}
%
%
with $A_m=k\gamma_{\alpha_m}d_m$ and $B_m=k\gamma_{\beta_m}d_m$. 

%
The   recurrence   relation  between  $S^{(m)}$  and  $S^{(m-1)}$  then  can  be
established  as  shown in Equation~\ref{eq:recursive_relation},  where it was used that
$\boldsymbol{G_{m}} = \boldsymbol{D_m}\boldsymbol{E_m^{-1}}$.
%

\begin{equation}
S^{(m)} 
= 
\boldsymbol{D_m}\boldsymbol{E_m^{-1}}S^{(m-1)} 
= 
\boldsymbol{G_{m}}S^{(m-1)}
\label{eq:recursive_relation}
\end{equation}



Recursively    applying    Equation~\ref{eq:recursive_relation}   leads   to  
Equation~\ref{eq:upper_layer_to_bottom_layer}.
%
Using   the  relation between displacement, stress response at $(m-1)^{th}$ interface $S^{(m-1)}$  and dilatational, rotational wave solutions $\Delta_m,\  \omega_m$,  Eq.
\ref{eq:boundary_bridge}  bridges  the  gap  between  the  upper boundary (i.e.,
response  at  ground  surface  $S^{(0)}$) and lower boundary (i.e., solutions of
wave  incident  layer  $\Delta_n$ and $\omega_n$), upon which specific boundary
conditions can be imposed.
%
\begin{equation}
S^{(n-1)} = \prod_{i=1}^{n-1}\boldsymbol{G_{i}}S^{(0)} 
\label{eq:upper_layer_to_bottom_layer}
\end{equation}
%

\begin{equation}
\begin{gathered}
S^{(0)} 
= 
\boldsymbol{L}[\Delta_n^{''}+\Delta_n^{'}, \Delta_n^{''}
-
\Delta_n^{'}, \omega_n^{''}-\omega_n^{'}, \omega_n^{''}+\omega_n^{'}]^{T} \\
%
\boldsymbol{L} 
= 
({\prod_{i=1}^{n-1}\boldsymbol{G_{i}}})^{-1}\boldsymbol{E_n}
\label{eq:boundary_bridge}
\end{gathered}
\end{equation}

The following  boundary  conditions  are  incorporated:  

\begin{enumerate}

  \item  At $n^{th}$ layer, the incident  in-plane  $P$  and  $SV$  wave 
  potential magnitude $\phi_n^{'}$ and $\Psi_n^{'}$ are given as $K_1$ and
  $K_2$;
%

  \item  At  the  ground  surface ($z=0$), the traction is free, i.e., the third
  and fourth component of surface response vector $S^{(0)}$ are 0.
%

\end{enumerate}


Therefore,  the  reflected  dilatational  wave  magnitude  and  rotational  wave
magnitude  can be solved using Equation~\ref{eq:n_layer_magnitude}, where
$\Delta_n^{'}$ is $-K_1\omega^2/\alpha_n^2$ and $\omega_n^{'}$ is
$K_2w^2/(2\beta_n^2)$.

\vspace{-0.5cm}
\begin{equation}
\resizebox{0.91\hsize}{!}{
$\begin{bmatrix}
\Delta_n^{''}\\\omega_n^{''}
\end{bmatrix}
=
\begin{bmatrix}
L_{31} + L_{32} & L_{33} + L_{34} \\
L_{41} + L_{42} & L_{43} + L_{44}
\end{bmatrix}^{-1}
\begin{bmatrix}
(L_{32}-L_{31})\Delta_{n}^{'} + (L_{33}-L_{34})\omega_{n}^{'}\\
(L_{42}-L_{41})\Delta_{n}^{'} + (L_{43}-L_{44})\omega_{n}^{'}
\end{bmatrix}$}
\label{eq:n_layer_magnitude}
\end{equation}

\vspace{0.2cm}
%
Finally,  recurrence  relation, given by   Equation~\ref{eq:trace_back_relation}

\begin{equation}
\begin{bmatrix}
\Delta_{m-1}^{''}+\Delta_{m-1}^{'}\\
\Delta_{m-1}^{''}-\Delta_{m-1}^{'}\\
\omega_{m-1}^{''}-\omega_{m-1}^{'}\\
\omega_{m-1}^{''}+\omega_{m-1}^{'}
\end{bmatrix}
=
\boldsymbol{D_{m-1}^{-1}E_m}
\begin{bmatrix}
\Delta_{m}^{''}+\Delta_{m}^{'}\\
\Delta_{m}^{''}-\Delta_{m}^{'}\\
\omega_{m}^{''}-\omega_{m}^{'}\\
\omega_{m}^{''}+\omega_{m}^{'}
\end{bmatrix}
\label{eq:trace_back_relation}
\end{equation}
%


\noindent
can be used to
trace back dilatational wave magnitude $\Delta_m$ and rotational wave magnitudes
$\omega_m$ for the rest $n-1$ layers.
%
Based  on solution for dilatational and rotational magnitudes for each layer, the
complete  displacement  and  stress  field  can  be  easily computed, using  
Equations~(\ref{eq:stress_displacement_relation_begin}) - 
(\ref{eq:stress_displacement_relation_end}).
%

In addition, viscosity can also be included with slight modification.
%
Considering  Kelvin-Voight viscoelastic material~\citep{Chirictua2008}, viscosity
can  be  handled  with  complex Lam{\'e} modulus and wave velocities as shown in Eq.
\ref{eq:viscosity modification}, where $\xi$ is the damping ratio.

\begin{equation}
G^{*} 
= 
G(1+2\xi i) \ \ \ \beta_m^{*} \simeq \beta_m (1+\xi i) \ \ \ \alpha_m^{*} \simeq \alpha_m(1+\xi i)
\label{eq:viscosity modification}
\end{equation}   










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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Domain Reduction Method} 
\label{The_Domain_Reduction_Method}

Domain    Reduction    Method    (DRM)    was    originally developed for
studying      local      topography      effects      on     seismic    
motions~\citep{Bielak2001,Yoshimura2001}, while earlier
work~\citep{Bielak1984,Cremonini1988} did note soil-structure  interaction 
modeling as the ultimate goal.
%
%It was then adopted to solve SSI with great success \cite{Jeremic2008d,Abell2016a,Solberg2016}. 
%
In  the context of DRM, engineering system is discretized using the finite
element method over interior domain $\Omega$, within boundary $\Gamma$,
containing local SSI system and reduced exterior domain $\Omega^{+}$, outside of
boundary $\Gamma$.
%
The nodes of the finite element model are then placed in three categories:
interior nodes, boundary nodes
between domains $\Omega$ and $\Omega^{+}$, on the boundary $\Gamma$, and
exterior nodes in exterior
domain $\Omega^{+}$.
%
%
Corresponding  nodal  displacements are denoted as $u_{i}$, $u_{b}$ and $u_{e}$,
for interior, boundary and exterior nodes, respectively.
%
Boundary  nodes  and  their  connected  exterior  nodes  form  a single layer of
elements, called DRM layer, surrounding the interior SSI domain.
%
The power  of DRM lies in the analytical formulation of effective seismic forces
$P^{eff}$, given by the Equation~\ref{eq:original_DRM}.
%
\begin{equation}
P^{eff} = 
\begin{Bmatrix}
P_i^{eff} \\
P_b^{eff} \\
P_e^{eff}
\end{Bmatrix}
= 
\begin{Bmatrix}
0 \\
-M_{be}^{\Omega^{+}}\ddot{u}^{0}_{e} - K_{be}^{\Omega^{+}} u_e^0 \\
M_{eb}^{\Omega^{+}} \ddot{u}^{0}_{b} + K_{eb}^{\Omega^{+}} u_b^0
\end{Bmatrix}
\label{eq:original_DRM}
\end{equation}
%



\noindent
Effective seismic forces $P^{eff}$ represent a dynamically consistent replacement 
for seismic  forces at the hypocenter. 
%
Effective seismic forces $P^{eff}$ are applied to the DRM layer, and produce the free field motions in a domain without local SSI system.
%
The effective forces are developed from free field seismic motions, hence for
free field finite element models, there are no seismic motions leaving the system.
%
When the structure is present, during SSI analysis the only outgoing motions are
related to the radiation damping of structural motions.




%
From  Eq. \ref{eq:original_DRM}, only free field motions ($u_e^{0}$, $u_b^0$) at
nodes of DRM layer and element mass and stiffness matrix ($M_{be}^{\Omega^{+}}$,
$K_{be}^{\Omega^{+}}$)  of  DRM layer are required to calculate effective forces
$P^{eff}$.
%
Free field motions developed in the previous section are used in creation of the effective seismic forces as per Equation~\ref{eq:original_DRM}.
%



Presented approach, using analytic solution for free field 3 component (3C) seismic motions, that feature both body and surface waves, is more efficient and straightforward than conventional substructure method.
%
In  addition  to  free field motions, substructure method requires to solve
foundation wave scattering and impedance function, both of which are challenging
tasks.
%
It is noted that very few  specific  shapes of foundation, e.g., circular and
rectangular shape, embedded in simplified ground conditions have been studied using sub-structuring
method~\citep{Luco1974,Wong1985,Crouse1990,Luco1976,Luco1982,De1995,Luco1987,Fu2017}.
%
For the presented approach, free  field  motions under inclined incident plane  waves are solved using wave potential formulation.
%
Both  wave  scattering  and  dynamic  SSI  are automatically handled by the time
domain FEM analysis that is dynamically loaded with DRM effective earthquake forces.
%
In addition, developed Wave Potential Formulation -- Domain Reduction Method (WPF-DRM) offers advantages for solving locally inhomogeneous and nonlinear SSI problems under inclined seismic excitations~\citep{Jeremic2004d,Wang2017,Sinha2017b}.
%









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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Illustrative Examples}
\label{sec:numerical_study}


Presented  WPF-DRM method is implemented in the Real-ESSI
Simulator~\citep{Real_ESSI_Simulator}. Described examples and publicly available
executables for the Real ESSI sequential and parallel programs are available
through Real ESSI Simulator web site \url{http://real-essi.info/}.
%
All  numerical examples presented here are analyzed using Real-ESSI Simulator
version 20.01, in parallel computing mode, on UC Davis and Amazon Web Services
parallel computers.
%









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Free Field Modeling and Verification}
\label{sec:numerical_study_verification}



Free  field  response of layered ground excited by an inclined incident seismic
wave is used to illustrate and verify developed methodology.
%
Analytic solutions based on  Thomson-Haskell propagation matrix
technique~\citep{Thomson1950,Haskell1953,Silva1976} are used for verification.
% %
% With  effective  forces  applied  on a DRM layer, the free field motions should be
% recovered  by  FEM  analysis  (elastic resumption behavior) and should match the
% analytical solution.
% %
% In  addition,  there  should  not  be any scattered wave outgoing from DRM
% layer~\citep{Bielak2001}.
% %

A finite  element  model for the free field, that is $300$m wide and $200$m
deep, consisting of three layers, as described in
Table~\ref{table:layering_property}, is used.
%
\begin{table}[!htbp]
\flushleft
\caption{Properties of layers:
thickness $d$, 
density $\rho$, 
shear wave velocity $V_s$,
compressional wave velocity $V_p$ and
Poisson's ratio $\nu$.}
\label{table:layering_property}
\begin{tabular}{cccccc}
\hline
%Layer ID & Thickness $d$ [m]     & density $\rho$ [$kg/m^3$]   & $V_s$ [$m/s$]   & $V_p$ [$m/s$]  & Poisson's ratio $\nu$ \\ \hline
Layer & $d$ [m]  & $\rho$ [$kg/m^3$] & $V_s$ [$m/s$] & $V_p$ [$m/s$] & $\nu$ \\ \hline
1     & 50       & 2100              & 500           & 816.5         & 0.2 \\
2     & 100      & 2300              & 750           & 1403.1        & 0.3 \\
3     & $\infty$ & 2500              & 1000          & 2081.7        & 0.35\\ \hline
\end{tabular}
\end{table}


It is noted that dimension of analyzed model is  $300$m $\times$ $200$m, however
there exist additional finite elements outside this domain, for the DRM layers,
as well as additional higher Rayleigh damping layers outside to damp out any
outgoing waves. It is also noted that theoretically there should be no waves
propagating outside of the DRM layer for a free field response. Additional
damping layers are added in order to accommodate further, non-free field model
expansions and additions.
%
Finite element size is set to $5$m, and with 10 finite elements per wave length,
this mesh can accurately propagate waves of up to $f=10$Hz, for surface soil
with shear wave velocity of $V_s=500 {\rm m/s^2}$, as per
\citet{Lysmer1969,Watanabe2015}.
%





A number of monochromatic, single frequency  plane  SV  wave, represented by a
cosine function, with variable inclinations $\theta=10^{o}, 45^{o}, 60^{o}, 80^{o}$ 
and variable frequencies, $f=1.0, 2.5, 5.0, 10.0$Hz,
are applied to the layered
ground model using developed methodology.
%
It is noted that inclination angle $\theta$ is measured between a wave
propagation direction vector and vertical axes.
%
The  incident SV wave magnitude from the depth is $0.06$m and is kept the same
for all the analyzed cases.
%
%
Free  field  motions  are developed and introduced into the model through WPF-DRM.
%
Figure~\ref{fig:free_field_different_inclination} shows snapshots of wave
displacements in the model, for a wave frequency of $f=5$Hz, for different
input plane wave inclinations, $\theta=10^{o}, 45^{o}, 60^{o}, 80^{o}$.
%
%
%
\begin{figure}[!htbp]
 \flushleft
 \subfloat[$f=5Hz\ \theta=10^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_10_degree_snippet.pdf}}
 \subfloat[$f=5Hz\ \theta=45^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_45_degree_snippet.pdf}} 
 \\
 \subfloat[$f=5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_60_degree_snippet.pdf}}
 \subfloat[$f=5Hz\ \theta=80^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_80_degree_snippet.pdf}}
 \\
 \caption{\label{fig:free_field_different_inclination}
 Displacement magnitudes for a free  field  response  under  incident SV wave,
 frequency $f=5$Hz, with different incident wave inclinations: (a)
 $\theta=10^{\circ}$ (b) $\theta=45^{\circ}$ (c) $\theta=60^{\circ}$ (d)
 incident angle $\theta=80^{\circ}$.}
\end{figure}










Figure~\ref{fig:free_field_different_frequency} shows snapshots of wave
displacements in the model, for a wave that is inclined at $\theta=60^{o}$, 
for variable input plane wave frequencies $f=1.0, 2.5, 5.0, 10.0$Hz.
%
%
%
\begin{figure}[!htbp]
 \flushleft
 \subfloat[$f=1Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.48\columnwidth]{Figures/free_field_1_frequency_snippet.pdf}}
 \subfloat[$f=2.5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_2_frequency_snippet.pdf}}
 \\
 \subfloat[$f=5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_5_frequency_snippet.pdf}}
 \subfloat[$f=10Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/free_field_10_frequency_snippet.pdf}}
 \\
 \caption{\label{fig:free_field_different_frequency}
Displacement magnitudes for a free  field  response  under  incident  SV wave at
an angle of $\theta=60^{o}$,
with different frequencies: (a) $f=1.0$Hz (b) $f=2.5$Hz (c) $f=5.0$Hz (d)
$f=10.0$Hz.}
\end{figure}




Few notes are in order upon visual inspection of results in
Figures~\ref{fig:free_field_different_inclination} and
\ref{fig:free_field_different_frequency}.
%



The  outgoing  waves in exterior region, outside DRM layer, are negligibly
small,
almost zero for all the cases. This is indeed expected, as it follows from the
theory of the domain reduction method \citep{Bielak2001,Yoshimura2001}, whereby
the so called residual field ($w_e$) should be non-existent for free field
motions, that were used to develop effective DRM forces.


%
Comparing  free field responses for SV wave with  different incident angles,
Figure \ref{fig:free_field_different_inclination}, the $\theta=10^{\circ}$ case
behaves very similar to 1D vertically propagating motion field that is commonly
used in engineering practice.
% 
It is noted, however that there are still vertical motions at the surface due to
such almost vertical SV wave interacting with the free surface.
%
For cases where wave inclination is more significant,  for   $\theta=45^{\circ}$  and
$\theta=60^{\circ}$, significant surface motions are observed, with pronounced
vertical and horizontal motions.
%
%  They  might  not  be simplified as 1D wave field any more because of the spatial
%  variation.
%  %
When  the  incident wave inclination  is  $\theta=80^{\circ}$,  seismic wave
propagates almost horizontally without generating significant surface motions.
%
It is also noted that the  displacement  magnitude  of  the seismic wave field for 
wave inclination case $\theta=80^{\circ}$ is much smaller than for the other cases.
%
This  is reasonable considering the site amplification for other free field cases
comes, in part, from the impedance contrast of vertical wave propagation.
%






Results, snapshots of displacement field  magnitudes for wave fields of different       
frequencies are shown in Figure~\ref{fig:free_field_different_frequency}  for
seismic motion inclined SV wave field at $\theta=60^{\circ}$.
%
% The horizontal wave length for $1Hz$ is $1155$m.
% %
It is noted that layer boundaries, impedance contrasts, are at $-50$m, and at
$-150$m.
%
Those layer boundaries can be visually identified from distribution of waves
through model depth with positive and negative interference reflected and
refracted waves within different layers of the domain.




% 
% Therefore,   the  free  field  response  is  approximately  uniform  within  the
% relatively small simulation region ($300$m wide).
% %
% While  for  the  cases  of  $5Hz$  and  $10Hz$,  significant surface motions are
% observed with much smaller wave length.
% %
% Variations  of  site  amplification  among  four  cases  of  different  incident
% frequencies are also evident.
% %
% This  is  because  the  site  amplification  depends  on  the difference between
% frequency  of  incident  excitation  and  fundamental  frequency  of the layered
% system.
% %

% \begin{figure}[!htbp]
%  \flushleft
%  \subfloat[$f=5Hz\ \theta=10^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/10degree_magnitude_verification.pdf}}
%  \subfloat[$f=5Hz\ \theta=45^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/45degree_magnitude_verification.pdf}}
%  \\
%  \subfloat[$f=5Hz\ \theta=60^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/60degree_magnitude_verification.pdf}}
%  \subfloat[$f=5Hz\ \theta=80^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/80degree_magnitude_verification.pdf}}
%  \\
%  \caption{\label{fig:free_field_verification_different_inclination.}
% Verification of free field modeling under incident SV wave with different inclinations: (a) incident angle $10^{\circ}$  (b) incident angle $45^{\circ}$ (c) incident angle $60^{\circ}$ (d) incident angle $80^{\circ}$}
% \end{figure}

% \begin{figure}[!htbp]
%  \flushleft
%  \subfloat[$f=1Hz\ \theta=60^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/1Hz_magnitude_verification.pdf}}
%  \subfloat[$f=2.5Hz\ \theta=60^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/2Hz_magnitude_verification.pdf}}
%  \\
%  \subfloat[$f=5Hz\ \theta=60^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/5Hz_magnitude_verification.pdf}}
%  \subfloat[$f=10Hz\ \theta=60^{o}$]{
%  \includegraphics[width=0.55\columnwidth]{Figures/10Hz_magnitude_verification.pdf}}
%  \\
%  \caption{\label{fig:free_field_verification_different_frequency.}
% Verification of free field modeling under incident SV wave with different frequencies: (a) frequency $1Hz$  (b) frequency $2.5Hz$ (c) frequency $5Hz$ (d) frequency $10Hz$}
% \end{figure}



Figures~\ref{fig:free_field_verification_different_inclination}       and
\ref{fig:free_field_verification_different_frequency} compare simulated free
field   horizontal  and  vertical  displacement  magnitudes against corresponding
analytical solutions along the depth. 
%
%
\begin{figure}[!htbp]
 \flushleft
 \subfloat[$f=5Hz\ \theta=10^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_10degree_new.pdf}}
 \subfloat[$f=5Hz\ \theta=45^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_45degree_new.pdf}}
 \\
 \subfloat[$f=5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_5Hz_new.pdf}}
 \subfloat[$f=5Hz\ \theta=80^{o}$]{
 \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_80degree_new.pdf}}
 \\
 \caption{\label{fig:free_field_verification_different_inclination}
Verification  of  free  field  modeling  under  incident  SV wave with different
incident angles $\theta$:  
(a) $\theta= 10^{\circ}$, 
(b) $\theta= 45^{\circ}$,
(c) $\theta= 60^{\circ}$ 
(d) $\theta= 80^{\circ}$.}
\end{figure}
%
It is noted that acceleration magnitudes can be obtained by
multiplying displacement magnitudes with $w^2$.
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Free field interference characteristic in terms of acceleration %%%%%%%% 
%%%% We found it is better to show interference in terms of acceleration %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% \begin{figure}[!htbp]
% \flushleft
% \subfloat[$f=5Hz\ \theta=10^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_10degree_new_Acc.pdf}}
% \subfloat[$f=5Hz\ \theta=45^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_45degree_new_Acc.pdf}}
% \\
% \subfloat[$f=5Hz\ \theta=60^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_5Hz_new_Acc.pdf}}
% \subfloat[$f=5Hz\ \theta=80^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_80degree_new_Acc.pdf}}
% \\
% \caption{\label{fig:free_field_verification_different_inclination}
% Verification  of  free  field  modeling  under  incident  SV wave with different
% inclinations:  (a)  incident  angle $10^{\circ}$ (b) incident angle $45^{\circ}$
% (c) incident angle $60^{\circ}$ (d) incident angle $80^{\circ}$.}
% \end{figure}

%
\begin{figure}[!htbp]
\flushleft
\subfloat[$f=1Hz\ \theta=60^{o}$]{
\includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_1Hz_new.pdf}}
\subfloat[$f=2.5Hz\ \theta=60^{o}$]{
\includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_2Hz_new.pdf}}
\\
\subfloat[$f=5Hz\ \theta=60^{o}$]{
\includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_5Hz_new.pdf}}
\subfloat[$f=10Hz\ \theta=60^{o}$]{
\includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_10Hz_new.pdf}}
\\
\caption{\label{fig:free_field_verification_different_frequency}
Verification  of  free  field  modeling  under  incident  SV wave with different
frequencies $f$:  
(a) $ f = 1.0$Hz, 
(b) $ f = 2.5$Hz, 
(c) $ f = 5.0$Hz, 
(d) $ f = 10.0$Hz.}
\end{figure}
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Free field interference characteristic in terms of acceleration %%%%%%%% 
%%%% We found it is better to show interference in terms of acceleration %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \begin{figure}[!htbp]
% \flushleft
% \subfloat[$f=1Hz\ \theta=60^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_1Hz_new_Acc.pdf}}
% \subfloat[$f=2.5Hz\ \theta=60^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_2Hz_new_Acc.pdf}}
% \\
% \subfloat[$f=5Hz\ \theta=60^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_5Hz_new_Acc_frequency.pdf}}
% \subfloat[$f=10Hz\ \theta=60^{o}$]{
% \includegraphics[width=0.50\columnwidth]{Figures/3C_depth_variation_10Hz_new_Acc.pdf}}
% \\
% \caption{\label{fig:free_field_verification_different_frequency}
% Verification  of  free  field  modeling  under  incident  SV wave with
% different
% frequencies:  (a)  frequency $1Hz$ (b) frequency $2.5Hz$ (c) frequency $5Hz$ (d)
% frequency $10Hz$.}
% \end{figure}
% 




%
Very  good  agreement is observed between results given by WPF-DRM simulation and
analytical solutions.
%  %
%  This  verifies  and  confirms the capability of WPF-DRM to recover the free field
%  response with known motions only at the DRM layer.
%  %
%
%
%
Several interesting observations can also be made:

\begin{enumerate}


 \item Along with the increase
in frequencies, the vertical wave length becomes shorter, and that results in more
wave interferences along the depth.
%

 \item  The  existence  of  layers and interfaces at $z=-50$m and $z=-150$m
complicates  the spatial variation of wave field along the depth, especially for
higher frequencies, $f=5Hz$ and $10Hz$.
%
The response   curves  at  depths $0\sim50$m  and $50\sim150$m
are quite different in both amplitude and variation pattern.
%

 \item  From Fig. \ref{fig:free_field_verification_different_inclination}, it can be
seen  that  inclination angle   of  input  SV  wave  also  plays  a crucial role in the
interference characteristic of inclined wave field.
%
Periodic  peaks  and  troughs  shown in the case of $10^{\circ}$ inclination are
typical  interference  characteristics of 1D homogeneous, vertically propagating
wave field.
% 
However,  the  interference  characteristics  given  by  other wave inclinations
show significant differences.
%
These  different  variation patterns along the depth, that  might not make much
difference for shallow founded surface structures, can result in very different
seismic response for deeply embedded structures.
%



\end{enumerate}






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Deeply Embedded Soil-Structure Model}
\label{Deeply_Embedded_Soil-Structure_Model}


Deeply  embedded  structural  model, a model of a Small Modular Reactor (SMR) is
analyzed and used to illustrate developed methodology.
%
%
The  FEM model  of an SMR  structure  embedded in layered
ground is shown in Figure~\ref{fig_full_model_SMR}(a).
%
%
\begin{figure}[!htbp]
\centering
\subfloat[FEM model of SSI system with embedded SMR]{
\includegraphics[width=0.50\textwidth]{Figures/full_model_SMR.pdf}\enspace
\includegraphics[width=0.50\textwidth]{Figures/half_model_SMR_annotation_crop.pdf}}\\
\subfloat[Representative points configuration]{
\includegraphics[width=0.50\textwidth]{Figures/Points_configuration.pdf}}
\caption{\label{fig_full_model_SMR} FEM model of embedded SMR and representative points.}
\end{figure}
%
%
The  embedment   depth is $36$m, while the height of SMR structure above ground is
$14$m.
%
Eleven    representative    points,    point    A    to    point   K   in
Figure~\ref{fig_full_model_SMR}(b), are selected to monitor the dynamic response
of SMR.
%
%More detailed description of geometry of this model SMR  is  documented
%by~\citet{Solberg2016}.
%
The  layered  ground  parameters  are the same as those used in free field study
given in Table \ref{table:layering_property}.
%



To  proper  model wave propagation, the finite element size and time step should
be carefully controlled to reduce discretization errors.
%
For  linear displacement approximation within finite element, in this case
eight-node brick elements, at least 10
nodes per wavelength should be used \citep{Watanabe2015}.
%
The time  step  length $\Delta  t$ is limited by Courant-Friedrichs-Lewy
condition~\citep{Courant1989} for stability.
%
In  addition,  following  requirement  needs to be met to accurately capture the
propagation of wave front~\citep{Jeremic2008d}, where $\Delta h$ is the mesh
size and $v$ is the highest wave velocity.
%
\begin{equation}
\Delta t < \frac{\Delta h}{v}
\end{equation}


In  this study, eight-node brick element with $4$m mesh size is used for spatial
discretization.
%
The  maximum frequency the model can propagate is about $12.5Hz$ considering the
minimum elastic shear wave velocity $500m/s$.
%
Time step  is chosen as $\Delta t = 0.005$s.
%
Newmark  time  integration  method~\citet{Newmark1959}  is  used, with small
amount of numerical, algorithmic damping that  is used  to damp out unrealistic
high frequency responses
introduced by
spatial discretization \citet{local-87}.
%
Gradually  increasing  Rayleigh  damping (7\%, 15\% and 30\%) is assigned to the
inner, middle and exterior part of the absorbing layers, outside of the DRM
layer, to prevent reflection of radiated outgoing
waves~\citep{Jeremic2008d,Abell2016a}.




% %
% For  sensitivity study, the same monochromatic SV waves at different frequencies
% and  inclinations  as  those  used  in  the free field study are incident to the
% layered SSI system of SMR.
% %
% For  each  numerical  case,  two  loading  stages are simulated: Self weight and
% dynamic loading with effective DRM forces.
% %
% The  simulation  results  and analyses of dynamic response of SMR are elaborated
% below.








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{SMR Excited with  Inclined SV Waves}
\label{SMR_Excited_with_Inclined_SV_Waves}



Deeply  embedded  SMR  structure  is  excited  with  inclined  plane  waves,  at
inclination angles of
$\theta  =  10^{\circ},  45^{\circ},  60^{\circ}$ and $80^{\circ}$. Seismic wave
frequency used for this set of numerical test was set at $f=5Hz$. As described in
table~\ref{table:layering_property}  on  page~\pageref{table:layering_property},
shear  wave velocities of top $50$m layer is $V_s=500$m/s while the lower layer
is  $100$m  think and has a shear wave velocoty of $V_s=750$m/s. Due to presence
of  layers,  seismic  wave  field  close  to the surface is made up Rayleigh and
Stoneley  waves~\citep{Aki2002,Semblat2009}. It might thus be difficult to
separate  influence  of  these  different surface waves the response of the SMR.
For       example,   in   Figure~\ref{fig:free_field_different_inclination}     
on  page~\pageref{fig:free_field_different_inclination}, that shows displacement
magnitudes  at  certain  time, for different inclination of incident plane wave,
Stoneley  wave  is apparent close to depth of $50$m. In addition,  Rayleigh wave
is also apparent close to free field surface. 
%
Those  wave  fields,  when  applied  to the SMR SSI system, produce response, at
location  of  point  A\footnote{Location  of  point  A  is  in the middle of SMR
structure,   where  center  of  the  free  field  model  would  be,  please  see
Figure~\ref{fig_full_model_SMR}  on  page~\pageref{fig_full_model_SMR}.}  on  SMR
structure,
%
%excited by an incident SV wave with
%frequency $f=5Hz$ and variable inclination angles, $\theta = 10^{\circ},
%45^{\circ}, 60^{\circ}$ and $80^{\circ}$ are shown in
as        shown        in       Figures~\ref{fig_point_A_inclinations}       and
\ref{fig_point_A_inclinations_acceleration}. 
%
%
%
\begin{figure}[!htb]
\centering
\vspace*{-2cm}
\subfloat[Horizontal displacement]{
\includegraphics[width=0.9\textwidth]{Figures/inclination_ux.pdf}}\\
\subfloat[Vertical displacement]{
\includegraphics[width=0.9\textwidth]{Figures/inclination_uz.pdf}}
\caption{\label{fig_point_A_inclinations}
Displacement  response  of  point A within  embedded SMR, excited by an inclined SV
wave with $f=5Hz$ and
different inclination angles, $\theta = 10^{\circ}, 45^{\circ}, 60^{\circ}$ and
$80^{\circ}$: (a) horizontal displacement (b) vertical displacement.}
\end{figure}
%
\begin{figure}[!htb]
\centering
\vspace*{-2cm}
\subfloat[Horizontal acceleration]{
\includegraphics[width=0.9\textwidth]{Figures/inclination_Ax.pdf}}\\
\subfloat[Vertical acceleration]{
\includegraphics[width=0.9\textwidth]{Figures/inclination_Az.pdf}}
\caption{\label{fig_point_A_inclinations_acceleration}
Acceleration  response of point A within  embedded SMR, excited by an  inclined
SV wave with $f=5Hz$ and
different inclination angles, $\theta = 10^{\circ}, 45^{\circ}, 60^{\circ}$ and
$80^{\circ}$: (a) horizontal acceleration (b) vertical acceleration.}
\end{figure}
%


It is noted that corresponding  free  field  motions  at  the  same location are
also plotted for comparison.
%
Variations  of  displacement  magnitudes  caused  by  different  inclinations of
incident   SV   wave   are  quite  noticeable  for  vertical  displacements  and
accelerations, while influence on horizontal displacements and accelerations is less severe.
%
% %
% Among  these four cases, SSI effects are not significant for horizontal response
% in  the sense that structural horizontal displacements are almost identical with
% corresponding free field motions at the same location.
% %
% In  contrast,  striking  differences  can be observed for vertical displacements
% between SMR and free field model.
% %
The    reduction of vertical displacement and accelerations that is  observed in
all  the  four  cases,  is  consistent  with  the concept of ``base averaging'',
``ironing out'' of seismic motions by~\citet{Housner1957}.
%
The    most  significant  reduction occurs for the  case  of incident wave at an
angle  $\theta=45^{\circ}$  while  little  reduction  is  seen  in  the  case of
$\theta=80^{\circ}$.
%

%%%% === Commented this discussions about V/H by hexiang ======

% Furthermore, the ratio of vertical displacement to horizontal displacement $Disp_V/Disp_H$
% significantly varies with the inclination of SV wave.
% %
% The  ratio $Disp_V/Disp_H$  is greater than $1.0$ in the case of SV wave inclined at
% $\theta=45^{\circ}$ in
% contrast to a
% very small ratio $Disp_V/Disp_H < 1/10$ for SV wave inclined at $\theta=80^{\circ}$.
% %

%OVDE
%OVDE

% {\large \bf  V/H is for accelerations, and  not for displacements, in fact it is for 
% response spectra, as far as I can tell form those papers by Bozorgnia and others!

% So we should get more results in terms of accelerations! 

% for example figure 4, 5, 9

% Also results for figure 11 12, 13 and 14 are for $\theta = 60^{\ast}$ however
% reviewers will question this as they mostly believe that waves are very close to
% vertical, so let us follow rule of thumb that Norm gave us that he heard from
% Walter Silva, waves should be $\theta = 15^{\ast}$ from vertical. 


% So it is good to
% show all these other examples, other angles, but let us use more results for
% $\theta = 15^{\ast}$ for V/H...

% We should not focus to much on V/H here, we'll get another paper on that...

% }

%%%%===================================================
%%%%===== Hexiang's reply =============================
%%%%===================================================

%%% 
%%% Yes, V/H is defined based on response spectra in engineering pratices. 
%%% As you suggetsed, we might get another paper to investigate V/H by running
%%% the cases of $\theta = 15^{\ast}$.  
%%% Have removed these comments here for V/H. 
%%%
%%%
%%% For the acceleration representation, figure 4 and 5 are interference characteristic 
%%% curves that are physically defined in terms of ground movement (i.e., displacement),
%%% so it is better to represent with displacement. Added note in the text said the 
%%% acceleration interference curve can be obtained by scaling displacement interference curve
%%% by a factor of w^2, where w is the angular frequency. 
%%%

% Therefore,  constant  $V/H$  value  $1/2  \sim  2/3$  that  is  commonly  adopted in
% engineering  practices might be inappropriate~\citep{Bozorgnia2004b}, especially
% in  near  source  regions where the seismic motions feature the relatively large
% inclinations  and  for  those structures that are susceptible to vertical ground
% motions~\citep{Papazoglou1996}.
% 



% Propose to remove it---  {\Large Hexiang:} this figure below is not quite intuitive, there is too much info and it is hard to 
% Propose to remove it--- make  distinction  between points! I think what we can remove it as it is really
% Propose to remove it--- complicated to gain undestanding of what is going on...
% Propose to remove it--- 
% Propose to remove it--- 
% Propose to remove it--- The  spatial  variation  of  displacement  magnitudes for all the representative
% Propose to remove it--- points                      are                     shown                     in
% Propose to remove it--- Figure~\ref{fig_spatial_distribution_magnitude_SMR_inclination}.
% Propose to remove it--- %
% Propose to remove it--- %
% Propose to remove it--- %
% Propose to remove it--- \begin{figure}[!htbp]
% Propose to remove it--- \centering
% Propose to remove it--- \subfloat[Horizontal displacement]{
% Propose to remove it--- \includegraphics[width=0.45\textwidth]{Figures/spatial_mag_distribution_Ux_inclination_new.pdf}}
% Propose to remove it--- \subfloat[Vertical displacement]{
% Propose to remove it--- \includegraphics[width=0.63\textwidth]{Figures/spatial_mag_distribution_Uz_inclination.pdf}}
% Propose to remove it--- \caption{\label{fig_spatial_distribution_magnitude_SMR_inclination}      Spatial
% Propose to remove it--- variation  of  vertical  displacement  magnitudes  of  representative points for
% Propose to remove it--- different inclinations.}
% Propose to remove it--- \end{figure}
% Propose to remove it--- %
% Propose to remove it--- %
% Propose to remove it--- %
% Propose to remove it--- In  general,  horizontal  displacement magnitude decreases from top to bottom of
% Propose to remove it--- the structure.
% Propose to remove it--- %
% Propose to remove it--- Significantly scattered horizontal displacement magnitudes are observed for four
% Propose to remove it--- different angles of incidence.
% Propose to remove it--- %
% Propose to remove it--- For  representative  points  on/above  the ground surface (points $A-D$, $I-K$),
% Propose to remove it--- maximum  horizontal  displacement  magnitudes given by $\theta = 60^{\circ}$ are
% Propose to remove it--- almost  $50\%$ greater than that of $\theta=10^{\circ}$ and approximately double
% Propose to remove it--- the horizontal responses given by $\theta=45^{\circ}$ and $\theta=80^{\circ}$.
% Propose to remove it--- %
% Propose to remove it--- From   Figure  \ref{fig_spatial_distribution_magnitude_SMR_inclination}  (b),  the
% Propose to remove it--- distinction   among   structural   vertical  displacement  caused  by  different
% Propose to remove it--- inclinations   is   still   very   noticeable,  though  not  so  significant  as
% Propose to remove it--- corresponding free field responses because of averaging effect from SSI.
% Propose to remove it--- %


The    deformed    shapes   of  SMR  at $t=0.4s$ for four scenarios are shown in
Figure~\ref{fig_inclination_SMR_deformation}.
%
\begin{figure}[!htb]
\centering
\includegraphics[width=1.0\textwidth]{Figures/SMR_inclination_snippet_t04.png}\enspace
\caption{\label{fig_inclination_SMR_deformation}  The  deformed shapes of SMR at
$t=0.4s$  for  incident  SV  wave  at different inclinations $\theta=10^{\circ},
45^{\circ}, 60^{\circ}$ and $80^{\circ}$.}
\end{figure}
%
In    the    cases  of  seismic  waves  at  inclinations $\theta=45^{\circ}$ and
$\theta=60^{\circ}$,  rocking  responses  of SMR are quite evident when compared
with  the  cases  of  almost vertical wave propagation ($\theta=10^{\circ}$) and
almost horizontal wave propagation ($\theta=80^{\circ}$).








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{SMR Excited with Variable Frequency Inclined SV Waves}
\label{SMR_Excited_with_Variable_Frequenct_Inclined_SV_Waves}


Keeping   incidence  angle constant, at $\theta=  60^{\circ}$, dynamic responses
of    an  SMR   under  different  frequencies  of  SV wave ($f=f=1Hz, 2.5Hz, 5Hz
{\rm and} 10Hz$) is investigated next.
%
Figures~\ref{fig_point_A_frequencies}              and                          
\ref{fig_point_A_frequencies_acceleration},  show  displacement and acceleration
responses at point A of SMR model.
%
\begin{figure}[!htbp]
\centering
\vspace*{-2cm}
\subfloat[Horizontal displacement]{
\includegraphics[width=0.9\textwidth]{Figures/frequency_ux.pdf}}\\
\subfloat[Vertical displacement]{
\includegraphics[width=0.9\textwidth]{Figures/frequency_uz.pdf}}
\caption{\label{fig_point_A_frequencies}
Displacement  response  of  point  A for scenarios with different frequencies of
incident SV wave: (a) Horizontal displacement (b) Vertical displacement.}
\end{figure}

\begin{figure}[!htbp]
\centering
\vspace*{-2cm}
\subfloat[Horizontal acceleration]{
\includegraphics[width=0.9\textwidth]{Figures/frequency_Ax.pdf}}\\
\subfloat[Vertical acceleration]{
\includegraphics[width=0.9\textwidth]{Figures/frequency_Az.pdf}}
\caption{\label{fig_point_A_frequencies_acceleration}
Acceleration  response  of  point  A for scenarios with different frequencies of
incident SV wave: (a) Horizontal acceleration (b) Vertical acceleration.}
\end{figure}

It  is noted that, again, free field response at the location of point A is also
shown for comparison purposes.
%
%
Significantly  variation in displacement and acceleration responses are produced
by incident SV wave at different frequencies.
%
The    largest   horizontal  displacement magnitude $0.30$m is observed  for the
case  of  frequency  of  $f=2.5Hz$  while  the  smallest horizontal magnitude of
$0.047$m for $f=1Hz$.
%
The  vertical  displacement  responses  varies  from  $0.02$m  for  $f=2.5Hz$ to
$0.085$m for $f=10Hz$.
%
SSI  effects  are almost negligible in the case of $f=1Hz$ due to long horizontal wave
length  of $1154$m. This observation follows similar observation made many years
ago by~\citet{Housner1957} for large stiff buildings.
%
Both  horizontal  and  vertical  displacements of SMR overlap with corresponding
free field response for $f=1Hz$.
%
Along  with  the  increase  of  incident  frequency,  SSI  effects  become  more
significant, especially for the vertical components of displacement and acceleration.
%
In the cases of $f= 2.5Hz$ and $f=5Hz$, horizontal response of SMR is still very
close  to  its free field counterpart, for both displacements and accelerations,
however the reduction of vertical response of SMR becomes more significant  for frequency
of $f=5Hz$,
% 
For relatively high frequency of $f=10Hz$, both horizontal and vertical response
of   SMR   are   significantly  different  from  free  field  modeling  in  both
displacements and accelerations.




The    spatial variation of displacements at the surface of free field model and
at  the  same  location  within SMR model, along the horizontal line through SMR
(i.e.   $x  \in  [-75$m,  $75$m$],  y  =  0m,  z=0$m),  at  $t=3.5$s  are  shown
in~Figure~\ref{fig_SSI_horizontal_spatial_distribution}.
% 
It is noted that SMR structure occupies space for $x\in [-15m, 15m]$, where flat
trace of displacements within a stiff structure is observed.
%
\begin{figure}[!htbp]
 \flushleft
 \subfloat[$f=1Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/1Hz_SSI_horizontal_comparison.pdf}}
 \subfloat[$f=2.5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/2Hz_SSI_horizontal_comparison.pdf}}
 \\
 \subfloat[$f=5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/5Hz_SSI_horizontal_comparison.pdf}}
 \subfloat[$f=10Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/10Hz_SSI_horizontal_comparison.pdf}}
 \\
 \caption{\label{fig_SSI_horizontal_spatial_distribution}
%
Spatial  variation  of  displacement  along  the horizontal axis at $t=3.5s$ for different incident wave frequencies
 (a) f$=1Hz$ (b) f$=2.5Hz$ (c) f$=5Hz$ (d) f$=10Hz$.}
\end{figure}
%
The  base slab averaging  is  observed for higher frequency, shorter wave length cases  of $f=5$Hz and
$f=10$Hz, while  it  is  almost negligible for incident waves at frequencies of
$f=1$Hz  or  $f=2.5$Hz  due  to the wavelength being longer that object size for
those low frequencies.
%
% For the possible amplification on the displacement magnitude. 
%




% OVDE
% OVDE
% OVDE
% OVDE
% OVDE
% OVDE




Similar spatial variation of displacement along the transverse axis (i.e. $x=0m,
y\in       [-75m,       75m],       z=0$m)       is      shown      in      
Figure~\ref{fig_SSI_transverse_spatial_distribution}.
%
%
%
\begin{figure}[!htbp]
 \flushleft
 \subfloat[$f=1Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/1Hz_SSI_transverse_comparison.pdf}}
 \subfloat[$f=2.5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/2Hz_SSI_transverse_comparison.pdf}}
 \\
 \subfloat[$f=5Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/5Hz_SSI_transverse_comparison.pdf}}
 \subfloat[$f=10Hz\ \theta=60^{o}$]{
 \includegraphics[width=0.5\columnwidth]{Figures/10Hz_SSI_transverse_comparison.pdf}}
 \\
 \caption{\label{fig_SSI_transverse_spatial_distribution}
Spatial  variation  of  displacement  along  the transverse axis at $t=3.5s$ for different
incident  wave frequency  (a) f$=1Hz$ (b) f$=2.5Hz$ (c) f$=5Hz$ (d) f$=10Hz$.}
\end{figure}
%
%
Since  the incident SV wave propagates within the XZ plane, uniform distribution
of both horizontal and vertical free field response along the transverse axis (Y
axis)             is        expected        and        presented        in      
Figure~\ref{fig_SSI_transverse_spatial_distribution}.
%
However,  the  existence  of SMR alters the original uniform distribution, and a
wave field in this, out plane of polarization direction.
%
Significant    wave  field  disturbance      effects  can be observed within the
structure  part  ($y  \in [-15m, 15m]$) in the cases of medium ($f=5Hz$) to high
frequency ($f=10Hz$).
%
In  other words, 3C dynamic response of soils surrounding the structure has been
induced  from  2C  excitation by an SV wave due to SSI and transverse wave field
disturbance effects.
%



Another            important           observation           from           Fig.
\ref{fig_SSI_transverse_spatial_distribution}(d)  is  that, although the reduction
of    displacement  amplitude  is  observed  within  the structure, in locations
where $y \in [-15m, 15m]$, near field motions close to the structure can be
amplified, for
example, motion within region $y\in \pm[25m, 50m]$ in this case.
%
This    implies    that there are potentially significant structure-soil-structure
dynamic effects for closely spaced  structures.
% to be constructed in the near
%field  of  current  SMR,  then  its  final  dynamic  response  might not be only
%influenced by the base averaging of itself.
% %
% The  amplified excitations caused by neighboring SSI system should also be taken
% into consideration.
% %
% The  claimed  benefits from base averaging should be utilized with care based on
% systematic analysis of structure-soil-structure interaction (SSSI).
% %
% Detailed  studies  on  this aspect is beyond the scope of this paper and will be
% reported in the future.














% Propose to remove it--- 
% Propose to remove it--- \begin{figure}[!htb]
% Propose to remove it--- \centering
% Propose to remove it--- \subfloat[Horizontal displacement]{
% Propose to remove it--- \includegraphics[width=0.40\textwidth]{Figures/spatial_mag_distribution_Ux_frequency_new.pdf}}
% Propose to remove it--- \subfloat[Vertical displacement]{
% Propose to remove it--- \includegraphics[width=0.60\textwidth]{Figures/spatial_mag_distribution_Uz_frequency.pdf}}
% Propose to remove it--- \caption{\label{fig_spatial_distribution_magnitude_SMR_frequency}        Spatial
% Propose to remove it--- variation  of  displacement  magnitudes of representative points for incident SV
% Propose to remove it--- wave with different frequencies.}
% Propose to remove it--- \end{figure}
% Propose to remove it--- 
% Propose to remove it--- 
% Propose to remove it--- 
% Propose to remove it--- The  spatial  variation  of  displacement magnitudes of representative points is
% Propose to remove it--- shown in Fig. \ref{fig_spatial_distribution_magnitude_SMR_frequency}.
% Propose to remove it--- %
% Propose to remove it--- The  large  scattering  among  results  of  different scenarios demonstrates the
% Propose to remove it--- important role of incident wave frequency (wavelength) in dynamic SSI.
% Propose to remove it--- %
% Propose to remove it--- Maximum horizontal displacement magnitude occurs with frequency $f=2.5Hz$, while
% Propose to remove it--- corresponding   vertical  displacement  magnitude  is  the  least  for  most  of
% Propose to remove it--- representative points.
% Propose to remove it--- %
% Propose to remove it--- High  frequency  excitation $f=10Hz$ gives rise to the largest vertical response
% Propose to remove it--- for many representative points (point $C-G$ and $I-K$).
% Propose to remove it--- %
% Propose to remove it--- It  is  also noted that, though magnitudes of free field motion are identical at
% Propose to remove it--- the  same  elevation, corresponding structural response (e.g., response of point
% Propose to remove it--- $A$, $B$ and $C$) shows noticeable spatial variation in the cases of $f=5Hz$ and
% Propose to remove it--- $10Hz$.
% Propose to remove it--- %
% Propose to remove it--- This implies that amplified structural rocking response, greater than free field
% Propose to remove it--- cord  rotation,  can  be  brought  by  SSI when the wavelength becomes short and
% Propose to remove it--- compatible to the dimension of structural object.
% Propose to remove it--- %



The  deformed  shapes  of  SMR  for  four frequency scenarios  at $t=0.3s$ with different
frequencies are shown in Fig. \ref{fig_frequency_SMR_deformation}.
%
%
\begin{figure}[!htb]
\centering
\includegraphics[width=1.0\textwidth]{Figures/SMR_frequency_snippet_t03.png}\enspace
\caption{\label{fig_frequency_SMR_deformation}  The  deformed  shapes  of SMR at
$t=0.3s$ for four scenarios.}
\end{figure}
%
%
The    aforementioned        wave   field  disturbance      effects  are clearly
visible for the low wave length, high frequency case of $f=10Hz$.
%
The  existence  of  local  structure  has  significantly  altered the near field
seismic  wave  due  to  strong  SSI effect, since wave lengths are shorter than
the dominant dimension of the structure.
% %
% In  this  case,  free  field  motion  can not be directly imposed onto the local
% structure.
% %
% Realistic  structural  response  can  be  captured through integral SSI analysis
% using  WPF-DRM,  which inputs the inclined wave field in a dynamically consistent
% way.
% %







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Summary}
\label{Summary}


Presented was Wave-Potential-Formulation  (WPF) -- Domain Reduction Method (DRM)
approach,  called  WPF-DRM,  for  solving  Earthquake Soil Structure Interaction
(ESSI) problems in layered ground excited by inclined incident seismic waves.
%
%  Assumptions adopted in many previous studies, for examples, rigid foundation and
%  homogeneous ground are removed in the proposed methodology.
%  %
%  Compared   with   conventional  substructure  method,  difficulties  of  solving
%  foundation wave scattering and impedance function are also circumvented.
%
Developed  WPF-DRM  methodology  removes a need for many simplifying assumptions
that are used in ESSI analysis, for example rigid foundation and
homogeneous ground  assumption.
%
In  addition,  difficulties    of    solving  for foundation wave scattering and
impedance function are also circumvented.
%
Most importantly, developed WPF-DRM method can be used with nonlinear, inelastic
soil, interface and structural material behavior.
%
WPF-DRM  is  verified  through recoverability test (i.e., resumption behavior) of
free field motions in a layered ground under incident SV waves.
%



Application  of  WPF-DRM  is  illustrated  by  analyzing  a  problems of an ESSI
response of a deeply embedded structure, a small modular reactor (SMR). 
%
Focus was on analyzing influence of a number of differently inclined plane waves
and a number of different wave frequencies, wave lengths.
%
It is noted that free   field  responses  for  incident  SV  waves  of  varying 
frequencies and inclinations show significant differences between free field and SSI.
%
For free field response, surface  rolling movement pattern, Rayleigh waves are captured.
%
This  is  different  from  typically assumed, vertically propagating wave field,
and  differences in SSI behavior are quite significant especially for medium and
high frequency inclined incident wave.
% 
% 
%  To  illustrate  WPF-DRM,  dynamic  response  of  SMR deeply embedded in a layered
%  ground under incident SV wave is modeled.
%
For  sensitivity  study, a monochromatic SSI response of SMR under incident SV wave
with different frequencies and inclinations is analyzed.
%
It is found that SSI effects are more prominent considering seismic motions with
non-vertical incidence and relatively high frequency, low wave lengths.
%
The vertical structural response is significantly influenced by the inclinations
of incident wave.
%
The  vertical  structure  response  can  vary  by  a  factor  of 7 for different
inclinations.
%
Compared   with  almost  vertical  wave  field  ($\theta=10^{\circ}$  inclination)  and
almost horizontal  wave  field  ($\theta=80^{\circ}$  inclination),  more  significant  structural
rocking response is observed in the cases of inclination $\theta=45^{\circ}$ and
$\theta=60^{\circ}$.
%
The  structural  response is almost identical to corresponding free field motion
in the case of low frequency $f=1Hz$ and long wavelength $1155$m.
%
As  the  frequency increases, structural response is different from free field
counterpart because of ``base averaging'' of ``ironing out'' effects.
%
This is particularly significant  for  high  frequency  incident  wave  ($f=10Hz$) where wavelength is
comparable to structural dimension, with observation of significant reduction in structural response.
%
Observed  are  also  wave  field  disturbance      effects   in the sense that near field motion is
notably altered by the existence of embedded structure, for example, in the case
of $f=10Hz$.
%  %
%  Large  scattering  in  the  response  of representative points from case to case
%  urges  the  necessity for more realistic modeling of spatially varying, inclined
%  motion field for dynamic SSI analysis.
%
Presented  examples provide  evidence of significance of  modeling uncertainties
that are introduced by the assumption of uniform, vertically propagating wave field.
%

%  It  is  noted  that  WPF-DRM  can  also  be  used to study many other effects and
%  influence  mechanisms in SSI system under inclined incident seismic waves, such
%  as topography, nonlinearity and other geological heterogeneity.
%  % 
%  More future work would be done on these aspects.
%  %








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgments}

This  work  was  supported  by  the University of California, by private funding
sources and in small part by the US-DOE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliography{refmech}

\bibliographystyle{unsrtnat}
% \bibliographystyle{elsarticle-num}


% % 
% % 
% % 
% % 
% % 
%  MOVED back to the main text
% % 
% \appendix
% 
% \section*{Appendix}
% 
% \renewcommand{\theequation}{A.\arabic{equation}}
% 
% % \subsection*{List of symbols}
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection*{Layered wave transformation matrix}
% \label{Layered_wave_transformation_matrix}
% 
% \begin{equation}
% \begin{gathered}
% \gamma_{\alpha_m} = \sqrt{(c/\alpha_m)^2-1} \ \ when \ \ \alpha_m \leq c \\
% \gamma_{\alpha_m} = -i\sqrt{1-(c/\alpha_m)^2} \ \ when \ \ \alpha_m > c \\
% \gamma_{\beta_m} = \sqrt{(c/\beta_m)^2-1} \ \ when \ \ \beta_m \leq c \\
% \gamma_{\beta_m} = -i\sqrt{1-(c/\beta_m)^2} \ \ when \ \ \beta_m >c 
% \end{gathered}
% \label{eq:angle_formulation}
% \end{equation}
% 
% % \resizebox{\linewidth}{!}{
% % $\boldsymbol{E_m} =
% % \begin{bmatrix}
% % -(\alpha_m/c)^2 & 0 & -\theta_m \gamma_{\beta_m} & 0 \\
% % 0 & -(\alpha_m/c)^2\gamma_{\alpha_m} & 0 & \gamma_m \\
% % -\rho_m\alpha_m^{2}(\theta_m-1) & 0 & -\rho_mc^2\theta_m^2\gamma_{\beta_m} & 0 \\
% % 0 & \rho_m\alpha_m^2\theta_m\gamma_{\alpha_m} & 0 & -\rho_mc^2\theta_m(\theta_m-1)
% % \end{bmatrix}$}
% % \\
% 
% \begin{equation}
% \resizebox{0.89\linewidth}{!}{
% $\boldsymbol{E_m} =
% \begin{bmatrix}
% -(\alpha_m/c)^2 & 0 & -\theta_m \gamma_{\beta_m} & 0 \\
% 0 & -(\alpha_m/c)^2\gamma_{\alpha_m} & 0 & \gamma_m \\
% -\rho_m\alpha_m^{2}(\theta_m-1) & 0 & -\rho_mc^2\theta_m^2\gamma_{\beta_m} & 0 \\
% 0 & \rho_m\alpha_m^2\theta_m\gamma_{\alpha_m} & 0 & -\rho_mc^2\theta_m(\theta_m-1)
% \end{bmatrix}$}
% \label{eq:matrix_Em}
% \end{equation}
% 
% where $\theta_{m} = 2(\beta_m/c)^2$. 
% \\
% \begin{equation}
% \resizebox{0.89\linewidth}{!}{
% $\boldsymbol{D_m} =
% \begin{bmatrix}
% -(\alpha_m/c)^2cosA_m & i(\alpha_m/c)^2sinA_m & -\theta_m \gamma_{\beta_m}cosB_m & i\theta_m\gamma_{\beta_m}sinB_m \\
% i(\alpha_m/c)^2\gamma_{\alpha_m}sinA_m & -(\alpha_m/c)^2\gamma_{\alpha_m}cosA_m & -i\theta_msinB_m & \theta_mcosB_m \\
% -\rho_m\alpha_m^{2}(\theta_m-1)cosA_m & i\rho_m\alpha_m^2(\gamma_m-1)sinA_m & -\rho_mc^2\theta_m^2\gamma_{\beta_m}cosB_m & i\rho_mc^2\theta_m^2\gamma_{\beta_m}sinB_m \\
% -i\rho_m\alpha_m^{2}\theta_m\gamma_{\alpha_m}sinA_m & \rho_m\alpha_m^2\theta_m\gamma_{\alpha_m}cosA_m & i\rho_mc^2\theta_m(\theta_m-1)sinB_m & -\rho_mc^2\theta_m(\theta_m-1)cosB_m
% \end{bmatrix}$}
% \label{eq:matrix_Dm}
% \end{equation}
% \\
% where $A_m=k\gamma_{\alpha_m}d_m$ and $B_m=k\gamma_{\beta_m}d_m$. 
% 
% 
% 
% 






\end{document}
%\grid
