\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{elsarticleharv}
%\bibliographystyle{elsarticlenum}
%\usepackage[numbers,round,colon,authoryear]{natbib}
% \usepackage[number,sort&compress]{natbib}
%% Water mark
%\usepackage{type1cm}
%\usepackage{esopic}
%\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 WPFDRM, is verified using analytic
solution for a free field response of layered ground subjected to inclined
incident waves.
%
%Compared with conventional substructure method, WPFDRM circumvents the
%additional difficulties of solving foundation wave scattering and impedance
%function.
%
Developed WPFDRM 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
nearfault regions and for structures with largeplan 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 subproblems:
\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 nonvertically
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 inplane 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 halfspace 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 viscousspring artificial boundary put forward
% by~\citet{Liu2006c} and \citet{Du2006} can efficiently simulate the radiation
% damping and the elastic resumption behavior of semiinfinite medium.
% %
%
%
Using developed viscousspring artificial boundary, various SSI and
rockstructure 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 soilstructure 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 (WPFDRM) 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_SoilStructure_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 threestep substructure method, twostep WPFDRM 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 WPFDRM 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 frequencywave 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 soilstructure
system is then performed using effective earthquake forces that are applied to a
single layer of finite elements surrounding soilstructure system, so called
DRM layer.
%
The only waves that are radiated from the soilstructure 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 soilstructure 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 LeviCivita 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)^21} \ \ 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)^21} \ \ 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)^21} & \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)^21} & \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 postcritical P and SV waves travel together simultaneously along the free surface; Love wave comes from PostCritical 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(12\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)]\\
+[12(\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^{(m1)}
=
\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_m1) & 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_m1)
% \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_m1) & 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_m1)
\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_m1)cosA_m & i\rho_m\alpha_m^2(\gamma_m1)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_m1)sinB_m & \rho_mc^2\theta_m(\theta_m1)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^{(m1)}$ 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^{(m1)}
=
\boldsymbol{G_{m}}S^{(m1)}
\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 $(m1)^{th}$ interface $S^{(m1)}$ 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^{(n1)} = \prod_{i=1}^{n1}\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}^{n1}\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 inplane $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_{m1}^{''}+\Delta_{m1}^{'}\\
\Delta_{m1}^{''}\Delta_{m1}^{'}\\
\omega_{m1}^{''}\omega_{m1}^{'}\\
\omega_{m1}^{''}+\omega_{m1}^{'}
\end{bmatrix}
=
\boldsymbol{D_{m1}^{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 $n1$ 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 KelvinVoight 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 soilstructure 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 substructuring
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 (WPFDRM) 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 WPFDRM method is implemented in the RealESSI
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://realessi.info/}.
%
All numerical examples presented here are analyzed using RealESSI 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 ThomsonHaskell 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, nonfree 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 WPFDRM.
%
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 nonexistent 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 WPFDRM simulation and
analytical solutions.
% %
% This verifies and confirms the capability of WPFDRM 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 SoilStructure Model}
\label{Deeply_Embedded_SoilStructure_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
eightnode brick elements, at least 10
nodes per wavelength should be used \citep{Watanabe2015}.
%
The time step length $\Delta t$ is limited by CourantFriedrichsLewy
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, eightnode 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{local87}.
%
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 $AD$, $IK$),
% 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 structuresoilstructure
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 structuresoilstructure 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 $CG$ and $IK$).
% 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 WPFDRM, which inputs the inclined wave field in a dynamically consistent
% way.
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Summary}
\label{Summary}
Presented was WavePotentialFormulation (WPF)  Domain Reduction Method (DRM)
approach, called WPFDRM, 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 WPFDRM 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 WPFDRM method can be used with nonlinear, inelastic
soil, interface and structural material behavior.
%
WPFDRM is verified through recoverability test (i.e., resumption behavior) of
free field motions in a layered ground under incident SV waves.
%
Application of WPFDRM 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 WPFDRM, 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
nonvertical 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 WPFDRM 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 USDOE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{refmech}
\bibliographystyle{unsrtnat}
% \bibliographystyle{elsarticlenum}
% %
% %
% %
% %
% %
% 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)^21} \ \ 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)^21} \ \ 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_m1) & 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_m1)
% % \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_m1) & 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_m1)
% \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_m1)cosA_m & i\rho_m\alpha_m^2(\gamma_m1)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_m1)sinB_m & \rho_mc^2\theta_m(\theta_m1)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