% For Earthquake Engineering and Structural Dynamics ==============================================
\documentclass[times]{eqeauth}
\usepackage{moreverb}
\usepackage[colorlinks,bookmarksopen,bookmarksnumbered,citecolor=red,urlcolor=red]{hyperref}
\newcommand\BibTeX{{\rmfamily B\kern-.05em \textsc{i\kern-.025em b}\kern-.08em
T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}
\def\volumeyear{2016}
%==================================================================================================

\usepackage[T1]{fontenc} % Use 8-bit encoding that has 256 glyphs
\usepackage[utf8]{inputenc} % Required for including letters with accents
% \usepackage[final]{graphicx} % Required for including images
\usepackage{enumitem} % Required for manipulating the whitespace between and within lists
\usepackage{subcaption} % Required for creating figures with multiple parts (subfigures)
\usepackage{amsmath,amssymb,amsthm} % For including math equations, theorems, symbols, etc

\usepackage{cite}
% \usepackage{natbib}

% \linespread{1.5}



\usepackage{myListings}
\usepackage{mymath}
\usepackage{myunits}
\usepackage{pstricks}


\newcommand{\ttod}{3$\times$1-D }

\usepackage[normalem]{ulem}

\definecolor{Gray}{gray}{0.90}   
\definecolor{josecommentcolor}{HTML}{35764D}   
\newcommand{\citationneeded}{{\color{red}[citation needed]}}
\newcommand{\done}{ {\color{blue} \Large \checkmark}}


%Uncomment to include comments
\newcommand{\jcomment}[1]{ {\color{josecommentcolor} \textbf{JA: #1 }}}
% \newcommand{\proposedchange}[2]{\sout{#1} {\color{blue} #2}}
% \newcommand{\proposedchange}[2]{{\color{blue} #2}}
\newcommand{\proposedchange}[2]{ #2}



\begin{document}

%----------------------------------------------------------------------------------------
% HEADERS
%----------------------------------------------------------------------------------------


\runningheads{J.~A.~Abell; B.~Jeremi\'c}{ESSI Modeling of NPPs for Near-field Events}

% \title{The Importance of Capturing the 3-D Character of Seismic Input for Earthquake-Soil-Structure Interaction Analysis of Nuclear Power Plants}

\title{Earthquake Soil Structure Interaction of Nuclear Power Plants, differences in response to  3-D, 3$\times$1-D, and 1-D excitations.}

\author{Jos{\'e} A. Abell$^{1}$\corrauth, 
        Neboj{\v s}a Orbovi{\'c}$^{2}$,
        David B. McCallen$^{3}$,  
        Boris Jeremi\'c$^{3,4}$}

\address{
$^1$ Facultad de Ingeniería y Ciencias Aplicadas, Universidad de los Andes, Santiago, Chile \\
$^2$ Canadian Nuclear Safety Commission, Ottawa, ON. Canada \\
$^3$ Lawrence Berkeley National Lab, Berkeley, California, USA \\
$^4$ Department of Civil and Environmental Engineering, University of California, Davis, California, USA 
}
\corraddr{José Abell, Facultad de Ingeniería y Ciencias Aplicadas, Universidad de los Andes, Mons. Álvaro del Portillo 12.455, Las Condes, Santiago, Chile. 7620001}

\begin{abstract}
In soil-structure interaction (SSI) modeling of systems subjected to earthquake motions it is classically assumed that the in-coming wave-field, produced by the earthquake, is unidimensional and vertically propagating. This work explores the validity of this assumption by performing earthquake-soil-structure interaction (ESSI) modeling, including explicit modeling of sources, seismic wave propagation, site and structure. The domain reduction method, DRM, is used to couple seismic (near-field) simulations with local SSI response.  The response of a generic Nuclear Power Plant (NPP) model computed using full ESSI simulations is compared with the current state-of-the-art method of de-convolving in depth the (simulated) free-field motions, recorded at the site of interest,  and assuming that the earthquake wave-field is spatially unidimensional.  Results show that the 1-D wave-field assumption does not hold in general. It is shown that the way in which full 3-D analysis results differ from those which assume a 1-D wave-field is dependent on fault-to-site geometry and motion frequency-content. It is argued that this is especially important for certain classes of soil-structure systems of which nuclear power plants subjected to near-field earthquakes are an example.

\end{abstract}

\keywords{Domain Reduction Method; Soil-Structure Interaction; Wave Propagation}

\maketitle

\vspace{-6pt}



\section{Introduction}


In soil-structure interaction (SSI) modeling of systems subjected to earthquake motions it is usually assumed, often implicitly by just using a particular software, that the in-coming earthquake wave-field is a unidimensional, vertically propagating plane-wave. Its been recommended for 40 years now `\emph{that soil-structure interaction studies should not be limited to seismic excitations with vertical incidence}'\cite{Luco1978}. Still, recent guidelines like the (2012) NEHRP Soil-Structure-Interaction for Building Structures recommendations\cite{NEHRP2012}, implicitly assume a vertically incident plane-wave input motion for virtually all modeling recommendations and example applications. The plane seismic wave-field assumption is useful because it tremendously simplifies the boundary conditions of the problem, and because it enables the direct use of recorded motions in seismic response analysis. On the other hand it introduces a modeling idealization, the implications of which are still poorly understood and still an area of active research. 

Conceptually, assuming a 1-D seismic wave-field implies that horizontal components of surface motion are produced by vertically-propagating horizontally-polarized shear waves while the vertical component is produced exclusively by vertically-propagating laterally-constrained compressional waves. It is known that the vertical component of surface motions is produced by a combination of inclined compressional and shear-wave as well as surface waves. This complexity in the earthquake wave field and its potential to significantly affect the response of structures subjected to it was recognized by Housner as early as 1956\cite{Housner1956}. He argued that, what we now call, wave-passage effects could reduce demands on structures with large in plan extents and rigid foundations (this is also known as the `base averaging effect'). Later, the incidence angle of plane earthquake wave-fields has been linked to amplified rocking as well as torsional response of structures\cite{Luco1976,Luco1982,Trifunac2015,Trifunac2016b,Trifunac2016a,Mavroeidis2017}. It has also been shown that the same wave-passage effect, linked before to response reductions due to base-averaging, can amplify structural response with respect to the vertically incident case\cite{Trifunac2011} for long structures. These extremes, response reduction in one case and amplification in another, hint that the problem of seismic SSI is probably best studied case-by-case: a specific structure on a specific site experiencing an earthquake from a specific source. 

Early on, response of structures to an obliquely incident plane body wave-field or surface wave fields including SSI was addressed using analytical methods\cite{Trifunac1974,Luco1978,Luco1980,Luco1982}. With these methods a basic understanding of the possible effects of wave-field obliquity was established. It was recognized that buildings with large in-plan extents are sensitive to wave-passage effects and can experience enhanced rocking and torsional response. Furthermore, high frequency components of structural response are more sensitive to the three-dimensional characteristics of wave propagation than low frequency ones. More insight has been gathered on the subject recently by Trifunac et al. \cite{Trifunac2015,Trifunac2016a,Trifunac2016a}, where a detailed 2-D non-linear numerical model subject to different phases of plane waves was used to establish the importance of wave incidence angle and phase type in the response of buildings. Such studies further confirm the inadequacy of assuming vertical propagation as input motions.

In the light of this, it is apparent that nuclear power plants (NPPs) have characteristics which make them especially sensitive to the three-dimensional character of seismic wave-fields. First, they are stiff structures usually with short fundamental periods. Second, NPPs have large lengths in plan and are usually founded on stiff soils of high strength, implying that several wavelengths of seismic motion will `fit' within the area of the NPPs at frequencies of input of engineering interest. Therefore, NPPs can transmit high-frequency motions to its components and equipment very efficiently. Third, the contents of NPPs are sensitive to accelerations (such as piping, heavy equipment and their connections) and to displacements such as fluid sloshing and differential movement of interconnected components.  

% The three-dimensional character of the incoming seismic wave-field results in different parts of the NPP foundation being out-of-phase (sometimes referred to as kinematic interaction in the SSI literature). Quantification of effects such as base averaging, often leveraged for design but never fully justified, depend on being able to capture this kinematic interaction accurately, which cannot be done if the seismic input method is incapable of producing these out-of-phase motions in the first place. 

Near-field earthquake events, that is events which occur near an NPP site (typical distance from site to source less than 20 kilometers), pose a special kind of threat to NPPs.  Despite the fact that NPPs would not likely be emplaced near a potentially active fault, NPPs have been built unknowingly near active faults that have been discovered after construction is completed \cite{Sykes2008,Hardebeck2010}. Data for events generated at these nearby faults is usually scarce to non-existent, such that the only feasible way to obtain detailed, three-dimensional ground motions is through numerical simulation. For nearby events, the wave-front in the near-field shows the curvature of initial propagation, such that a plane-wave approximation of this surface might not be satisfactory. Additionally, the high frequency components of motion are less damped because there has not been enough propagation distance to produce significant attenuation. Therefore, small magnitude nearby events can produce high demands on a structure and its contents for some performance metrics, which are hard to anticipate, and there are reasons to suspect that plane-wave approximations of seismic field--even inclined ones--are not suitable.  Advantageously, modeling near-field events is computationally less expensive than far-field events due to the size of the domain that needs to be included in the wave propagation computation, well within the reach of modern small-scale HPC clusters or cloud-based resources.  


Therefore, in the present work, the adequacy of vertically incident plane-wave approximations to seismic wave-fields from near-field sources is evaluated in terms of the differences in response produced on a model of a generic nuclear power plant. This is accomplished by using seismic simulations of regional near-field wave propagation coupled with a detailed site and structure finite element model through the domain-reduction method.

The domain reduction method (DRM), presented in the pioneering work of Bielak et al.\cite{Bielak2003}, allows coupling of realistic seismic waveform simulations with time-history analysis of a soil-structure system in such a way that the full three-dimensional character of seismic-waves can be propagated into a finite element model without further assumptions. Its usage requires detailed knowledge of the incident seismic-wave field on all points located along a boundary of finite-elements which encompasses the soil-structure system with an adequate spatial resolution (less than one meter for high-frequency applications).  Such a dense characterization of the in-coming wave field at the boundary nodes can only be obtained, as was envisioned by Bielak et. al, by simulation of the propagation phenomenon from source to site. This can come at high computational cost.

% Validated, in this case, means using a model and simulation of real source rupture, crustal propagation and local site conditions whose response has been shown to explain instrumental measures, and extrapolate these results to the locations of DRM nodes. 

% Even if the data necessary for using the DRM correctly is obtainable, storing and accessing DRM motion data is a challenge for a realistic case. For example, a small model for SSI might have on the order of tens of thousands of nodes along the DRM boundary. For each node, 3 components of displacement and acceleration time histories need to be provided. A seismic record would have a duration of significant shaking anywhere from a few tens of seconds to hundreds of seconds for large subduction-zone events. Finally, the sampling rates needed to appropriately represent the involved waveforms are in the order of 100s of Hertz. This results in datasets with sizes on the order of gigabytes.

\proposedchange{[Details on data storage and access was removed]}{}
% So if there are $N_{\mathrm{nodes}}=10,000$ nodes, at which 3 components of motions (accelerations and displacements) need to be stored as a double-precision number, which time-history lasts 30 seconds at a sampling rate of 100Hz, the total amount of data needed to store and access is:
% \begin{align}
% \mathrm{DataSize} &= 2 \times 3 \bfunit{components}{node} \times 10,000\bunit{nodes} \times 30 \bfunit{sec}{component} \times \hdots \nonumber \\
%  &\hdots \times 100 \bfunit{datapoint}{sec} \times 8 \bfunit{bytes}{datapoint} \nonumber  \\
%  &= 1,440,000,000 \bunit{bytes} = 1.440 \bunit{GB} 
% \end{align}

% For this amount of data, an inadequate storage and access pattern can seriously slow down computations especially in parallel processing mode. Furthermore, this number grows with the inverse-square of the discretization length (number of DRM-boundary nodes grows with the inverse-square of the grid spacing length). This challenge has provided another barrier to the adoption of the method, which made it hard apply at the time of proposal of the DRM, but are now quite easy to deal with. 

% Even if the data necessary for correctly using the DRM is obtainable, storing and accessing DRM motion data is a challenge for a realistic case. Even relatively simple models with node counts within a hundred thousand can require datasets on the order of a few gigabytes. For this amount of data, an inadequate storage and access pattern can seriously slow down computations, especially in parallel processing mode. The challenge posed by these technical issues has resulted in another barrier to the adoption of the method, which its application hard at the time the DRM was originally formulated, but are now easier to deal with. 

To get around this, approaches to using DRM for SSI simulations in the past have \proposedchange{either assumed}{simplified the problem to}  a 1-D wave-field \proposedchange{again to match some}{matching a}  seismic record. Refer, for example, to Psarropolous (2007)\cite{Psarropoulos2007}, Jeremi\'c (2009)\cite{Jeremic2009}, Kontoe (2012)\cite{Kontoe2012}, Tripe (2013)\cite{Tripe2013}, Tafazzoli (2013)\cite{Tafazzoli2013}, Zhong (2014)\cite{Zhong2014} \proposedchange{Isbiliroglu (2015)$[10]$,}{, and Solberg}\cite{Solberg2016}. Alternatively, some studies do compute fully three-dimensional DRM motions from a wave-propagation direct simulation approach but with a limited bandwidth. For example, Bielak et al. \cite{Bielak2003} and Yoshimura et. al. \cite{Yoshimura2001}. Usage of DRM motions developed using a fully 3D source and heterogeneous  models capable of resolving high frequency motions have been limited by computational demands.  \proposedchange{}{With recent advancements in High Performance Computing (HPC), the ability to exploit the full potential of the DRM is continuously increasing as a result of the capability for executing ever larger, higher resolution regional scale ground motion simulations that can model extended faults and resolve frequencies of engineering interest. The work of Isbiliroglu et. al. (2015) \cite{Isbiliroglu2015} is an example of utilization of DRM in conjunction with a 3-D extended fault model in a spatially heterogeneous medium with an application to the quantification of the ef ects of clusters of simplified buildings on ground motion as well as those of soil-structure interaction in individual structures and of structure-soil-structure interaction between buildings, for frequencies up to 5 Hz.
}

A more `classical' approach to seismic SSI involves inputing the (1-D) wave-field as a stress input at a given depth under the model and capturing out-going motions with Lysmer-type boundaries \cite{kuhlemeyer1973,Kouroussis2011}. Notable examples of this approach include: Elgamal et al (2008)\cite{Elgamal2008}, and Ostadan et al (2008) \cite{Ostadan2008}, to name a few. 




%For example, the recent near-field earthquake ($M_w = 5.8$ at a distance of $16 \kmts$) that hit the North Anna NPP produced virtually no damage to the NPP.


In this scenario, physical modeling of the complete ESSI problem is possible using DRM to couple a detailed local model of site and structure with high-performance, validated models of earthquake sources and regional wave propagation properties. This is both practical and necessary to show adequate seismic response of structure, especially NPPs, to these events. 

The main idea in this article is to determine, using DRM, if the response of a structure excited by a plane and vertically propagating input wave-field is different from the one excited by the full 3-D wave field, even when these fields are identical at the point of emplacement of the structure. For this purpose, it will suffice to use a linear model of the NPP and site as well as simple point source for the earthquake because the differences, it is postulated, will arise even in this basic setup as a basic feature of the physical phenomenon at hand.

%The difficulties that the 1-D assumption encounters when predicting non-linear site-specific, i.e. free-field, response has already been well established (for example in \cite{Rathje2015} and \cite{Regnier2016}) and it can only be expected that it will be worse for cases involving response of non-linear soil-structure systems. 

A simplified NPP model was created in the Real ESSI Simulator System \cite{Real_ESSI_Simulator} developed at the University of California at Davis and Lawrence Berkeley National Lab. Earthquake motions will be input into the system using equivalent forces obtained via the Domain Reduction Method. These forces were obtained by evaluating the response of the seismic wave propagation problem at the nodes located at the DRM boundary, for any given source. Earthquakes will be modeled as dip-slip point sources, at varying depths and with an adequate frequency content.  For the motions produced at each source depth, three types of DRM input motions were developed: `Full 3-D' motions which correspond to the direct input of the motions produced by the seismic simulations into the model, `1-D'  motions which are a one-component, vertically propagating plane-wave approximation to the Full 3-D motions, and `\ttod' motions which are a three-component vertically propagating plane-wave approximation to the Full 3-D motions. Both plane wave approximations were obtained by performing 1-D motion de-convolutions of the recorded motion at the NPP site.

% Motions due to a `Full 3-D' wave field, eg. the results from the 3-D seismic simulations are used 

From the simulated 3-D seismic wave field, ie. the `free-field' surface motions computed at the NPP site using the seismic simulation,  vertically incident plane wave-field approximations are developed by de-convolving each component of motion independently. De-convolutions are carried out in the SHAKE91\cite{SHAKE91} program for 1-D site response analysis (see considerations by Mejía and Dawson \cite{Mejia2006}). For horizontal components, a vertically incident, horizontally polarized pure shear wave-field is assumed, using shear-wave speed and density assumed for the site. For the vertical component, a laterally-constrained 
vertically-propagating compression wave is assumed, hence the constrained modulus of elasticty and 
other properties. When only the horizontal component is used (as in many analysis cases) the resulting wave-field is called, herein, a `1-D' wave. When all three components are de-convolved and input simultaneously the resulting wave-field is called `\ttod' wave since 3 components are used, but they only depend on the depth and time. Finally, using the complete seismic wave is termed `Full 3-D' wave-field. 

This approach differs from previous approaches to SSI analysis in several ways. First, motions come from fully three-dimensional source, are propagated and coupled perfectly into the local FEM model using DRM. This means that no-plane wave assumption is made and all wave phases are inherently considered. Second, frequencies up to $10\Hz$ are resolved, which, while possibly insufficient to completely characterize NPP hazard, can be accomplished with modest computational resources and is enough to prove the point of the article. Last, the goal is to compare  locally-equivalent vertically incident plane-wave approximations of wave field, a state of practice, with the Full 3-D approach to establish the consequences of such assumptions. The combination of all these modeling features serving this goal is the main novelty of this article.

The rest of the article develops as follows. Section 2 elaborates on the methods used to generate the full 3-D wave field from point source earthquake excitations and following through the propagation of the site, as well as the details regarding the deconvolution process used to develop 1-D and \ttod motions. Section 3 explains the `local' model used to represent the generic nuclear power plant and the local site used herein. Section 4 presents the results obtained in terms of NPP response at different points. The response from a 3-D wave-field is compared with 1-D wave-field, and is followed by comparing full 3-D with \ttod. Section 5 presents the main analysis and discussion of the results in terms of the different modeling assumptions. Finally, section 6 presents the main conclusions of this work.


\section{Generating DRM Motions}

In this study, earthquake motions are generated using a fourth-order finite-difference code SW4\cite{Sjogreen2011} which solves for the fully 3-D elastodynamics equations. The domain considered is a $10\kmts \times 5\kmts \times 5 \kmts$ box, with the NPP site located at the middle of its top surface. Figure \ref{fig:sw4_geometry} (left) shows this general setup. 


\begin{figure}[h]
\centering
    \includegraphics[width=\textwidth]{./img1/SW4_geom.pdf}
    \caption{Geometry of local crust at NPP site (left), and location of source earthquake relative to NPP site (right). }
    \label{fig:sw4_geometry}
\end{figure}


A horizontally layered elastic half-space is used to model the crust in the region around the NPP site. In full 3-D modeling with finite-differences, the model must contain both the source and the site. Ten homogeneous layers with shear wave speeds from $V_s=500 \mtsps$  at the top to $V_s=2500\mtsps$ at the bottom of the model represent the variation of crustal properties within the first $1 \kmts$ of depth. Properties at greater depths are assumed constant. \proposedchange{The compressional wave speed is set to twice the shear wave speed, while the density is varied from  $\rho = 2000 \kgmmm $ at the surface to $\rho = 2350 \kgmmm $ at $z = 1 \kmts$ and constant th
ereafter.}{Table~\ref{tab:properties} reports these and other properties used in all layers of regional model.} The chosen setup is meant to represent a generic site for a nuclear power facility, \proposedchange{usually located on sites featuring}{which usually feature} rock or high strength soils. 

\begin{table}
\centering
\renewcommand{\arraystretch}{1.0}% Tighter
\caption{Properties of the layers used in the regional model.}
\label{tab:properties}
\begin{tabular}{ccccc}
\hline
$V_p$ ($\mathrm{m/s}$) &
$V_s$ ($\mathrm{m/s}$) &
$\rho $ ($\mathrm{m/s}$) &
$z_1$ ($\mathrm{m}$) &
$z_2$ ($\mathrm{m}$)  \\ \hline
 1000 & 500 & 2000 & 0 &  50\\
 1894 & 947 & 2018 & 50 & 100 \\
 2265 & 1132 & 2035 & 100 & 200 \\
 2789 & 1394 & 2070 & 200 & 300 \\
 3191 & 1595 & 2105 & 300 & 400 \\
 3530 & 1765 & 2140 & 400 & 500 \\
 3828 & 1914 & 2175 & 500 & 600 \\
 4098 & 2049 & 2210 & 600 & 700 \\
 4347 & 2173 & 2245 & 700 & 800 \\
 4578 & 2289 & 2280 & 800 & 900 \\
 4795 & 2397 & 2315 & 900 & 1000 \\
 5000 & 2500 & 2350 & 1000 & $\infty$ \\ \hline \\
\end{tabular}

$V_p$: speed of primary waves, $V_s$: speed of secondary waves, $\rho$ density, $z_1$ depth of the top of the layer, $z_2$ depth of the bottom of the layer. 
\end{table}


Reverse faulting earthquakes, modeled using point-source double-couple mechanisms, are placed at a distance of $2.5 \kmts$ away from the NPP site in the $x$ direction at source depths ($z_S$) of $550 \mts$, $850\mts$, and $1200 \mts$. This is depicted in Figure~\ref{fig:sw4_geometry} (right).  Large events might be modeled as a superposition of these types of sources with properly scaled magnitudes and delays given by the rupture propagation front. At least for a linear analysis, differences arising from comparing full 3-D wave-field with equivalent 1-D or \ttod ones would arise in even the case of having a point-source. 

% Capturing an adequate frequency content for the resulting motions is very important for a realistic analysis using point sources. Brune \cite{Brune1970,Brune1971} derived the shape of the source time-function based on physical considerations of the available stress to drive the near-fault acceleration of a rupturing fault. A smoothed version of the Brune source-time-function, with a corner-frequency of  was used herein

% This source time-functional shape, shown in Equation \ref{eq:05_brune_source_time_function}, used to drive the elasto-dynamic equations and results in a Fourier spectrum of simulated motions which matches the $\omega^{-2}$ decay of Fourier amplitudes which has been observed in recorded motions.
% \begin{equation}
% g(t,\,t_0\,,\omega_0) = \left\lbrace 
% \begin{array}{ll}
% 0 & t < t_0 \\
% 1 - \expon{-\omega_0 (t - t_0) }(1 + \omega_0 (t-t_0)) & t \geq t_0
% \end{array}
% \right.
% \label{eq:05_brune_source_time_function}
% \end{equation}
% In Equation \ref{eq:05_brune_source_time_function} $t_0$ is the rupture start time and $\omega_0 = 2 \pi f_0$ and $f_0$ is the so-called `source corner frequency' which is related to the magnitude of the event and the `stress-drop' $\Delta \sigma$ which is the change in shear stress across the fault which drives the acceleration on the fault. The source corner frequency is given by:
% \begin{equation}
% f_0 = 4.9 \times 10^{-9} V_s \pare{\dfrac{\Delta \sigma }{M_0}}^{1/3}
% \label{eq:05_brune_source_corner_frequency_SI}
% \end{equation}
% Where $f_0$ is in Hertz, $V_s$ in $\mtsps$,  $\Delta \sigma$ is in Pascals, and $M_0$ is the source moment in $\bunit{N \cdot m}$. 
% \begin{figure}[tb]
%     \centering
%     \includegraphics[width=0.7\textwidth]{./img2/corner_freq_with_magnitude.pdf}
%     \caption{Corner frequency vs. moment magnitude for different values of the stress drop.  }
%     \label{fig:05_corner_freq_with_magnitude}
% \end{figure}



% % Re word this. 

% % freq = 18
% % f0 = 2.86
% % fmax = 11.45
% % h = 4m  
% % lambdamax = 48m
% % fmax_mesh = 500 m/s / 48 = 10.41666667 Hz



% Figure \ref{fig:05_corner_freq_with_magnitude} shows plots of the corner frequency with magnitude for a Brune-type source. In this study it is desired to resolve motions up to a maximum frequency of $10 \Hz$. To produce comparable 3-D motions for sources at different depths, the amplitudes of motions at the site are normalized to produce a $1\bunit{cm}$ displacement in the horizontal $X$ direction.  Even though this study fixes this parameter, for real applications it is important to keep in mind the relation between stress-drop, corner-frequency and magnitude in the context of the known local geology and tectonic setting. Realistic sources with large magnitudes will have to be modeled as an ensemble of point sources, each paying special attention to this detail. The point here is to show that any differences in input wave assumptions stem right from the 3-D nature of the problem in the simplest of cases.  

% To avoid spurious high-frequency artifacts due to the discontinuous second derivative of the Brune source time-function at $t=t_0$, a smoothed version of the Brune source time-function is used:
% \begin{equation}
% g(t,\,t_0\,\omega_0) = \left\lbrace 
% \begin{array}{ll}
% 0 & t < t_0 \\
% 1 - \expon{-\omega_0 (t - t_0) } 
% \left[
%     1 + \omega_0 (t - t_0 ) + 
%     \dfrac{ (\omega_0 (t - t_0 ))^2}{2} + 
% \right. & \hdots \\
% \left.
%    \hdots -\dfrac{3 (\omega_0 (t - t_0 ))^3}{2\tau_0} + 
%     \dfrac{3 (\omega_0 (t - t_0 ))^4 }{2{\tau_0}^2}
%     -\dfrac{ (\omega_0 (t - t_0 ))^5}{2{\tau_0}^3}
% \right]
% & 0 < \omega_0 (t-t_0) < \tau_0 \\
% 1 - \expon{-\omega_0 (t - t_0)} (1 + \omega_0 (t-t_0)) & \omega_0(t - t_0) \geq \tau_0
% \end{array}
% \right.
% \end{equation}
% where $\tau_0 = 2.31$ is a constant chosen to achieve a balance between smoothness and closeness to the original Brune function. Figure \ref{fig:05_brune_sources} shows the time and frequency plots of the Brune source and its smoothed version for the corner frequency $f_0 = 10.0 \Hz$. The maximum frequency $f_{\max}$, a measure of the bandwidth of the source-time-function, can be defined as the frequency at which the amplitude of the Fourier spectrum of the source-time-function drops below 5\% of its peak value. Numerical evaluation of this definition for $f_0 = 10.0\Hz $ yields a value of $f_{\max}$ of $43.4\Hz$ for the Brune source while for the \texttt{BruneSmoothed} it is $43.2\Hz$. The approximate rule of thumb $f_{\max} \approx 4 f_0$ can be used to evaluate the maximum frequency.

 % With this in mind, if the highest resolvable frequency is to be $10\Hz$, then the maximum corner frequency for the Brune source is $f_0 = 2.5 \Hz$, which is used in this work.

% \begin{figure}[tb]
%     \centering
%     \includegraphics[width=\textwidth]{./img2/brune.pdf}
%     \caption{Time and frequency plots for a Brune and BruneSmoothed source time-functions both for $f_0 = \frac{\omega_0}{2\pi} = 10\Hz$. }
%     \label{fig:05_brune_sources}
% \end{figure}



% \begin{figure}[tb]
%     \centering
%     \includegraphics[width=\textwidth]{./img2/3dsite_symm.png}
%     \caption{FEM mesh of the site. }
%     \label{fig:05_mesh}
% \end{figure}



\proposedchange{[Details of Brune source model....]}{Capturing an adequate frequency content for the resulting motions is very important for a realistic analysis using point sources. Brune \cite{Brune1970,Brune1971} derived the shape of the source time-function based on physical considerations of the available stress to drive the near-fault acceleration of a rupturing fault. A smoothed version of the Brune source-time-function, with a corner-frequency of $f_0 = 10\Hz$ was used herein}. The seimic source model therefore consists of single double-couple point-sources at three source depths ($z_s = \lbrace 550 \mts,\,850\mts,\,1200\mts\rbrace$ ). Faulting mechanism assumes a reverse-fault at 45\dgr dipping angle towards the direction of the NPP. The simulated displacement response of the elastic medium is recorded at the site of the geometric center of NPP foundation, which is then used to produce equivalent 1-D and \ttod wave-fields for comparison. Motions are also recorded at the locations of the FEM mesh nodes on the DRM boundary, and used to develop DRM forces. Because the seismic model grid points and FEM mesh nodal points don't exactly match in space as a result of different discretization sizes, a nearest-neighbor approach is used to assign motions from the seismic simulation into each DRM node. Free field studies using the FEM to predict site motions from the DRM input, reported elsewhere\cite{phdAbell2016}, show good agreement between target motions at the site computed with DRM. All seismic simulations were carried out using the finite-difference code SW4, on the Lawrence Berkeley National Laboratory National Energy Research Scientific Computing Center's (NERSC) `Edison' supercomputer using 512 processors running for about 5 hours. 

In order to remove the effect of geometric attenuation, and focus the work on the source-to-site geometry effect, all motions were scaled to obtain a $1 \cmts$ maximum displacement in the horizontal component, other components being scaled accordingly. This allows a fair comparison of the relative importance of the presence of 3-D effects (such as surface waves), at least in the linear analysis case, for a given stress drop. The scaling and overall treatment of the fault rupture specification are not meant to be representative of a specific real rupture, rather its meant to emphasize the point that differences between a fully 3-D analysis and equivalent 1-D or \ttod analyses arise from the 3-D nature of problem even in the most basic circumstances. 


% Two faulting styles are explored: dip-slip (with a $45\dgr$ dip angle) and strike slip. 

Figure \ref{fig:input_motions} shows the displacement and accelerations recorded at the NPP site for the different source depths after scaling. \proposedchange{}{Motions occurring at frequencies above $f_{\max}$, unresolved frequencies, appear as high-frequency noise and were filtered out by a forward-backward fourth order Butterworth filter before application to the local model.}   Note the enhanced presence of surface waves for the shallow-source cases, which manifests as a delayed wave (after $t=2.0\seg$). 

Figure \ref{fig:input_motions_fft} shows the absolute value of the discrete Fourier transforms  for the same records shown in Figure \ref{fig:input_motions}. Note the high frequency decay rate of the spectrum, starting at the corner frequency, of two orders of magnitude per decade $\omega^{-2}$, characteristic of the Brune source time function. 


\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./img2/legend.pdf}

\includegraphics[width=0.49\textwidth]{./img2/eq_disp_dip.pdf}
\includegraphics[width=0.49\textwidth]{./img2/eq_acce_dip.pdf}

% \includegraphics[width=0.49\textwidth]{./img2/eq_disp_strike.pdf}
% \includegraphics[width=0.49\textwidth]{./img2/eq_acce_strike.pdf}

\caption{Synthetic records produced by point sources at varying depths, recorded at the site of the NPP. Records are scaled to have $1 \cmts$  of maximum horizontal displacement. }
\label{fig:input_motions}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./img2/legend.pdf}


\includegraphics[width=0.49\textwidth]{./img2/eq_fft_X.pdf}
\includegraphics[width=0.49\textwidth]{./img2/eq_fft_Z.pdf}
% \includegraphics[width=0.49\textwidth]{./img2/eq_fft_Y.pdf}

\caption{Discrete Fourier transforms of the \proposedchange{}{(unfiltered)} displacement records shown in Figure~\ref{fig:input_motions}.}
\label{fig:input_motions_fft}
\end{figure}

% Motions due to a `Full 3-D' wave field, eg. the results from the 3-D seismic simulations are used 

From the simulated 3-D seismic wave field, ie. the `free-field' surface motions computed at the NPP site using the seismic simulation,  vertically incident plane wave-field approximations are developed by de-convolving each component of motion independently. De-convolutions are carried out with a customized version of the SHAKE91\cite{SHAKE91} program for 1-D site response analysis (see considerations by Mejía and Dawson \cite{Mejia2006}). SHAKE91 is a 1-D equivalent linear site response analysis program, as such it only requires a single wave speed (or modulus) and density to specify propagation properties. For horizontal components the SHAKE91 output is interpreted as a vertically incident horizontally-polarized pure shear wave-field, therefore the shear-wave speed and density at the site is used. For the vertical component, SHAKE91 output is interpreted as a laterally-constrained vertically-propagating compression waves, which requires that the speed of 
constrained compression waves and density be used as main parameters.  When only the horizontal component is used the resulting wave-field is called, herein, a `1-D' wave. When all three components of de-convolved motion are input simultaneously the resulting wave-field is called `\ttod' wave since 3 components are used, but they only depend on the depth and time. Finally, using the complete seismic wave is termed `Full 3-D' wave-field. 

\proposedchange{}{Because the ground motions were generated using a smooth source-time function, a simple layered crust with heterogeneous isotropic elastic material and no damping, and because care was taken to filter-out unresolved frequencies, the ground motions at the surface will be highly coherent when time-lag is removed. This is important in the interpretation of the results, where the only source of difference is the 3-D character of the `Full 3-D' wave-field compared to the plane-wave approximation of this motion.}

\section{NPP Modeling}


Figure \ref{fig:05_npp_schematic} shows a drawing of a generic nuclear power plant used as an example structure in this work. Shown are plan and elevation views featuring both the containment as well as the auxiliary building. The structure has a 1:1 in-plan aspect ratio. Post-tensioned and reinforced concrete walls with thicknesses ranging from $0.4\mts$ to $1.6\mts$ are spaced every $12.5\mts$ within the auxiliary building. The containment building consists of a $40\mts$ diameter and $40\mts$ high cylinder capped by a dome, with $1.6\mts$ of wall thickness. The NPP  has a $3.5\mts$ thick continuous foundation slab. 

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{./img2/npp/Model_NPP_02_example_032.pdf}
    \caption{Nuclear power plant geometry}
    \label{fig:05_npp_schematic}
\end{figure}

% \section{Finite-element mesh of the NPP and neighboring soil}

The \texttt{gmsh} \cite{Geuzaine2009} meshing program was used to construct a three-dimensional geometrical model of the NPP and produce a mesh composed of hexahedral finite elements for soil and foundation slab and quadrilateral shell elements for the walls and slabs of the auxiliary building as well as the containment building. 


Figure \ref{fig:05_mesh} shows an overall and an internal view of the generated mesh. The model uses shell elements for containment and auxiliary building and brick elements for foundation slab and soil. The soil volume is separated in an internal domain, a DRM boundary (one layer), and an absorption layer of elements outside of the DRM boundary, which is used to damp any out-going wave motion. The absolution layer is important, because it captures the waves radiated outward from the domain due to the motion of the NPP. These waves are not canceled by DRM because the input motions used to derive DRM forces only include free-field motions. \proposedchange{}{The DRM layer is placed at a distance of $40 \mts $ away from the NPP edges. Therefore, the internal soil box has dimensions $180 \times 180 \times 40$ meters. }

The discretization of the model is controlled by the discretization of the free-field, which is controlled, in turn, by the desired accuracy of wave propagation. Accurately resolving the wave amplitudes for a given maximum frequency imposes restrictions on the spatio-temporal discretization which are more stringent than just meeting the Nyquist criterion\cite{lysmer1969}. For finite-elements with linear displacement interpolation functions involving linear material and geometrical response, it has been shown many times (see Watanabe et. al. \cite{Watanabe2016} for a recent discussion) that at least ten elements per shortest wavelength are required to accurately resolve displacement amplitudes. In addition to that, the time step is limited by the Courant--Friedrichs--Lewy condition\cite{courant1956}, e.g. the time-step must be chosen such that the distance traveled by the fastest phase is less than the smallest spacing between grid nodes for all points in the mesh.  In the present case, considering a target maximum resolvable frequency of $10\Hz$ and the selected surface shear wave velocity, a discretization of size $\ud x = 5\mts$ is needed. Considering the aspects presented in the previous chapter and other geometrical constraints of the mesh, a mesh size of $\ud x = 3.2 \mts$ is selected for the soil. Within the DRM model, which has a uniformly-spaced mesh, the fastest phase is the P-wave with $V_p = 1000\mtsps$, therefore the time step is chosen to meet $\ud t < \ud x / V_p   = 3.2 \mts / 1000\mtsps = 0.0032 \seg$. In practice the time-step of the SW4 output ($0.00083\seg$) is used. 

The resulting mesh consists of $106,165$ nodes, of which $9,728$ define 6 degrees-of-freedom (three translations and three rotations) and $96,437$ define 3 degrees-of-freedom. A total of $100,736$ elements are defined, consisting of $31,824$ shells and ($68,912$) eight-node bricks. $37,632$ of the continuum elements will contain material representing the soil, while the rest use materials representing the DRM layer, absorption layer and foundation slab. 

It is important to mention that, in the finite element model, all shear walls are extended into the foundation slab by adding additional embedded elements within it. This is important for the correct determination of out-of-plane stiffness of the wall/foundation system which would, otherwise, present an unrealistic hinged condition at the base in that direction. 

The high-performance linear quadrilateral ANDES shell finite-element formulation (\cite{Bergan85,Alvin92,Felippa92a,Felippa92b}) implemented in Real ESSI is used herein . This element formulation is able to capture the vibrational dynamics of shells with high accuracy at the cost of forfeiting the concept of displacement interpolation functions. Given the thickness in shear walls for the auxiliary and containment buildings, and because these are often built using post-tensioned concrete, a linear model for concrete behavior suffices. 

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{./img2/npp/fem_mesh.png}
    % \includegraphics[width=0.49\textwidth]{./img2/npp/fig1.png}
     % \includegraphics[width=0.49\textwidth]{./img2/npp/fig3.png}
    % \includegraphics[width=0.49\textwidth]{./img2/npp/fig2.png}
    % \includegraphics[width=0.49\textwidth]{./img2/npp/fig4.png}
    \caption{FEM mesh of the site and NPP featuring, in orange, the layer of DRM elements. }
    \label{fig:05_mesh}
\end{figure}


\begin{figure}[tb]
    \centering
     \includegraphics[width=0.99\textwidth]{./img2/npp/auxilliary_modes.pdf}
     \includegraphics[width=0.99\textwidth]{./img2/npp/containment_modes.pdf}
    % \includegraphics[width=0.49\textwidth]{./img2/npp/fig2.png}
    % \includegraphics[width=0.49\textwidth]{./img2/npp/fig4.png}
    \caption{ Undeformed shape and first six fixed-base eigenfrequencies and eigenmodes for auxiliary building \textbf{(top row)} and containment building \textbf{(bottom row)}.  }
    \label{fig:eigenmodes}
\end{figure}




A Rayleigh damping model is used to assign viscous damping to all parts of the domain. The damping ratio is set at $20\%$ of critical for elements located at the DRM absorption layer, \proposedchange{}{which is meant to model the radiation damping of the soil-structure system. For the soil elements the damping is set to $5\%$ with the idea that the soil will dissipate some additional energy due to localized hysteretic behavior.  Finally, all concrete elements were assumed to remain  linear and show reduced damping (cracking of concrete is avoided in NPP design), therefore a  $2\%$ critical damping ratio was used for structural components and foundations.}

\proposedchange{}{
The only physical role attributed to damping in this work is to absorb the out-going waves from the portion of motions not canceled by the DRM layer. Damping elsewhere in the model plays a cosmetic role, because damping is expected to occur. Indeed the assignment of different damping values to different parts of the domain results in a non-classical damping matrix, which makes the assessment and tuning of a global effective damping for a model of this size a hard task. Furthermore, the use of a mass-proportional damping term in the Rayleigh model is only justifiable in the case of the DRM absorption layer (similar to Lysmer dashpots) and, possibly, within the soil were it to be below the phreatic level (as soils in NPP sites usually are). No claim is made herein that damping terms other than those associated with the DRM layer play a physical role. 
 }


\begin{figure}
\centering
\includegraphics[width=0.32\textwidth]{./img2/npp/positions_top_of_foundation.pdf}
\includegraphics[width=0.32\textwidth]{./img2/npp/positions_top_of_containment.pdf}
\includegraphics[width=0.32\textwidth]{./img2/npp/positions_nw_aux_bldg.pdf}
\caption{Positions on the NPP FEM model where the response is reported.}
\label{fig:positions}
\end{figure}



\proposedchange{}{For reference, Figure~\ref{fig:eigenmodes} shows the first five fixed-base eigenmodes and eigen frequencies for both the auxiliary building and the containment building. The first two translational eigen-modes of the auxiliary have a frequency of $9.6\Hz$, while the same modes for the containment building are at $6.1 \Hz$. This result is within the expected range for both buildings. Of course, the dynamics of the response of the NPP is governed not by the fixed-base modes but by the natural frequencies of the NPP including foundation and site, which have lower frequencies than the fixed-based modes, as will be shown. }

After post-processing the results from the SW4 simulations, nodal displacements were used to obtain nodal accelerations. Data and metadata necessary for DRM computations was stored in an HDF5 dataset designed for high-performance read operations for Real ESSI. Simulations for the nuclear power plant and site were executed using DRM forces computed from these datasets. The simulations were run on 192 processors of the NERSC `Edison' cluster for approximately 4.5 hours. 

\section{Results}



Figure~\ref{fig:positions} shows the locations on the finite-element model for the nuclear power plant where displacement responses were calculated, namely at the bottom of foundation, top of containment building and north-west corner of the auxiliary building. Figures~\ref{fig_z550_3dvs1d} through \ref{fig_z1200_3dvs1d} report results comparing full 3-D analysis with equivalent 1-D analysis, that is, considering only the horizontal component for the 1-D deconvolution of input motions produced by sources at the different depths considered. Next, figures \ref{fig_z550_3dvs3t1d} through \ref{fig_z1200_3dvs3t1d} report the response at the top of foundation and top of the containment building, this time comparing between full 3-D analysis and \ttod analysis, with input motions proceeding from the same source depths as before. For all cases the displacement and acceleration time history responses are shown. 

Figure \ref{fig_3dvs1d_fft} presents the comparison of the Fourier response spectra for the motions recorded at the top of the containment building for sources of varying depths, using full 3-D and 1-D analysis.  Figure \ref{fig_3dvs3t1d_fft} presents the same plots, now comparing with \ttod analysis. 

Finally, figure \ref{fig_cuts_comparison} shows deformed shapes of NPP and site at an instant in time due to a $z_s = 550m$ fault using the different types of input motions, allowing to assess the overall differences in response throughout the building. 


% \clearpage\centering

% \subsection{ $z_s= 1200m$ Full 3-D vs. 1-D}



\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_9/node_73_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_9/node_73_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_9/node_145_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_9/node_145_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_9/node_733_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_9/node_733_disp.pdf}

% \includegraphics[width=0.99\textwidth]{./results/3d_vs_1d_9/node_733_fft.pdf}
\caption{Response of the NPP to the motions produced by a source at $z_s = 550m$ comparing Full 3-D analysis with equivalent 1-D analysis. }
\label{fig_z550_3dvs1d}
\end{figure}



\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_6/node_73_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_6/node_73_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_6/node_145_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_6/node_145_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_6/node_733_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_6/node_733_disp.pdf}

% \includegraphics[width=0.99\textwidth]{./results/3d_vs_1d_6/node_733_fft.pdf}
\caption{Response of the NPP to the motions produced by a source at $z_s = 850m$ comparing Full 3-D analysis with equivalent 1-D analysis. }
\label{fig_z850_3dvs1d}
\end{figure}







\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_3/node_73_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_3/node_73_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_3/node_145_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_3/node_145_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_3/node_733_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3d_vs_1d_3/node_733_disp.pdf}

% \includegraphics[width=0.99\textwidth]{./results/3d_vs_1d_3/node_733_fft.pdf}
\caption{Response of the NPP to the motions produced by a source at $z_s = 1200m$ comparing Full 3-D analysis with equivalent 1-D analysis. }
\label{fig_z1200_3dvs1d}
\end{figure}



\begin{figure}
\centering
\includegraphics[width=0.99\textwidth]{./results/3d_vs_1d_9/node_733_fft.pdf}
\includegraphics[width=0.99\textwidth]{./results/3d_vs_1d_6/node_733_fft.pdf}
\includegraphics[width=0.99\textwidth]{./results/3d_vs_1d_3/node_733_fft.pdf}
\caption{Fourier response of the NPP to the motions produced by a source at $z_s = 550m$ (top), $z_s = 850m$ (middle), $z_s = 1200m$ (bottom) comparing Full 3-D analysis with equivalent 1-D analysis. }
\label{fig_3dvs1d_fft}
\end{figure}




% \clearpage

% \subsection{ $z_s= 550m$ Full 3-D vs. \ttod}



\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./results/9/node_73_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/9/node_73_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/9/node_145_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/9/node_145_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/9/node_733_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/9/node_733_disp.pdf}

% \includegraphics[width=0.99\textwidth]{./results/9/node_733_fft.pdf}
\caption{Response of the NPP to the motions produced by a source at $z_s = 550$ comparing Full 3-D analysis with equivalent \ttod analysis. }
\label{fig_z550_3dvs3t1d}
\end{figure}



% \clearpage
% \subsection{ $z_s= 850m$ Full 3-D vs. \ttod}


\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./results/6/node_73_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/6/node_73_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/6/node_145_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/6/node_145_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/6/node_733_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/6/node_733_disp.pdf}

% \includegraphics[width=0.99\textwidth]{./results/6/node_733_fft.pdf}
\caption{Response of the NPP to the motions produced by a source at $z_s = 850$ comparing Full 3-D analysis with equivalent \ttod analysis. }
\label{fig_z850_3dvs3t1d}
\end{figure}





% \clearpage

% \subsection{ $z_s= 1200m$ Full 3-D vs. \ttod}


\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{./results/3/node_73_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3/node_73_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3/node_145_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3/node_145_disp.pdf}

\includegraphics[width=0.49\textwidth]{./results/3/node_733_acce.pdf}
\includegraphics[width=0.49\textwidth]{./results/3/node_733_disp.pdf}

% \includegraphics[width=0.99\textwidth]{./results/3/node_733_fft.pdf}
\caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
\label{fig_z1200_3dvs3t1d}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=0.99\textwidth]{./results/9/node_733_fft.pdf}
\includegraphics[width=0.99\textwidth]{./results/6/node_733_fft.pdf}
\includegraphics[width=0.99\textwidth]{./results/3/node_733_fft.pdf}
\caption{Fourier response of the NPP to the motions produced by a source at $z_s = 550m$ (top), $z_s = 850m$ (middle), $z_s = 1200m$ (bottom) comparing Full 3-D analysis with equivalent \ttod analysis. }
\label{fig_3dvs3t1d_fft}
\end{figure}



\begin{figure}

\centering

\includegraphics[width=0.70\textwidth]{./img_cross_sections/cut_case9_3d.png}

\includegraphics[width=0.70\textwidth]{./img_cross_sections/cut_case9_1d.png}

\includegraphics[width=0.70\textwidth]{./img_cross_sections/cut_case9_3x1d.png}

\caption{Deformed shape of the NPP and site (exaggerated) at time $t=1.61\seg$ for source at $z_s = 550m$. Top shows full 3-D analysis, while middle and bottom show equivalent 1-D  and \ttod analyses.  Note that the deformation at the far edges (free-field) are nearly identical. }

\label{fig_cuts_comparison}

\end{figure}


\section{Discussion}

A simple calculation using Snell's law,\proposedchange{}{which relates the ratio of the wave speeds with the ratio of the sines of the incidence angles ($V_{s1}/V_{s2} = \sin \theta_1 / \sin \theta_2$)}, shows that for a 45\dgr dipping fault, which is a source of S-waves inclined at $45\dgr$, the incidence angle measured from the vertical at the surface would be $\theta = 10.6\dgr$, $8.9\dgr$ and $8.1\dgr$ for the faults located at $z_s=550\mts$, $850\mts$, and $1200\mts$ respectively. Therefore, the wave-field has an approximately vertical angle of incidence at the site. Using a typical argument in SSI research and practice, even for this simplified numerical experiment, the wave-field is an almost vertically propagating S-wave. This would be then used to argue in favor of using a 1-D or \ttod approach to modeling the seismic wave-field. In this study P-wave speeds and S-wave speeds are in a 2:1 ratio, so incidence angles are affected in the same way for both phases. In general, different wave phases will reflect and diffract depending on their particular speed ratio. 
% \begin{align}
% \dfrac{\sin \theta_1 }{\sin \theta_2 } &= \dfrac{V_{s1}}{V_{s2}} 
% \label{snell}
% \end{align}
From the results presented in the previous section, it can be readily inferred that the responses of the nuclear power plant to a full 3-D, \ttod and 1-D wave-field are different, even if the latter two were developed using a state-of-practice approach of de-convolving the motion in depth. Notably, the deconvolution method matches perfectly the recorded motion at the site of the NPP when only free field motions are considered. This result reduces the credibility of the argument that a near-vertical incidence angle, computed via Snell's law, suffices to establish the adequacy of a 1-D or \ttod SSI analysis for the purposes of establishing the earthquake-soil-structure response of structures.

\proposedchange{}{Comparison with the results by Luco et. al\cite{Luco1980} in the case of close sources shows that wave passage effects are more pronounced (i.e. there seems to be more de-amplification) than the ones obtained by those authors in the case of very near sources. Direct comparison, however, is not possible due to different layering structure and the consideration of a perfectly rigid superstructure in the case of Luco et. al. } 

An overarching theme across all comparisons herein is that wave-fields derived from de-convolution (1-D or \ttod) are associated with increased peak responses for both displacements and accelerations. This is because these wave-fields deliver all the wave energy in phase, at the same time, across the span of the foundation. Conversely, in the full 3-D case the motion starts at one foundation edge and traverses the whole foundation in a phased manner. This is due to the wave-field being inclined as opposed to purely one-dimensional. In the uni-dimensional cases the seismic motions propagate into the NPP structure and constructively interfere inside it. The delayed delivery of energy in the 3-D case means that the motions do not constructively (nor destructively) interfere inside the structure, resulting in reduced net response.  This effect can be simplistically identified as the concept of kinematic interaction and is one of the rationale for the base-averaging effect (the other being motion incoherence). In this particular case, responses are reduced by about 50\% from those computed with de-convolved (1-D or \ttod) motions. The authors would like to stress that this is not a generalizable result.  

Nevertheless, this aspect of response is very much a frequency-dependent one. Indeed, looking at discreet Fourier spectra across all examples (figures \ref{fig_3dvs1d_fft} and \ref{fig_3dvs3t1d_fft}) it is apparent that the 3-D, 1-D and \ttod methods are coincident for low-frequency input except for the deep source where the 1-D method produces different results. Low frequencies are associated with long wavelengths. For example, for a frequency of $0.5\Hz$ a corresponding S-wave will have around $1\kmts$ of length, which exceeds the size of the NPP foundation slab. Therefore, for these frequencies the site motion approximates that of a rigid body and it can be said that a 1-D wave-field will be a good approximation. This is also in-line with the results reported in reference \cite{Luco1980} where it was noted that the wave-passage effect got worse with increasing frequency and decreasing distance to the fault.


At the high frequency end of the Fourier response, it is seen that the 1-D and \ttod wave-fields tend to deliver more energy, although not much can be said for frequencies higher that $5\Hz$ from these results. 

Still with reference to figures \ref{fig_3dvs1d_fft} and \ref{fig_3dvs3t1d_fft}, notable peaks in the Fourier spectra occur at around $0.6 \Hz$, $1.3 \Hz$, and about $3\Hz$. The first two are identified as source-related, because they can be identified in the free-field Fourier responses. The last peak is inferred to be the structure-site natural mode of vibration as measured at the top of the containment building, \proposedchange{}{implying that the fundamental frequency of the containment-building is halved when SSI is included}. Frequency-domain response is also over-estimated by both 1-D and \ttod input modes, with a strong dependence on source depth. 

It is interesting to note the counter-intuitive effect that omitting the vertical component (1-D case) might actually increase the response in some places of the structure. Compare, for example, Figure \ref{fig_z1200_3dvs1d} with Figure \ref{fig_z1200_3dvs3t1d}, the peak horizontal response in terms of both displacements and accelerations at the base of the foundation is in fact increased by removing the vertical component. Looking at the discrete Fourier spectrum for the horizontal response in this case shows that even lower-frequency components are affected in this case. 

Due to the symmetry of this example, response at places like the foundation bottom and top-of containment building show a de-coupling of vertical and horizontal components, as can be seen as a flat response in the response at these points seen in figures \ref{fig_z550_3dvs1d}, \ref{fig_z850_3dvs1d}, and \ref{fig_z1200_3dvs1d}.  This is not the general case and is why looking at a different location, like the north-west corner of the auxiliary building, helps in understanding how vertical and horizontal components might interact. 

Also noteworthy is the fact that vertical components of structural response for 1-D and \ttod input seem to match the response due to a full 3-D wave better than the horizontal components. This phenomenon can be partially explained by the anisotropy of effective shear wave speed present in this example. At $2.5\Hz$ the implied wavelength for the top layer is $200\mts$, but that layer is only $50 \mts$ deep. Conceptually, this means that a wave at this frequency traveling along the layer, horizontally, can develop this wavelength but a vertically propagating wave will have a longer effective wavelength because it will mobilize deeper layers with higher wave speeds. Therefore, it is expected that this is a frequency-dependent phenomenon and the results would change at higher frequencies.

In the cross-sections shown in Figure \ref{fig_cuts_comparison},  showing displaced shape (heavily exaggerated), it is notable how the response at the DRM boundary is nearly identical in all three cases shown (3-D, 1-D and \ttod), except in that the 1-D case does not include the vertical component by definition. 


% \begin{itemize}
% \item Results are not the same. 
% \item 3$\times$1-D analysis tends to deliver increased energy due to plane wave (arrives all at once).
% \item Greater peak responses for 3$\times$1-D 
% \item Identical responses for long periods (long wavelengths are essentially good for plane wave approx. site moves like rigid body at these wavelengths.)
% \item Mid frequencies show increased response for 3$\times$1-D.
% \item Higher frequencies are unclear. Might be very mesh dependent. 
% \item Vertical component is overall similar for both cases. Probably not true for embedded structures. (Base averaging occurs for horizontal components but not vertical!)
% \end{itemize}


% From Figure XXX through YYY it can be concluded that, first, results from full 3-D analysis differ from those that assume a 1-D seismic wave field (\ttod). For the chosen source time-function it is seen that full 3-D analysis yields reduced peak values when compared with \ttod. The reason for this is that in \ttod all the seismic energy arrives at the base of the NPP at the same time, delivering a single coherent pulse into the NPP. Conversely, since the waves arrive at an angle when doing full 3-D, energy is delivered not at the same time but with delays. This lends strength to the idea of `base averaging' effect so commonly cited, but it is cautioned here that this is a frequency-dependent effect and case-by-case assessment with more realistic motions is necessary if this wants to be used to reduce seismic demands. 

% Frequency-domain response results from Figures WWW through ZZZ show some 


\section{Concluding remarks}

The goal of earthquake soil-structure interaction (ESSI) modeling is to explicitly account for all processes which contribute to the seismic response of a soil-structure system, from source fault rupture, through crustal wave propagation,  to soil-structural response and, finally, components. The domain reduction method, by Bielak at. al., is the definitive tool that enables this process. Before the DRM, assumptions like the vertically incident plane-wave approximation were necessary to arrive at an answer, but introduced modeling uncertainties which could not be evaluated because they could not be separated from other sources of uncertainty. 

% The domain reduction method is an effective tool for modeling and simulation of earthquake-soil-structure interaction problems, providing opportunities to model all the complexities of the problem along the way. Using this method proves to be a fully rational framework for the assessment of the safety of infrastructure in the sense that no assumptions are needed to simplify the input of earthquake motions into the model, which can be incorporated in all their complexity. Earthquake scenarios can be postulated and evaluated specifically for a particular structure. This is crucial in evaluating different modeling assumptions especially when data is scarce and of questionable quality. 


This work presents a step towards evaluating the modeling uncertainty introduced by the plane-wave approximations in the seismic response of soil-structure systems. Equivalent 1-D and \ttod wave-field approximations were developed from 3-D wave-fields, obtained through numerical simulation, and the response of a soil-structure system to each of these was evaluated. Differences between the modeling approaches are not just in magnitude of the peak response, but in overall time-history shape and frequency content. By comparing the response of the structure to 3-D wave-fields with the response to point-wise equivalent plane-wave approximations of the same motions the relevance of full 3-D modeling is established. \proposedchange{}{This conclusion coincides and extends the assessment of the United States Nuclear Regulatory Commission via document NUREG/CR-3805\cite{CR3805} on characterization of ground motions, which appraises the importance of incidence angle (equivalently termed apparent velocity) and presence of surface waves in deciding whether a more complex approach is needed for ground motion modeling. }

%I would say that in this case we got conservative result using 3D simulations. However, based on sensitivity studies which are not presented in this paper, the results showed high sensitivity to model parameters. Therefore, the level of conservatism or un-conservatism can be  assessed only on a case by case basis.


Despite the fact that the results shown herein seem to imply that assuming 1-D or \ttod wave-fields is conservative with respect to the full 3-D wave-field, it is very important to note that this is not a generalizable conclusion. Indeed, the results are very sensitive to all model parameters: structure, site and input motions. Therefore, the level of conservatism or un-conservatism can be  assessed only on a case by case basis with a full 3-D analysis. \proposedchange{Due to this, the authors argue that performing ESSI analysis, given the right conditions, is a more rational and complete approach in that it involves less assumptions.}{The role of regional heterogeneity, basin effects, non-linear site response, sloping site or presence of topography, among other complicating factors, in this context, has yet to be quantified. Due to the factibility of modeling these and other effects within the DRM framework without incurring in additional assumptions, the authors argue that performing ESSI analysis, given the right conditions, is a more rational and complete approach}. 



A maximum resolved frequency of $10\Hz$ was achieved in this study, which is insufficient to properly characterize all important aspects of NPP response pertaining to safety of its contents, but is enough to provide evidence that assuming a 1-D or \ttod wave comes with significant changes in response amplitude, shape and frequency content. Achieving higher frequencies is an important goal in order to completely assess hazard to NPPs due to nearby faults. Modeling for higher frequencies, becomes a significantly harder task both computationally and from a modeling standpoint. First, computational cost grows cubicly with maximum frequency such that achieving, say $30\Hz$, is possible for a small near-field model like the one in this article if a larger HPC clusters is available. Regional-scale models are out of the question at this moment for these kinds of frequencies, though programs like SW4 are getting ready for the next generation of exascale machines\cite{peterssen2017} that will tackle this problem. Second, and very importantly, modeling higher frequencies requires an adequate frequency-dependent representation of both anelastic attenuation (inherent damping) and scattering effects.  Future work will look at extended sources, high frequency motions and the role of non-linearity, but the intent of this paper was to develop fundamental understanding of the mechanics and which can be built upon to include these and many other features.

Realistic seismic modeling of the input DRM motions requires that complex, extended faults with multiple sub-fault segments, most commonly represented by linear superposition of point sources, be used (a concrete example on how to setup such a scenario can be found in \cite{Isbiliroglu2015}). In this way it is possible to model important extended-fault effects such as time-varying frequency content, pulse directivity and longer duration motions with more motion-reversal cycles. These effects become important when evaluating non-linear effects, especially in the presence of hysteresis. It is shown here, though, that the response discrepancies between models excited by full 3-D seismic fields and plane-wave approximations of these arise from the very element used to construct extended sources: the point source. Since large sources are compositions of point sources, these 3-D effects will be present for extended sources too. 


% {
% \color{red}
% Need to extend this to now extend this to more complex, realistic scenarios with distributed faults etc.
% }

% The authors promote the idea that evaluation of earthquake safety of infrastructure using explicit modeling of earthquake event, seismic wave propagation, site response and soil-structure interaction based on the DRM is the most rational approach currently available. In particular it leads to results different from those obtained by using 1-D wave-field simplifications. It should be regarded as another tool in the engineer's toolkit to evaluate safety of civil infrastructure. Its place in a modeling hierarchy is that of the most advanced and detailed analysis method currently available, and its modeling/simulation cost reflects this fact. 

For concrete applications, it is important that the earthquake models used to derive DRM input motions should also be able to reproduce, within uncertainties, the most important aspects of the motions observed at the seismic networks within the modeling region. In addition, it important to note that the greatest flexibility of the \proposedchange{framework}{DRM framework applied herein}, by virtue of its completely physical formulation, is that it can also be used to postulate and evaluate new scenarios, which are plausible in the seismic setting but have never been observed and recorded. 

% In the light of this work, a new question arises: if the seismic safety of an infrastructure object is shown based on a plane-wave approximation of the seismic wave-field, is it also safe in the case of general 3-D input?

% Especially in applications involving non-linearity of the hysteric type where there will be a large dependence in the energy dissipated with the duration and frequency content of the simulated records. The influence of the 1-D family of assumptions in this type of models is yet to be determined, but the present work based on simple linear theory presents solid grounds to claim that differences can be expected to arise. 




% Full validation of a DRM-based simulation requires, first, a validated seismic and crustal velocity structure model, an attenuation model for the crust ($Q$ factors), and an instrumented structure 

% \begin{itemize}
%     \item Full 3D analysis is possible. This is how its done. 
%     \item Results are not obvious. Frequency, geometry, site, etc. dependence. 
%     \item It is necessary for important structures that are sensitive to wave passage, surface waves, etc. Classic NPPs, SMRs, bridges, tall structure, tank sloshing. 
%     \item Can lead to increase/decrease in seismic demands or performance. Unpredictable.
%     \item 3$\times$1-D analysis tends to deliver increased energy. Full 3-D might be the way to prove that demands are actually lower. 
%     \item Difficulty of these analysis. Requisites for generating good DRM motions. 
%     \item Importance of good (validated) seismic rupture models. 
%     \item Importance of 3-D crust models (increased durations, scattering, etc.)
%     \item Importance of validation of non-linear 3D site response. 
%     \item Differences should be greater for non-linear analysis (over-predicting plastic deformations) but need more realistic source models (more cycles, extended fault effects)

% \end{itemize}


\ack 

The work presented in this paper was partially supported by the Canadian Nuclear Safety Commission (CNSC) and the United States Department of Energy (US DOE). This support is gratefully acknowledged.

\bibliographystyle{wileyj}

\bibliography{./_MyPapers-2016-NearFiieldESSINPPs} % The file containing the bibliography


% \section*{Appendix}

% % /home/jaabell/Papers/Journals/2016/NearFieldESSINPPs/results_old/2
% % /home/jaabell/Papers/Journals/2016/NearFieldESSINPPs/results_old/4
% % /home/jaabell/Papers/Journals/2016/NearFieldESSINPPs/results_old/5
% % /home/jaabell/Papers/Journals/2016/NearFieldESSINPPs/results_old/7
% % /home/jaabell/Papers/Journals/2016/NearFieldESSINPPs/results_old/10
% % /home/jaabell/Papers/Journals/2016/NearFieldESSINPPs/results_old/12
% Strike slip faults. 


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

% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/2/node_73_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/2/node_73_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/2/node_73_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_1}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/2/node_145_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/2/node_145_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/2/node_145_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_2}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/2/node_733_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/2/node_733_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/2/node_733_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_3}
% \end{figure}




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

% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/4/node_73_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/4/node_73_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/4/node_73_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_1}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/4/node_145_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/4/node_145_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/4/node_145_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_2}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/4/node_733_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/4/node_733_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/4/node_733_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_3}
% \end{figure}



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

% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/5/node_73_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/5/node_73_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/5/node_73_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_1}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/5/node_145_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/5/node_145_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/5/node_145_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_2}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/5/node_733_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/5/node_733_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/5/node_733_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_3}
% \end{figure}




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

% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/7/node_73_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/7/node_73_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/7/node_73_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_1}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/7/node_145_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/7/node_145_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/7/node_145_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_2}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/7/node_733_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/7/node_733_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/7/node_733_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_3}
% \end{figure}


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

% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/10/node_73_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/10/node_73_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/10/node_73_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_1}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/10/node_145_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/10/node_145_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/10/node_145_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_2}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/10/node_733_acce.pdf}
% \includegraphics[width=0.49\textwidth]{./results_old/10/node_733_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/10/node_733_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_3}
% \end{figure}



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

% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/12/node_73_acce.pdf}12\includegraphics[width=0.49\textwidth]{./results_old/2/node_73_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/12/node_73_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_1}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/12/node_145_acce.pdf}12\includegraphics[width=0.49\textwidth]{./results_old/2/node_145_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/12/node_145_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_2}
% \end{figure}


% \begin{figure}
% \centering
% \includegraphics[width=0.49\textwidth]{./results_old/12/node_733_acce.pdf}12\includegraphics[width=0.49\textwidth]{./results_old/2/node_733_disp.pdf}

% % \includegraphics[width=0.5\textwidth]{./results_old/12/node_733_fft.pdf}
% \caption{Response of the NPP to the motions produced by a source at $z_s = 1200$ comparing Full 3-D analysis with equivalent \ttod analysis. }
% \label{fig_z1200_3dvs3t1d_3}
% \end{figure}



\end{document}